Moving a body that a boolean made no longer hangs: a boolean between a cylinder and a tilted plane (an elliptical edge) froze the app at 100 % CPU, because its retries move both bodies and the result back
Moving a body transforms its intersection-curve edges, and truck-geometry 0.5.0 scales each curve's tolerance by the operator norm of the matrix. matext4cgmath 0.1.0 computes it from eigenvalues through a quartic whose Newton refinement loops while the residual is above a tolerance and never gives up. For a translation matrix all four eigenvalues are 1, and near that quadruple root the iteration can cycle forever (about one translation in eight of the boolean retries' shifts, by a scan of sizes 0.1-2000 mm). matext4cgmath is vendored with both Newton loops stopped after 100 steps (vendor/README.md, ATTRIBUTION.md). truck master no longer computes the norm in transform_by, so the vendored copy can go once a truck release after 0.5.0 is adopted. Tests: moving_a_body_with_an_intersection_curve_edge_ends (kernel) cuts a cylinder with a tilted box and moves the result by one of the shifts that hung (hangs on main). examples/regressions/oblique_cut_of_cylinder_hang.json is the variant of the oblique cut that hung the app; it hangs on main, and since the previous commit it passes without reaching the shifted retries, so it is replayed by the engine test as a scenario only. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
This commit is contained in:
CypherNaugh_0xandClaude Opus 5.5 authored and bherbruck committed
1 parent
f3a689f284
commit
5249ae84ef
14 files changed
+920
-3
No files matched your search
Vendored
+153
@@ -0,0 +1,153 @@
|
||||
use crate::*;
|
||||
|
||||
impl<F: BaseFloat> EigenValues for Matrix2<F> {
|
||||
type EigenValues = [Complex<F>; 2];
|
||||
fn eigenvalues(self) -> [Complex<F>; 2] {
|
||||
solver::solve_quadratic(-self.trace(), self.determinant())
|
||||
}
|
||||
}
|
||||
|
||||
impl<F: BaseFloat> EigenValues for Matrix3<F> {
|
||||
type EigenValues = [Complex<F>; 3];
|
||||
fn eigenvalues(self) -> [Complex<F>; 3] {
|
||||
let poly2 = self[0][0] * self[1][1] + self[0][0] * self[2][2] + self[1][1] * self[2][2]
|
||||
- self[0][1] * self[1][0]
|
||||
- self[1][2] * self[2][1]
|
||||
- self[0][2] * self[2][0];
|
||||
solver::solve_cubic(-self.trace(), poly2, -self.determinant())
|
||||
}
|
||||
}
|
||||
|
||||
impl<F: BaseFloat> EigenValues for Matrix4<F> {
|
||||
type EigenValues = [Complex<F>; 4];
|
||||
fn eigenvalues(self) -> [Complex<F>; 4] {
|
||||
let poly2 = self[0][0] * self[1][1]
|
||||
+ self[0][0] * self[2][2]
|
||||
+ self[0][0] * self[3][3]
|
||||
+ self[1][1] * self[2][2]
|
||||
+ self[1][1] * self[3][3]
|
||||
+ self[2][2] * self[3][3]
|
||||
- self[0][1] * self[1][0]
|
||||
- self[0][2] * self[2][0]
|
||||
- self[0][3] * self[3][0]
|
||||
- self[1][2] * self[2][1]
|
||||
- self[1][3] * self[3][1]
|
||||
- self[2][3] * self[3][2];
|
||||
let poly3 = (self[0][0] * self[1][1] + self[0][0] * self[2][2] + self[1][1] * self[2][2]
|
||||
- self[0][1] * self[1][0]
|
||||
- self[0][2] * self[2][0]
|
||||
- self[1][2] * self[2][1])
|
||||
* self[3][3]
|
||||
+ (self[2][0] * self[3][2] - self[2][2] * self[3][0] + self[1][0] * self[3][1]
|
||||
- self[1][1] * self[3][0])
|
||||
* self[0][3]
|
||||
+ (self[2][1] * self[3][2] - self[2][2] * self[3][1] + self[3][0] * self[0][1]
|
||||
- self[3][1] * self[0][0])
|
||||
* self[1][3]
|
||||
+ (self[3][1] * self[1][2] - self[3][2] * self[1][1] + self[3][0] * self[0][2]
|
||||
- self[3][2] * self[0][0])
|
||||
* self[2][3]
|
||||
+ Matrix3::from_cols(self[0].truncate(), self[1].truncate(), self[2].truncate())
|
||||
.determinant();
|
||||
solver::solve_quartic(-self.trace(), poly2, -poly3, self.determinant())
|
||||
}
|
||||
}
|
||||
|
||||
impl<F: BaseFloat> OperatorNorm for Matrix2<F> {
|
||||
#[inline]
|
||||
fn norm_l1(self) -> F {
|
||||
F::max(
|
||||
F::abs(self[0][0]) + F::abs(self[0][1]),
|
||||
F::abs(self[1][0]) + F::abs(self[1][1]),
|
||||
)
|
||||
}
|
||||
#[inline]
|
||||
fn norm_linf(self) -> F {
|
||||
F::max(
|
||||
F::abs(self[0][0]) + F::abs(self[1][0]),
|
||||
F::abs(self[0][1]) + F::abs(self[1][1]),
|
||||
)
|
||||
}
|
||||
#[inline]
|
||||
fn norm_l2(self) -> F {
|
||||
let eigens = (self.transpose() * self).eigenvalues();
|
||||
F::sqrt(F::max(eigens[0].re, eigens[1].re))
|
||||
}
|
||||
}
|
||||
|
||||
impl<F: BaseFloat> OperatorNorm for Matrix3<F> {
|
||||
#[inline]
|
||||
fn norm_l1(self) -> F {
|
||||
F::max(
|
||||
F::max(
|
||||
F::abs(self[0][0]) + F::abs(self[0][1]) + F::abs(self[0][2]),
|
||||
F::abs(self[1][0]) + F::abs(self[1][1]) + F::abs(self[1][2]),
|
||||
),
|
||||
F::abs(self[2][0]) + F::abs(self[2][1]) + F::abs(self[2][2]),
|
||||
)
|
||||
}
|
||||
#[inline]
|
||||
fn norm_linf(self) -> F {
|
||||
F::max(
|
||||
F::max(
|
||||
F::abs(self[0][0]) + F::abs(self[1][0]) + F::abs(self[2][0]),
|
||||
F::abs(self[0][1]) + F::abs(self[1][1]) + F::abs(self[2][1]),
|
||||
),
|
||||
F::abs(self[0][2]) + F::abs(self[1][2]) + F::abs(self[2][2]),
|
||||
)
|
||||
}
|
||||
#[inline]
|
||||
fn norm_l2(self) -> F {
|
||||
let eigens = (self.transpose() * self).eigenvalues();
|
||||
F::sqrt(F::max(F::max(eigens[0].re, eigens[1].re), eigens[2].re))
|
||||
}
|
||||
}
|
||||
|
||||
impl<F: BaseFloat> OperatorNorm for Matrix4<F> {
|
||||
#[inline]
|
||||
fn norm_l1(self) -> F {
|
||||
F::max(
|
||||
F::max(
|
||||
F::max(
|
||||
F::abs(self[0][0])
|
||||
+ F::abs(self[0][1])
|
||||
+ F::abs(self[0][2])
|
||||
+ F::abs(self[0][3]),
|
||||
F::abs(self[1][0])
|
||||
+ F::abs(self[1][1])
|
||||
+ F::abs(self[1][2])
|
||||
+ F::abs(self[1][3]),
|
||||
),
|
||||
F::abs(self[2][0]) + F::abs(self[2][1]) + F::abs(self[2][2]) + F::abs(self[2][3]),
|
||||
),
|
||||
F::abs(self[3][0]) + F::abs(self[3][1]) + F::abs(self[3][2]) + F::abs(self[3][3]),
|
||||
)
|
||||
}
|
||||
#[inline]
|
||||
fn norm_linf(self) -> F {
|
||||
F::max(
|
||||
F::max(
|
||||
F::max(
|
||||
F::abs(self[0][0])
|
||||
+ F::abs(self[1][0])
|
||||
+ F::abs(self[2][0])
|
||||
+ F::abs(self[3][0]),
|
||||
F::abs(self[0][1])
|
||||
+ F::abs(self[1][1])
|
||||
+ F::abs(self[2][1])
|
||||
+ F::abs(self[3][1]),
|
||||
),
|
||||
F::abs(self[0][2]) + F::abs(self[1][2]) + F::abs(self[2][2]) + F::abs(self[3][2]),
|
||||
),
|
||||
F::abs(self[0][3]) + F::abs(self[1][3]) + F::abs(self[2][3]) + F::abs(self[3][3]),
|
||||
)
|
||||
}
|
||||
#[inline]
|
||||
fn norm_l2(self) -> F {
|
||||
let eigens = (self.transpose() * self).eigenvalues();
|
||||
F::sqrt(F::max(
|
||||
F::max(F::max(eigens[0].re, eigens[1].re), eigens[2].re),
|
||||
eigens[3].re,
|
||||
))
|
||||
}
|
||||
}
|
||||
+106
@@ -0,0 +1,106 @@
|
||||
use crate::*;
|
||||
|
||||
impl<F: BaseFloat> Exponential for Matrix2<F> {}
|
||||
impl<F: BaseFloat> Exponential for Matrix3<F> {}
|
||||
impl<F: BaseFloat> Exponential for Matrix4<F> {}
|
||||
|
||||
impl<F: BaseFloat> Decomposition for Matrix2<F> {
|
||||
fn iwasawa_decomposition(self) -> Option<(Self, Self, Self)> {
|
||||
let v0 = self[0];
|
||||
let mag0 = dot(v0, v0);
|
||||
if mag0 <= F::zero() {
|
||||
return None;
|
||||
}
|
||||
let n0 = dot(self[1], v0) / mag0;
|
||||
let v1 = self[1] - v0 * n0;
|
||||
let mag1 = dot(v1, v1);
|
||||
if mag1 <= F::zero() {
|
||||
return None;
|
||||
}
|
||||
let (a0, a1) = (F::sqrt(mag0), F::sqrt(mag1));
|
||||
Some((
|
||||
Matrix2::from_cols(v0 / a0, v1 / a1),
|
||||
Matrix2::from_diagonal(Vector2::new(a0, a1)),
|
||||
Matrix2::new(F::one(), F::zero(), n0, F::one()),
|
||||
))
|
||||
}
|
||||
}
|
||||
|
||||
impl<F: BaseFloat> Decomposition for Matrix3<F> {
|
||||
fn iwasawa_decomposition(self) -> Option<(Self, Self, Self)> {
|
||||
let v0 = self[0];
|
||||
let mag0 = dot(v0, v0);
|
||||
if mag0 <= F::zero() {
|
||||
return None;
|
||||
}
|
||||
let n0 = dot(self[1], v0) / mag0;
|
||||
let v1 = self[1] - v0 * n0;
|
||||
let mag1 = dot(v1, v1);
|
||||
if mag1 <= F::zero() {
|
||||
return None;
|
||||
}
|
||||
let n1 = dot(self[2], v0) / mag0;
|
||||
let n2 = dot(self[2], v1) / mag1;
|
||||
let v2 = self[2] - v0 * n1 - v1 * n2;
|
||||
let mag2 = dot(v2, v2);
|
||||
if mag2 <= F::zero() {
|
||||
return None;
|
||||
}
|
||||
let (a0, a1, a2) = (F::sqrt(mag0), F::sqrt(mag1), F::sqrt(mag2));
|
||||
#[rustfmt::skip]
|
||||
let res = Some((
|
||||
Matrix3::from_cols(v0 / a0, v1 / a1, v2 / a2),
|
||||
Matrix3::from_diagonal(Vector3::new(a0, a1, a2)),
|
||||
Matrix3::new(
|
||||
F::one(), F::zero(), F::zero(),
|
||||
n0, F::one(), F::zero(),
|
||||
n1, n2, F::one(),
|
||||
),
|
||||
));
|
||||
res
|
||||
}
|
||||
}
|
||||
|
||||
impl<F: BaseFloat> Decomposition for Matrix4<F> {
|
||||
fn iwasawa_decomposition(self) -> Option<(Self, Self, Self)> {
|
||||
let v0 = self[0];
|
||||
let mag0 = dot(v0, v0);
|
||||
if mag0 <= F::zero() {
|
||||
return None;
|
||||
}
|
||||
let n0 = dot(self[1], v0) / mag0;
|
||||
let v1 = self[1] - v0 * n0;
|
||||
let mag1 = dot(v1, v1);
|
||||
if mag1 <= F::zero() {
|
||||
return None;
|
||||
}
|
||||
let n1 = dot(self[2], v0) / mag0;
|
||||
let n2 = dot(self[2], v1) / mag1;
|
||||
let v2 = self[2] - v0 * n1 - v1 * n2;
|
||||
let mag2 = dot(v2, v2);
|
||||
if mag2 <= F::zero() {
|
||||
return None;
|
||||
}
|
||||
let n3 = dot(self[3], v0) / mag0;
|
||||
let n4 = dot(self[3], v1) / mag1;
|
||||
let n5 = dot(self[3], v2) / mag2;
|
||||
let v3 = self[3] - v0 * n3 - v1 * n4 - v2 * n5;
|
||||
let mag3 = dot(v3, v3);
|
||||
if mag3 <= F::zero() {
|
||||
return None;
|
||||
}
|
||||
let (a0, a1, a2, a3) = (F::sqrt(mag0), F::sqrt(mag1), F::sqrt(mag2), F::sqrt(mag3));
|
||||
#[rustfmt::skip]
|
||||
let res = Some((
|
||||
Matrix4::from_cols(v0 / a0, v1 / a1, v2 / a2, v3 / a3),
|
||||
Matrix4::from_diagonal(Vector4::new(a0, a1, a2, a3)),
|
||||
Matrix4::new(
|
||||
F::one(), F::zero(), F::zero(), F::zero(),
|
||||
n0, F::one(), F::zero(), F::zero(),
|
||||
n1, n2, F::one(), F::zero(),
|
||||
n3, n4, n5, F::one(),
|
||||
),
|
||||
));
|
||||
res
|
||||
}
|
||||
}
|
||||
Vendored
+120
@@ -0,0 +1,120 @@
|
||||
//! Unofficial third-party cgmath extensions for calculate eigenvalues, operator norms and Iwasawa decomposition.
|
||||
|
||||
// SolveCraft: compiler warnings no longer fail the build.
|
||||
#![deny(clippy::all, rust_2018_idioms)]
|
||||
#![warn(
|
||||
missing_docs,
|
||||
missing_debug_implementations,
|
||||
trivial_casts,
|
||||
trivial_numeric_casts,
|
||||
unsafe_code,
|
||||
unstable_features,
|
||||
unused_import_braces,
|
||||
unused_qualifications
|
||||
)]
|
||||
#![allow(clippy::many_single_char_names)]
|
||||
|
||||
pub use cgmath;
|
||||
use cgmath::*;
|
||||
use num_complex::Complex;
|
||||
|
||||
mod eigens;
|
||||
mod exp_decomp;
|
||||
/// solvers for low dimensional algebraic equations.
|
||||
pub mod solver;
|
||||
|
||||
/// extension for eigen values
|
||||
pub trait EigenValues: VectorSpace {
|
||||
/// the type of the array of eigen values
|
||||
type EigenValues;
|
||||
/// calculate eigen values.
|
||||
///
|
||||
/// # Examples
|
||||
///
|
||||
/// ```
|
||||
/// use cgmath::*;
|
||||
/// use matext4cgmath::*;
|
||||
/// use num_complex::Complex;
|
||||
/// const EPS: f64 = 1.0e-10;
|
||||
///
|
||||
/// let mat = Matrix2::new(4.0, -2.0, 3.0, -1.0);
|
||||
/// let mut eigens = mat.eigenvalues();
|
||||
/// // Even in the case of real solutions, the order in the array is not guaranteed.
|
||||
/// eigens.sort_by(|x, y| x.re.partial_cmp(&y.re).unwrap());
|
||||
/// assert!(Complex::norm(eigens[0] - 1.0) < EPS);
|
||||
/// assert!(Complex::norm(eigens[1] - 2.0) < EPS);
|
||||
/// ```
|
||||
fn eigenvalues(self) -> Self::EigenValues;
|
||||
}
|
||||
|
||||
#[cfg_attr(doc, katexit::katexit)]
|
||||
/// [operator norms](https://en.wikipedia.org/wiki/Matrix_norm): $L^1$, $L^2$, and $L^\infty$.
|
||||
/// # Examples
|
||||
///
|
||||
/// ```
|
||||
/// use cgmath::*;
|
||||
/// use matext4cgmath::*;
|
||||
///
|
||||
/// let mat = Matrix2::new(1.0, 3.0, -2.0, 4.0);
|
||||
/// // L^1 operator norm is the maximum column absolute sumation.
|
||||
/// assert_eq!(mat.norm_l1(), 6.0);
|
||||
/// // L^2 operator norm is the maximum singular value.
|
||||
/// let ans_norm_l2 = (5.0 + f64::sqrt(5.0)) / f64::sqrt(2.0);
|
||||
/// assert!(f64::abs(mat.norm_l2() - ans_norm_l2) < 1.0e-10);
|
||||
/// // L^∞ operator norm is the maximum row absolute sumation.
|
||||
/// assert_eq!(mat.norm_linf(), 7.0);
|
||||
/// ```
|
||||
pub trait OperatorNorm: VectorSpace {
|
||||
/// operator norm for absolute value sumation: $L^1$.
|
||||
fn norm_l1(self) -> Self::Scalar;
|
||||
/// operator norm for Euclidean norm: $L^2$.
|
||||
fn norm_l2(self) -> Self::Scalar;
|
||||
/// operator norm for absolute value maximum: $L^{\infty}$.
|
||||
fn norm_linf(self) -> Self::Scalar;
|
||||
}
|
||||
|
||||
/// calculate exponential value
|
||||
pub trait Exponential: OperatorNorm + std::ops::AddAssign<Self> + One
|
||||
where
|
||||
Self::Scalar: BaseFloat,
|
||||
{
|
||||
/// calculate exponential
|
||||
///
|
||||
/// # Examples
|
||||
///
|
||||
/// ```
|
||||
/// use cgmath::*;
|
||||
/// use matext4cgmath::*;
|
||||
///
|
||||
/// let x = Matrix2::new(0.0, 1.0, -1.0, 0.0);
|
||||
/// let res = x.exp();
|
||||
/// let ans = Matrix2::from_angle(Rad(1.0));
|
||||
/// assert!((res - ans).norm_l1() < 1.0e-10);
|
||||
/// ```
|
||||
fn exp(self) -> Self {
|
||||
use num_traits::{Float, NumCast};
|
||||
let eps = <Self::Scalar as Float>::epsilon();
|
||||
let mut a = self;
|
||||
let mut res = Self::one();
|
||||
for i in 2_u16..=64 {
|
||||
res += a;
|
||||
a = a * self / <Self::Scalar as NumCast>::from(i).unwrap();
|
||||
if a.norm_linf() < eps * self.norm_linf() {
|
||||
break;
|
||||
}
|
||||
}
|
||||
res
|
||||
}
|
||||
}
|
||||
|
||||
#[cfg_attr(doc, katexit::katexit)]
|
||||
/// some decompositions of matrix
|
||||
pub trait Decomposition: VectorSpace {
|
||||
#[cfg_attr(doc, katexit::katexit)]
|
||||
/// Returns $(K, A, N)$ of the Iwasawa decomposition: $M = KAN$.
|
||||
///
|
||||
/// - $K$: orthonormal matrix
|
||||
/// - $A$: diagonal matrix
|
||||
/// - $N$: upper-half unipotent matrix
|
||||
fn iwasawa_decomposition(self) -> Option<(Self, Self, Self)>;
|
||||
}
|
||||
Vendored
+246
@@ -0,0 +1,246 @@
|
||||
use crate::*;
|
||||
|
||||
#[cfg_attr(doc, katexit::katexit)]
|
||||
/// Solve equation: $x^2 + ax + b = 0$.
|
||||
///
|
||||
/// # Examples
|
||||
///
|
||||
/// ```
|
||||
/// use matext4cgmath::solver;
|
||||
/// use num_complex::Complex;
|
||||
/// const EPS: f64 = 1.0e-10;
|
||||
///
|
||||
/// let res = solver::solve_quadratic(3.0, -4.0);
|
||||
/// // Real solutions are always sorted in order of solution size.
|
||||
/// assert!(Complex::norm(res[0] + 4.0) < EPS);
|
||||
/// assert!(Complex::norm(res[1] - 1.0) < EPS);
|
||||
/// ```
|
||||
pub fn solve_quadratic<F: BaseFloat>(a: F, b: F) -> [Complex<F>; 2] {
|
||||
let two = F::one() + F::one();
|
||||
let four = two + two;
|
||||
let det = a * a - four * b;
|
||||
match det >= F::zero() {
|
||||
true => {
|
||||
let h = F::sqrt(det);
|
||||
[
|
||||
Complex::new((-a - h) / two, F::zero()),
|
||||
Complex::new((-a + h) / two, F::zero()),
|
||||
]
|
||||
}
|
||||
false => {
|
||||
let h = F::sqrt(-det);
|
||||
[Complex::new(-a, h) / two, Complex::new(-a, -h) / two]
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
#[cfg_attr(doc, katexit::katexit)]
|
||||
/// Solve equation: $x^3 + px + q = 0$.
|
||||
///
|
||||
/// # Examples
|
||||
///
|
||||
/// ```
|
||||
/// use matext4cgmath::solver;
|
||||
/// use num_complex::Complex;
|
||||
/// const EPS: f64 = 1.0e-10;
|
||||
///
|
||||
/// let mut res = solver::pre_solve_cubic(-7.0, -6.0);
|
||||
/// // Even in the case of real solutions, the order in the array is not guaranteed.
|
||||
/// res.sort_by(|x, y| x.re.partial_cmp(&y.re).unwrap());
|
||||
/// let ans = [Complex::from(-2.0), Complex::from(-1.0), Complex::from(3.0)];
|
||||
/// res.iter().zip(ans).for_each(|(x, y)| {
|
||||
/// assert!(Complex::norm(x - y) < EPS);
|
||||
/// });
|
||||
/// ```
|
||||
pub fn pre_solve_cubic<F: BaseFloat>(p: F, q: F) -> [Complex<F>; 3] {
|
||||
let two = F::one() + F::one();
|
||||
let three = two + F::one();
|
||||
let sqrt3_2 = F::sqrt(three) / two;
|
||||
let omega = Complex::new(-F::one() / two, sqrt3_2);
|
||||
let omega2 = Complex::new(-F::one() / two, -sqrt3_2);
|
||||
let eps_2 = F::sqrt(F::epsilon());
|
||||
|
||||
let p_3 = p / three;
|
||||
let q_2 = q / two;
|
||||
let alpha2 = q_2 * q_2 + p_3 * p_3 * p_3;
|
||||
let (x, y) = match alpha2 >= F::zero() {
|
||||
true => {
|
||||
let alpha = F::sqrt(alpha2);
|
||||
let tmpx = -q_2 - alpha;
|
||||
let tmpy = -q_2 + alpha;
|
||||
(
|
||||
Complex::new(
|
||||
F::signum(tmpx) * F::powf(F::abs(tmpx), F::one() / three),
|
||||
F::zero(),
|
||||
),
|
||||
Complex::new(
|
||||
F::signum(tmpy) * F::powf(F::abs(tmpy), F::one() / three),
|
||||
F::zero(),
|
||||
),
|
||||
)
|
||||
}
|
||||
false => {
|
||||
let alphai = F::sqrt(-alpha2);
|
||||
(
|
||||
Complex::powf(Complex::new(-q_2, alphai), F::one() / three),
|
||||
Complex::powf(Complex::new(-q_2, -alphai), F::one() / three),
|
||||
)
|
||||
}
|
||||
};
|
||||
let mut res = [x + y, omega * x + omega2 * y, omega2 * x + omega * y];
|
||||
// precision by Newton method
|
||||
res.iter_mut().for_each(|x| {
|
||||
let mut f = *x * *x * *x + *x * p + q;
|
||||
let mut f_prime = *x * *x * three + p;
|
||||
// SolveCraft: at most 100 steps. Near a double root (the eigenvalues of a
|
||||
// translation matrix) the iteration can cycle forever without meeting the tolerance.
|
||||
let mut steps = 0;
|
||||
while f.norm() > eps_2 * f_prime.norm() && steps < 100 {
|
||||
steps += 1;
|
||||
if f_prime.norm() < eps_2 {
|
||||
return;
|
||||
}
|
||||
*x -= f / f_prime;
|
||||
f = *x * *x * *x + *x * p + q;
|
||||
f_prime = *x * *x * three + p;
|
||||
}
|
||||
});
|
||||
res
|
||||
}
|
||||
|
||||
#[cfg_attr(doc, katexit::katexit)]
|
||||
/// solve equation: $x^3 + ax^2 + bx + c = 0$.
|
||||
///
|
||||
/// # Examples
|
||||
///
|
||||
/// ```
|
||||
/// use matext4cgmath::solver;
|
||||
/// use num_complex::Complex;
|
||||
/// const EPS: f64 = 1.0e-10;
|
||||
///
|
||||
/// let mut res = solver::solve_cubic(-3.0, 0.0, 4.0);
|
||||
/// // Even in the case of real solutions, the order in the array is not guaranteed.
|
||||
/// res.sort_by(|x, y| x.re.partial_cmp(&y.re).unwrap());
|
||||
/// let ans = [Complex::from(-1.0), Complex::from(2.0), Complex::from(2.0)];
|
||||
/// res.iter().zip(ans).for_each(|(x, y)| {
|
||||
/// assert!(Complex::norm(x - y) < EPS);
|
||||
/// });
|
||||
/// ```
|
||||
pub fn solve_cubic<F: BaseFloat>(a: F, b: F, c: F) -> [Complex<F>; 3] {
|
||||
let two = F::one() + F::one();
|
||||
let three = two + F::one();
|
||||
let twenty_seven = three * three * three;
|
||||
let p = b - a * a / three;
|
||||
let q = c - a * b / three + two * a * a * a / twenty_seven;
|
||||
let mut res = pre_solve_cubic(p, q);
|
||||
res.iter_mut().for_each(|x| {
|
||||
*x -= a / three;
|
||||
});
|
||||
res
|
||||
}
|
||||
|
||||
#[cfg_attr(doc, katexit::katexit)]
|
||||
/// solve equation: $x^4 + px^2 + qx + r = 0$.
|
||||
///
|
||||
/// # Examples
|
||||
///
|
||||
/// ```
|
||||
/// use matext4cgmath::solver;
|
||||
/// use num_complex::Complex;
|
||||
/// const EPS: f64 = 1e-7;
|
||||
///
|
||||
/// let mut res = solver::pre_solve_quartic(-5.0, 0.0, 4.0);
|
||||
/// // Even in the case of real solutions, the order in the array is not guaranteed.
|
||||
/// res.sort_by(|x, y| x.re.partial_cmp(&y.re).unwrap());
|
||||
/// let ans = [Complex::from(-2.0), Complex::from(-1.0), Complex::from(1.0), Complex::from(2.0)];
|
||||
/// res.iter().zip(ans).for_each(|(x, y)| {
|
||||
/// assert!(Complex::norm(x - y) < EPS);
|
||||
/// });
|
||||
/// ```
|
||||
pub fn pre_solve_quartic<F: BaseFloat>(p: F, q: F, r: F) -> [Complex<F>; 4] {
|
||||
let one = F::one();
|
||||
let two = one + one;
|
||||
let four = two + two;
|
||||
let eps_2 = F::sqrt(F::epsilon());
|
||||
|
||||
let a = two * p;
|
||||
let b = p * p - four * r;
|
||||
let c = -q * q;
|
||||
let f = solve_cubic(a, b, c);
|
||||
let a = f[0].sqrt() / two;
|
||||
let b = f[1].sqrt() / two;
|
||||
let c = f[2].sqrt() / two;
|
||||
|
||||
let mut res = (0..8)
|
||||
.map(|i| {
|
||||
let a = a * F::powi(-F::one(), i % 2);
|
||||
let b = b * F::powi(-F::one(), (i / 2) % 2);
|
||||
let c = c * F::powi(-F::one(), (i / 4) % 2);
|
||||
[-a - b - c, -a + b + c, a - b + c, a + b - c]
|
||||
})
|
||||
.map(|x| {
|
||||
let f = x
|
||||
.iter()
|
||||
.map(|t| (t * t * t * t + t * t * p + t * q + r).norm_sqr())
|
||||
.max_by(|x, y| x.partial_cmp(y).unwrap())
|
||||
.unwrap();
|
||||
(x, f)
|
||||
})
|
||||
.min_by(|x, y| x.1.partial_cmp(&y.1).unwrap())
|
||||
.unwrap()
|
||||
.0;
|
||||
// refinement by Newton method
|
||||
res.iter_mut().for_each(|x| {
|
||||
let mut f = *x * *x * *x * *x + *x * *x * p + *x * q + r;
|
||||
let mut f_prime = *x * *x * *x * four + *x * p * two + q;
|
||||
// SolveCraft: at most 100 steps. Near a double root (the eigenvalues of a
|
||||
// translation matrix) the iteration can cycle forever without meeting the tolerance.
|
||||
let mut steps = 0;
|
||||
while f.norm() > eps_2 * f_prime.norm() && steps < 100 {
|
||||
steps += 1;
|
||||
if f_prime.norm() < eps_2 {
|
||||
return;
|
||||
}
|
||||
*x -= f / f_prime;
|
||||
f_prime = *x * *x * *x * four + *x * p * two + q;
|
||||
f = *x * *x * *x * *x + *x * *x * p + *x * q + r;
|
||||
}
|
||||
});
|
||||
res
|
||||
}
|
||||
|
||||
#[cfg_attr(doc, katexit::katexit)]
|
||||
/// solve equation: $x^4 + ax^3 + bx^2 + cx + d = 0$.
|
||||
///
|
||||
/// # Examples
|
||||
///
|
||||
/// ```
|
||||
/// use matext4cgmath::solver;
|
||||
/// use num_complex::Complex;
|
||||
/// const EPS: f64 = 1.0e-10;
|
||||
///
|
||||
/// let mut res = solver::solve_quartic(1.0, -7.0, -1.0, 6.0);
|
||||
/// // Even in the case of real solutions, the order in the array is not guaranteed.
|
||||
/// res.sort_by(|x, y| x.re.partial_cmp(&y.re).unwrap());
|
||||
/// let ans = [Complex::from(-3.0), Complex::from(-1.0), Complex::from(1.0), Complex::from(2.0)];
|
||||
/// res.iter().zip(ans).for_each(|(x, y)| {
|
||||
/// assert!(Complex::norm(x * x * x * x + x * x * x - 7.0 * x * x - x + 6.0) < EPS);
|
||||
/// assert!(Complex::norm(x - y) < EPS, "{x} {y}");
|
||||
/// });
|
||||
/// ```
|
||||
pub fn solve_quartic<F: BaseFloat>(a: F, b: F, c: F, d: F) -> [Complex<F>; 4] {
|
||||
let one = F::one();
|
||||
let two = one + one;
|
||||
let three = one + two;
|
||||
let four = two + two;
|
||||
let six = two * three;
|
||||
let eight = four + four;
|
||||
|
||||
let a_4 = a / four;
|
||||
let p = b - six * a_4 * a_4;
|
||||
let q = c - two * b * a_4 + eight * a_4 * a_4 * a_4;
|
||||
let r = d - c * a_4 + b * a_4 * a_4 - three * a_4 * a_4 * a_4 * a_4;
|
||||
let mut res = pre_solve_quartic(p, q, r);
|
||||
res.iter_mut().for_each(|x| *x -= a_4);
|
||||
res
|
||||
}
|
||||
Reference in new issue
Block a user