CypherNaugh_0xandClaude Opus 5.5
5249ae84ef
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>
121 lines
3.6 KiB
Rust
121 lines
3.6 KiB
Rust
//! 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)>;
|
|
}
|