From 5249ae84ef8abb6b5c9744549cbe6953d365aa3c Mon Sep 17 00:00:00 2001 From: CypherNaugh_0x <9931495+CypherNaught-0x@users.noreply.github.com> Date: Sat, 10 Oct 2026 23:31:46 +0200 Subject: [PATCH] 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 --- ATTRIBUTION.md | 1 + Cargo.lock | 2 - Cargo.toml | 1 + crates/engine/src/boolean_regression_tests.rs | 5 +- crates/kernel/src/tests.rs | 16 ++ .../oblique_cut_of_cylinder_hang.json | 15 ++ vendor/README.md | 12 + vendor/matext4cgmath/Cargo.toml | 40 +++ vendor/matext4cgmath/LICENSE | 201 ++++++++++++++ vendor/matext4cgmath/README.md | 5 + vendor/matext4cgmath/src/eigens.rs | 153 +++++++++++ vendor/matext4cgmath/src/exp_decomp.rs | 106 ++++++++ vendor/matext4cgmath/src/lib.rs | 120 +++++++++ vendor/matext4cgmath/src/solver.rs | 246 ++++++++++++++++++ 14 files changed, 920 insertions(+), 3 deletions(-) create mode 100644 examples/regressions/oblique_cut_of_cylinder_hang.json create mode 100644 vendor/matext4cgmath/Cargo.toml create mode 100644 vendor/matext4cgmath/LICENSE create mode 100644 vendor/matext4cgmath/README.md create mode 100644 vendor/matext4cgmath/src/eigens.rs create mode 100644 vendor/matext4cgmath/src/exp_decomp.rs create mode 100644 vendor/matext4cgmath/src/lib.rs create mode 100644 vendor/matext4cgmath/src/solver.rs diff --git a/ATTRIBUTION.md b/ATTRIBUTION.md index 4880cd91..92321a60 100644 --- a/ATTRIBUTION.md +++ b/ATTRIBUTION.md @@ -49,3 +49,4 @@ Generated-in-code assets are original and have no file to list: | vendor/truck-meshalgo | truck-meshalgo 0.4.0, https://github.com/ricosjp/truck (RICOS Co. Ltd.) | Apache-2.0 | meshing: parameter searches kept inside the surface domain (vendor/README.md) | | vendor/truck-stepio | truck-stepio 0.3.0, https://github.com/ricosjp/truck (RICOS Co. Ltd.) | Apache-2.0 | STEP writer: intersection-curve surfaces written at unique ids (vendor/README.md) | | vendor/truck-shapeops | truck-shapeops 0.4.0, https://github.com/ricosjp/truck (RICOS Co. Ltd.) | Apache-2.0 | robust meshing and nearest-parameter fallback in booleans; trace macro (vendor/README.md) | +| vendor/matext4cgmath | matext4cgmath 0.1.0, https://github.com/ricosjp/matext4cgmath (RICOS Co. Ltd.) | Apache-2.0 | root refinement stops after 100 Newton steps (vendor/README.md) | diff --git a/Cargo.lock b/Cargo.lock index 699ac874..3289c02f 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -1753,8 +1753,6 @@ checksum = "f9f8bd3e56ce4dfc153cf470fffbfa98c7620958b312ca5c3a4b8d5181fd13c6" [[package]] name = "matext4cgmath" version = "0.1.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "e6541e181de37f70f0aceb25441823a3d0efa9cc1d23f475a0d3926678949178" dependencies = [ "cgmath", "katexit", diff --git a/Cargo.toml b/Cargo.toml index 693cb1ce..f6190e3c 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -98,3 +98,4 @@ debug = 0 truck-shapeops = { path = "vendor/truck-shapeops" } truck-stepio = { path = "vendor/truck-stepio" } truck-meshalgo = { path = "vendor/truck-meshalgo" } +matext4cgmath = { path = "vendor/matext4cgmath" } diff --git a/crates/engine/src/boolean_regression_tests.rs b/crates/engine/src/boolean_regression_tests.rs index 54f0ee80..8d8492ca 100644 --- a/crates/engine/src/boolean_regression_tests.rs +++ b/crates/engine/src/boolean_regression_tests.rs @@ -35,7 +35,10 @@ fn shallow_loft_cut_from_the_top_face() { #[test] fn oblique_cut_through_a_small_arc_matches_a_split() { - let scripts = [include_str!("../../../examples/regressions/oblique_cut_of_cylinder.json")]; + let scripts = [ + include_str!("../../../examples/regressions/oblique_cut_of_cylinder.json"), + include_str!("../../../examples/regressions/oblique_cut_of_cylinder_hang.json"), + ]; // The same bar split by the prism's tilted top: the piece above it is the cut's result. let mut split = Session::default(); let script: Value = serde_json::from_str(scripts[0]).unwrap(); diff --git a/crates/kernel/src/tests.rs b/crates/kernel/src/tests.rs index 869a2639..2dbd1c58 100644 --- a/crates/kernel/src/tests.rs +++ b/crates/kernel/src/tests.rs @@ -258,6 +258,22 @@ fn loft_and_sweep() { assert!(rel(measure(&t).unwrap().volume, expect) < 1e-3, "{} vs {expect}", measure(&t).unwrap().volume); } +#[test] +fn moving_a_body_with_an_intersection_curve_edge_ends() { + // A cylinder cut by a tilted box keeps an elliptical intersection curve as an edge. Moving + // such a body scales the curve's tolerance by the matrix norm, whose root refinement cycled + // forever for about one translation in eight (this one included): the app hung. + let cyl = cylinder(Vec3::ZERO, Vec3::Z, 5.0, 40.0).unwrap(); + let block = box_solid(Vec3::new(-20.0, -20.0, 20.0), Vec3::new(20.0, 20.0, 60.0)).unwrap(); + let block = crate::ops::transform(&block, Vec3::ZERO, Vec3::new(0.0, 0.0, 20.0), Vec3::X, 0.4).unwrap(); + let cut = boolean(&cyl, &block, BoolOp::Cut).unwrap().unwrap(); + assert!(cut.solid.edge_iter().any(|e| matches!(e.curve(), truck_modeling::Curve::IntersectionCurve(_)))); + let t = Vec3::new(0.6911, -0.3337, -0.6127) * 2.334; + let moved = crate::ops::transform(&cut, t, Vec3::ZERO, Vec3::Z, 0.0).unwrap(); + let (v0, v1) = (measure(&cut).unwrap().volume, measure(&moved).unwrap().volume); + assert!(rel(v0, v1) < 1e-6, "{v0} vs {v1}"); +} + #[test] fn planar_booleans_with_coincident_faces() { // L-bracket: a wall flush with three sides of its base. diff --git a/examples/regressions/oblique_cut_of_cylinder_hang.json b/examples/regressions/oblique_cut_of_cylinder_hang.json new file mode 100644 index 00000000..76aea7db --- /dev/null +++ b/examples/regressions/oblique_cut_of_cylinder_hang.json @@ -0,0 +1,15 @@ +{ + "description": "The oblique cut of oblique_cut_of_cylinder.json with the prism running past the bar's end. Moving a boolean result back after a shifted retry computed a matrix norm whose root refinement never ended: the app hung at full CPU. Result: the same trimmed bar.", + "commands": [ + {"command": "sketch.create", "params": {"plane": "XZ", "offset": -462, "name": "Profile"}}, + {"command": "sketch.line", "params": {"points": [[147.4407, -849.95], [147.4407, -851.0], [151.3952, -851.0], [151.3952, -847.4709], [150.3952, -847.4709]]}}, + {"command": "sketch.arc.center_point", "params": {"center": [147.4407, -846.95], "start": "l1.start", "end": "l4.end"}}, + {"command": "sketch.finish", "params": {}}, + {"command": "solid.extrude", "params": {"sketch": "Profile", "profiles": "all", "direction": "positive", "start_offset": -1, "distance": 437.52, "operation": "new", "body_name": "Tool"}}, + {"command": "sketch.create", "params": {"plane": "YZ", "name": "Clip"}}, + {"command": "sketch.line", "params": {"points": [[-480, 856.9979999999999], [-432, 838.47], [-432, 820], [-480, 820]], "closed": true}}, + {"command": "sketch.finish", "params": {}}, + {"command": "solid.extrude", "params": {"sketch": "Clip", "profiles": "all", "start_offset": -10, "distance": 920, "operation": "new", "body_name": "ClipBody"}}, + {"command": "solid.combine", "params": {"target": "Tool", "tools": ["ClipBody"], "operation": "cut", "keep_tools": true}} + ] +} diff --git a/vendor/README.md b/vendor/README.md index 229a7e85..aa10eece 100644 --- a/vendor/README.md +++ b/vendor/README.md @@ -49,3 +49,15 @@ Changes: direction (long triangles from a trim loop to the grid's far ends cut through cylinders). - `lib.rs`: compiler warnings no longer fail the build. Example, test and bench targets are dropped (their sources are not vendored). + +## matext4cgmath 0.1.0 (Apache-2.0, © RICOS Co. Ltd., https://github.com/ricosjp/matext4cgmath) + +truck-geometry uses it for a matrix's operator norm (an intersection curve scales its tolerance +by it whenever it is moved). + +Changes: +- `solver.rs`: the Newton refinement of cubic and quartic roots stops after 100 steps. Near a + double root it could cycle forever: moving a body whose edges a boolean made (the retries + move both bodies and the result back) hung the app. +- `lib.rs`: compiler warnings no longer fail the build. Test targets are dropped (their sources + are not vendored). diff --git a/vendor/matext4cgmath/Cargo.toml b/vendor/matext4cgmath/Cargo.toml new file mode 100644 index 00000000..ff2370c9 --- /dev/null +++ b/vendor/matext4cgmath/Cargo.toml @@ -0,0 +1,40 @@ +# THIS FILE IS AUTOMATICALLY GENERATED BY CARGO +# +# When uploading crates to the registry Cargo will automatically +# "normalize" Cargo.toml files for maximal compatibility +# with all versions of Cargo and also rewrite `path` dependencies +# to registry (e.g., crates.io) dependencies. +# +# If you are reading this file be aware that the original Cargo.toml +# will likely look very different (and much more reasonable). +# See Cargo.toml.orig for the original contents. + +[package] +edition = "2021" +name = "matext4cgmath" +version = "0.1.0" +authors = ["Yoshinori Tanimura "] +description = "Unofficial third-party cgmath extensions for calculate eigenvalues, operator norms and Iwasawa decomposition." +homepage = "https://github.com/ricosjp/matext4cgmath" +keywords = [ + "cgmath", + "graphics", + "linear-algebra", +] +categories = ["graphics"] +license = "Apache-2.0" +repository = "https://github.com/ricosjp/matext4cgmath" +resolver = "2" +autotests = false +autoexamples = false +autobenches = false + +[dependencies.cgmath] +version = "0.18.0" + +[dependencies.katexit] +version = "0.1.2" + +[dependencies.num-complex] +version = "0.4.0" + diff --git a/vendor/matext4cgmath/LICENSE b/vendor/matext4cgmath/LICENSE new file mode 100644 index 00000000..1ad4a8f7 --- /dev/null +++ b/vendor/matext4cgmath/LICENSE @@ -0,0 +1,201 @@ + Apache License + Version 2.0, January 2004 + http://www.apache.org/licenses/ + + TERMS AND CONDITIONS FOR USE, REPRODUCTION, AND DISTRIBUTION + + 1. Definitions. + + "License" shall mean the terms and conditions for use, reproduction, + and distribution as defined by Sections 1 through 9 of this document. + + "Licensor" shall mean the copyright owner or entity authorized by + the copyright owner that is granting the License. + + "Legal Entity" shall mean the union of the acting entity and all + other entities that control, are controlled by, or are under common + control with that entity. For the purposes of this definition, + "control" means (i) the power, direct or indirect, to cause the + direction or management of such entity, whether by contract or + otherwise, or (ii) ownership of fifty percent (50%) or more of the + outstanding shares, or (iii) beneficial ownership of such entity. + + "You" (or "Your") shall mean an individual or Legal Entity + exercising permissions granted by this License. + + "Source" form shall mean the preferred form for making modifications, + including but not limited to software source code, documentation + source, and configuration files. + + "Object" form shall mean any form resulting from mechanical + transformation or translation of a Source form, including but + not limited to compiled object code, generated documentation, + and conversions to other media types. + + "Work" shall mean the work of authorship, whether in Source or + Object form, made available under the License, as indicated by a + copyright notice that is included in or attached to the work + (an example is provided in the Appendix below). + + "Derivative Works" shall mean any work, whether in Source or Object + form, that is based on (or derived from) the Work and for which the + editorial revisions, annotations, elaborations, or other modifications + represent, as a whole, an original work of authorship. For the purposes + of this License, Derivative Works shall not include works that remain + separable from, or merely link (or bind by name) to the interfaces of, + the Work and Derivative Works thereof. + + "Contribution" shall mean any work of authorship, including + the original version of the Work and any modifications or additions + to that Work or Derivative Works thereof, that is intentionally + submitted to Licensor for inclusion in the Work by the copyright owner + or by an individual or Legal Entity authorized to submit on behalf of + the copyright owner. For the purposes of this definition, "submitted" + means any form of electronic, verbal, or written communication sent + to the Licensor or its representatives, including but not limited to + communication on electronic mailing lists, source code control systems, + and issue tracking systems that are managed by, or on behalf of, the + Licensor for the purpose of discussing and improving the Work, but + excluding communication that is conspicuously marked or otherwise + designated in writing by the copyright owner as "Not a Contribution." + + "Contributor" shall mean Licensor and any individual or Legal Entity + on behalf of whom a Contribution has been received by Licensor and + subsequently incorporated within the Work. + + 2. Grant of Copyright License. Subject to the terms and conditions of + this License, each Contributor hereby grants to You a perpetual, + worldwide, non-exclusive, no-charge, royalty-free, irrevocable + copyright license to reproduce, prepare Derivative Works of, + publicly display, publicly perform, sublicense, and distribute the + Work and such Derivative Works in Source or Object form. + + 3. Grant of Patent License. Subject to the terms and conditions of + this License, each Contributor hereby grants to You a perpetual, + worldwide, non-exclusive, no-charge, royalty-free, irrevocable + (except as stated in this section) patent license to make, have made, + use, offer to sell, sell, import, and otherwise transfer the Work, + where such license applies only to those patent claims licensable + by such Contributor that are necessarily infringed by their + Contribution(s) alone or by combination of their Contribution(s) + with the Work to which such Contribution(s) was submitted. If You + institute patent litigation against any entity (including a + cross-claim or counterclaim in a lawsuit) alleging that the Work + or a Contribution incorporated within the Work constitutes direct + or contributory patent infringement, then any patent licenses + granted to You under this License for that Work shall terminate + as of the date such litigation is filed. + + 4. Redistribution. You may reproduce and distribute copies of the + Work or Derivative Works thereof in any medium, with or without + modifications, and in Source or Object form, provided that You + meet the following conditions: + + (a) You must give any other recipients of the Work or + Derivative Works a copy of this License; and + + (b) You must cause any modified files to carry prominent notices + stating that You changed the files; and + + (c) You must retain, in the Source form of any Derivative Works + that You distribute, all copyright, patent, trademark, and + attribution notices from the Source form of the Work, + excluding those notices that do not pertain to any part of + the Derivative Works; and + + (d) If the Work includes a "NOTICE" text file as part of its + distribution, then any Derivative Works that You distribute must + include a readable copy of the attribution notices contained + within such NOTICE file, excluding those notices that do not + pertain to any part of the Derivative Works, in at least one + of the following places: within a NOTICE text file distributed + as part of the Derivative Works; within the Source form or + documentation, if provided along with the Derivative Works; or, + within a display generated by the Derivative Works, if and + wherever such third-party notices normally appear. The contents + of the NOTICE file are for informational purposes only and + do not modify the License. You may add Your own attribution + notices within Derivative Works that You distribute, alongside + or as an addendum to the NOTICE text from the Work, provided + that such additional attribution notices cannot be construed + as modifying the License. + + You may add Your own copyright statement to Your modifications and + may provide additional or different license terms and conditions + for use, reproduction, or distribution of Your modifications, or + for any such Derivative Works as a whole, provided Your use, + reproduction, and distribution of the Work otherwise complies with + the conditions stated in this License. + + 5. Submission of Contributions. Unless You explicitly state otherwise, + any Contribution intentionally submitted for inclusion in the Work + by You to the Licensor shall be under the terms and conditions of + this License, without any additional terms or conditions. + Notwithstanding the above, nothing herein shall supersede or modify + the terms of any separate license agreement you may have executed + with Licensor regarding such Contributions. + + 6. Trademarks. This License does not grant permission to use the trade + names, trademarks, service marks, or product names of the Licensor, + except as required for reasonable and customary use in describing the + origin of the Work and reproducing the content of the NOTICE file. + + 7. Disclaimer of Warranty. Unless required by applicable law or + agreed to in writing, Licensor provides the Work (and each + Contributor provides its Contributions) on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or + implied, including, without limitation, any warranties or conditions + of TITLE, NON-INFRINGEMENT, MERCHANTABILITY, or FITNESS FOR A + PARTICULAR PURPOSE. You are solely responsible for determining the + appropriateness of using or redistributing the Work and assume any + risks associated with Your exercise of permissions under this License. + + 8. Limitation of Liability. In no event and under no legal theory, + whether in tort (including negligence), contract, or otherwise, + unless required by applicable law (such as deliberate and grossly + negligent acts) or agreed to in writing, shall any Contributor be + liable to You for damages, including any direct, indirect, special, + incidental, or consequential damages of any character arising as a + result of this License or out of the use or inability to use the + Work (including but not limited to damages for loss of goodwill, + work stoppage, computer failure or malfunction, or any and all + other commercial damages or losses), even if such Contributor + has been advised of the possibility of such damages. + + 9. Accepting Warranty or Additional Liability. While redistributing + the Work or Derivative Works thereof, You may choose to offer, + and charge a fee for, acceptance of support, warranty, indemnity, + or other liability obligations and/or rights consistent with this + License. However, in accepting such obligations, You may act only + on Your own behalf and on Your sole responsibility, not on behalf + of any other Contributor, and only if You agree to indemnify, + defend, and hold each Contributor harmless for any liability + incurred by, or claims asserted against, such Contributor by reason + of your accepting any such warranty or additional liability. + + END OF TERMS AND CONDITIONS + + APPENDIX: How to apply the Apache License to your work. + + To apply the Apache License to your work, attach the following + boilerplate notice, with the fields enclosed by brackets "[]" + replaced with your own identifying information. (Don't include + the brackets!) The text should be enclosed in the appropriate + comment syntax for the file format. We also recommend that a + file or class name and description of purpose be included on the + same "printed page" as the copyright notice for easier + identification within third-party archives. + + Copyright 2022 Research Institute for Computational Science Co. Ltd + + Licensed under the Apache License, Version 2.0 (the "License"); + you may not use this file except in compliance with the License. + You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + + Unless required by applicable law or agreed to in writing, software + distributed under the License is distributed on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + See the License for the specific language governing permissions and + limitations under the License. diff --git a/vendor/matext4cgmath/README.md b/vendor/matext4cgmath/README.md new file mode 100644 index 00000000..c2ce89d7 --- /dev/null +++ b/vendor/matext4cgmath/README.md @@ -0,0 +1,5 @@ +# matext4cgmath + +[![Crates.io](https://img.shields.io/crates/v/matext4cgmath.svg)](https://crates.io/crates/matext4cgmath) [![Docs.rs](https://docs.rs/matext4cgmath/badge.svg)](https://docs.rs/matext4cgmath) + +Unofficial third-party cgmath extensions for calculate eigenvalues, operator norms and Iwasawa decomposition. diff --git a/vendor/matext4cgmath/src/eigens.rs b/vendor/matext4cgmath/src/eigens.rs new file mode 100644 index 00000000..5f7db7e9 --- /dev/null +++ b/vendor/matext4cgmath/src/eigens.rs @@ -0,0 +1,153 @@ +use crate::*; + +impl EigenValues for Matrix2 { + type EigenValues = [Complex; 2]; + fn eigenvalues(self) -> [Complex; 2] { + solver::solve_quadratic(-self.trace(), self.determinant()) + } +} + +impl EigenValues for Matrix3 { + type EigenValues = [Complex; 3]; + fn eigenvalues(self) -> [Complex; 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 EigenValues for Matrix4 { + type EigenValues = [Complex; 4]; + fn eigenvalues(self) -> [Complex; 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 OperatorNorm for Matrix2 { + #[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 OperatorNorm for Matrix3 { + #[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 OperatorNorm for Matrix4 { + #[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, + )) + } +} diff --git a/vendor/matext4cgmath/src/exp_decomp.rs b/vendor/matext4cgmath/src/exp_decomp.rs new file mode 100644 index 00000000..20a90f21 --- /dev/null +++ b/vendor/matext4cgmath/src/exp_decomp.rs @@ -0,0 +1,106 @@ +use crate::*; + +impl Exponential for Matrix2 {} +impl Exponential for Matrix3 {} +impl Exponential for Matrix4 {} + +impl Decomposition for Matrix2 { + 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 Decomposition for Matrix3 { + 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 Decomposition for Matrix4 { + 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 + } +} diff --git a/vendor/matext4cgmath/src/lib.rs b/vendor/matext4cgmath/src/lib.rs new file mode 100644 index 00000000..08c8d2a2 --- /dev/null +++ b/vendor/matext4cgmath/src/lib.rs @@ -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 + 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 = ::epsilon(); + let mut a = self; + let mut res = Self::one(); + for i in 2_u16..=64 { + res += a; + a = a * self / ::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)>; +} diff --git a/vendor/matext4cgmath/src/solver.rs b/vendor/matext4cgmath/src/solver.rs new file mode 100644 index 00000000..8713ce3d --- /dev/null +++ b/vendor/matext4cgmath/src/solver.rs @@ -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(a: F, b: F) -> [Complex; 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(p: F, q: F) -> [Complex; 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(a: F, b: F, c: F) -> [Complex; 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(p: F, q: F, r: F) -> [Complex; 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(a: F, b: F, c: F, d: F) -> [Complex; 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 +}