| | | 1 | | //======================================================================= |
| | | 2 | | // CylinderPairValuePolynomial.cs |
| | | 3 | | //======================================================================= |
| | | 4 | | // MIT License, Copyright (c) 2024-present David Oravsky (mrdav30) |
| | | 5 | | // See LICENSE file in the project root for full license information. |
| | | 6 | | //======================================================================= |
| | | 7 | | |
| | | 8 | | using System; |
| | | 9 | | |
| | | 10 | | namespace FixedMathSharp.Geometry; |
| | | 11 | | |
| | | 12 | | /// <summary> |
| | | 13 | | /// Constructs the exact, unnormalized squared-distance discriminant for two rims. |
| | | 14 | | /// </summary> |
| | | 15 | | internal static class CylinderPairValuePolynomial |
| | | 16 | | { |
| | | 17 | | internal const int InvariantCount = 11; |
| | | 18 | | internal const int CoefficientCount = 9; |
| | | 19 | | internal const int ScratchCoefficientCount = 73; |
| | | 20 | | |
| | | 21 | | /// <summary> |
| | | 22 | | /// Invariant slots are u², w², b², ra², rb², c², c.u, c.w, b.u, |
| | | 23 | | /// b.w, b.c. Every scalar has <paramref name="words"/> little-endian |
| | | 24 | | /// magnitude words and a separate canonical sign. Returns the degree, |
| | | 25 | | /// or -1 for zero, of three times the cubic-pencil discriminant. |
| | | 26 | | /// </summary> |
| | | 27 | | /// <remarks> |
| | | 28 | | /// The caller supplies disjoint input, output, and scratch storage: eleven, |
| | | 29 | | /// nine, and seventy-three coefficient slots, respectively. The basis u,w |
| | | 30 | | /// is orthogonal in the first circle's plane; independently reduced basis |
| | | 31 | | /// vectors are allowed. Geometry admission and physical S scaling belong |
| | | 32 | | /// to the caller. No coefficient content or repeated factor is removed. |
| | | 33 | | /// |
| | | 34 | | /// Every intermediate must fit words*64 bits. A conservative sufficient |
| | | 35 | | /// bound for arbitrary input magnitudes below 2^b, b at least one, is |
| | | 36 | | /// words*64 at least 32*b+128. A caller may instead use a tighter proven |
| | | 37 | | /// geometry-specific bound, including convolution sum carries. Magnitude |
| | | 38 | | /// operations do not detect truncation. There is no nested stack scratch; |
| | | 39 | | /// construction storage need not remain live during root isolation. |
| | | 40 | | /// An optional direction adds seventy-two scratch slots and includes its |
| | | 41 | | /// input magnitudes in b. The raw scale is W=3*discriminant for both value |
| | | 42 | | /// and derivative: allow up to two additional coefficient bits relative |
| | | 43 | | /// to a bound for the unscaled discriminant. Default direction storage |
| | | 44 | | /// performs no derivative convolutions. |
| | | 45 | | /// Scratch slots use min(words, ceil((32*b+128)/64)) words; input and |
| | | 46 | | /// retained output slots keep the caller's original words stride. |
| | | 47 | | /// </remarks> |
| | | 48 | | internal static int Build(ReadOnlySpan<ulong> invariants, |
| | | 49 | | ReadOnlySpan<sbyte> invariantSigns, int words, |
| | | 50 | | Span<ulong> coefficients, Span<sbyte> signs, |
| | | 51 | | Span<ulong> scratch, Span<sbyte> scratchSigns, |
| | | 52 | | CylinderPairValueDerivative derivative = default) |
| | | 53 | | { |
| | 6404 | 54 | | bool differentiate = !derivative.Invariants.IsEmpty; |
| | 6404 | 55 | | int inputBits = 1; |
| | 153696 | 56 | | for (int index = 0; index < InvariantCount; index++) |
| | 70444 | 57 | | inputBits = Math.Max(inputBits, |
| | 70444 | 58 | | WideArithmetic.GetMagnitudeBitLength(invariants.Slice(index * words, words))); |
| | 6404 | 59 | | if (differentiate) |
| | 71792 | 60 | | for (int index = 0; index < CylinderPairValueDerivative.InvariantCount; index++) |
| | 30768 | 61 | | inputBits = Math.Max(inputBits, |
| | 30768 | 62 | | WideArithmetic.GetMagnitudeBitLength(derivative.Invariants.Slice(index * words, words))); |
| | | 63 | | |
| | | 64 | | // Both the generic bound and the caller's geometry-specific capacity |
| | | 65 | | // bound every intermediate. Their minimum avoids scanning padded words |
| | | 66 | | // without changing retained coefficients or truncating dual products. |
| | 6404 | 67 | | int workWords = Math.Min(words, (32 * inputBits + 128 + 63) / 64); |
| | 6404 | 68 | | scratch = scratch.Slice(0, ScratchCoefficientCount * workWords); |
| | 6404 | 69 | | scratchSigns = scratchSigns.Slice(0, ScratchCoefficientCount); |
| | 6404 | 70 | | scratch.Clear(); |
| | 6404 | 71 | | scratchSigns.Clear(); |
| | 153696 | 72 | | for (int index = 0; index < InvariantCount; index++) |
| | 70444 | 73 | | invariants.Slice(index * words, workWords) |
| | 70444 | 74 | | .CopyTo(scratch.Slice(index * workWords, workWords)); |
| | 6404 | 75 | | invariantSigns.Slice(0, InvariantCount).CopyTo(scratchSigns); |
| | 6404 | 76 | | if (differentiate) |
| | | 77 | | { |
| | 5128 | 78 | | derivative.Scratch.Slice(0, CylinderPairValueDerivative.ScratchCoefficientCount * workWords).Clear(); |
| | 5128 | 79 | | derivative.ScratchSigns.Slice(0, CylinderPairValueDerivative.ScratchCoefficientCount).Clear(); |
| | 71792 | 80 | | for (int index = 0; index < CylinderPairValueDerivative.InvariantCount; index++) |
| | | 81 | | { |
| | 30768 | 82 | | int target = index < 5 ? index + 3 : 10; |
| | 30768 | 83 | | derivative.Invariants.Slice(index * words, workWords) |
| | 30768 | 84 | | .CopyTo(derivative.Scratch.Slice(target * workWords, workWords)); |
| | 30768 | 85 | | derivative.ScratchSigns[target] = derivative.InvariantSigns[index]; |
| | | 86 | | } |
| | | 87 | | } |
| | | 88 | | |
| | 6404 | 89 | | ConstructionScratch workspace = new(scratch, scratchSigns, workWords, derivative); |
| | 6404 | 90 | | Polynomial g = workspace.Slice(0, 1); |
| | 6404 | 91 | | Polynomial h = workspace.Slice(1, 1); |
| | 6404 | 92 | | Polynomial delta = workspace.Slice(2, 1); |
| | 6404 | 93 | | Polynomial rho = workspace.Slice(3, 1); |
| | 6404 | 94 | | Polynomial tau = workspace.Slice(4, 1); |
| | 6404 | 95 | | Polynomial cSquared = workspace.Slice(5, 1); |
| | 6404 | 96 | | Polynomial c0 = workspace.Slice(6, 1); |
| | 6404 | 97 | | Polynomial c1 = workspace.Slice(7, 1); |
| | 6404 | 98 | | Polynomial b0 = workspace.Slice(8, 1); |
| | 6404 | 99 | | Polynomial b1 = workspace.Slice(9, 1); |
| | 6404 | 100 | | Polynomial j = workspace.Slice(10, 1); |
| | | 101 | | |
| | 6404 | 102 | | Polynomial e = workspace.Slice(11, 1); |
| | 6404 | 103 | | Polynomial bigK = workspace.Slice(12, 1); |
| | 6404 | 104 | | Polynomial s0 = workspace.Slice(13, 1); |
| | 6404 | 105 | | Polynomial s1 = workspace.Slice(14, 1); |
| | 6404 | 106 | | Polynomial s2 = workspace.Slice(15, 1); |
| | 6404 | 107 | | Polynomial s3 = workspace.Slice(16, 1); |
| | 6404 | 108 | | Polynomial s4 = workspace.Slice(17, 1); |
| | 6404 | 109 | | Polynomial k = workspace.Slice(18, 2); |
| | 6404 | 110 | | Polynomial t0 = workspace.Slice(20, 2); |
| | 6404 | 111 | | Polynomial t1 = workspace.Slice(22, 2); |
| | 6404 | 112 | | Polynomial bigL = workspace.Slice(24, 3); |
| | 6404 | 113 | | Polynomial d = workspace.Slice(27, 3); |
| | 6404 | 114 | | Polynomial a = workspace.Slice(30, 1); |
| | 6404 | 115 | | Polynomial b = workspace.Slice(31, 3); |
| | 6404 | 116 | | Polynomial c = workspace.Slice(34, 3); |
| | 6404 | 117 | | Polynomial constant = workspace.Slice(37, 2); |
| | 6404 | 118 | | Polynomial p = workspace.Slice(39, 5); |
| | 6404 | 119 | | Polynomial r = workspace.Slice(44, 5); |
| | 6404 | 120 | | Polynomial n = workspace.Slice(49, 5); |
| | 6404 | 121 | | Polynomial temporary0 = workspace.Slice(54, 9); |
| | 6404 | 122 | | Polynomial temporary1 = workspace.Slice(63, 9); |
| | 6404 | 123 | | Span<ulong> product = scratch.Slice(72 * workWords, workWords); |
| | | 124 | | |
| | | 125 | | // e=(c0*b1-c1*b0)^2. |
| | 6404 | 126 | | Multiply(s0, c0, b1, product); |
| | 6404 | 127 | | Multiply(s1, c1, b0, product); |
| | 6404 | 128 | | Add(s0, s1, -1); |
| | 6404 | 129 | | Multiply(e, s0, s0, product); |
| | | 130 | | |
| | | 131 | | // K=delta*(g*c1^2+h*c0^2)+tau*(g*b1^2+h*b0^2). |
| | 6404 | 132 | | Multiply(s0, c1, c1, product); |
| | 6404 | 133 | | Multiply(s1, g, s0, product); |
| | 6404 | 134 | | Multiply(s0, c0, c0, product); |
| | 6404 | 135 | | Multiply(s2, h, s0, product); |
| | 6404 | 136 | | Add(s1, s2, 1); |
| | 6404 | 137 | | Multiply(s3, delta, s1, product); |
| | 6404 | 138 | | Multiply(s0, b1, b1, product); |
| | 6404 | 139 | | Multiply(s1, g, s0, product); |
| | 6404 | 140 | | Multiply(s0, b0, b0, product); |
| | 6404 | 141 | | Multiply(s2, h, s0, product); |
| | 6404 | 142 | | Add(s1, s2, 1); |
| | 6404 | 143 | | Multiply(s4, tau, s1, product); |
| | 6404 | 144 | | Add(bigK, s3, 1); |
| | 6404 | 145 | | Add(bigK, s4, 1); |
| | | 146 | | |
| | | 147 | | // k=c^2+rho-tau-S. |
| | 6404 | 148 | | Add(k, cSquared, 1); |
| | 6404 | 149 | | Add(k, rho, 1); |
| | 6404 | 150 | | Add(k, tau, -1); |
| | 6404 | 151 | | k.Magnitudes[workWords] = 1UL; |
| | 6404 | 152 | | k.Signs[1] = -1; |
| | | 153 | | |
| | | 154 | | // L=g*(2*j*c1-k*b1)^2+h*(2*j*c0-k*b0)^2. |
| | 6404 | 155 | | Multiply(s0, j, c1, product); |
| | 6404 | 156 | | Add(t0, s0, 2); |
| | 6404 | 157 | | Multiply(t1, k, b1, product); |
| | 6404 | 158 | | Add(t0, t1, -1); |
| | 6404 | 159 | | Multiply(temporary0, t0, t0, product); |
| | 6404 | 160 | | Multiply(bigL, g, temporary0, product); |
| | 6404 | 161 | | t0.Clear(); |
| | 6404 | 162 | | Multiply(s0, j, c0, product); |
| | 6404 | 163 | | Add(t0, s0, 2); |
| | 6404 | 164 | | Multiply(t1, k, b0, product); |
| | 6404 | 165 | | Add(t0, t1, -1); |
| | 6404 | 166 | | Multiply(temporary0, t0, t0, product); |
| | 6404 | 167 | | Multiply(temporary1, h, temporary0, product); |
| | 6404 | 168 | | Add(bigL, temporary1, 1); |
| | | 169 | | |
| | | 170 | | // Tzz=d=delta*(k^2-4*tau*S)+4*tau*j^2. |
| | 6404 | 171 | | Multiply(temporary0, k, k, product); |
| | 6404 | 172 | | Multiply(d, delta, temporary0, product); |
| | 6404 | 173 | | Multiply(s0, delta, tau, product); |
| | 6404 | 174 | | Add(d, s0, -4, 1); |
| | 6404 | 175 | | Multiply(s0, j, j, product); |
| | 6404 | 176 | | Multiply(s1, tau, s0, product); |
| | 6404 | 177 | | Add(d, s1, 4); |
| | | 178 | | |
| | | 179 | | // det(T-lambda*C)=A*lambda^3+B*lambda^2+C1*lambda+D. |
| | 6404 | 180 | | Multiply(s0, g, h, product); |
| | 6404 | 181 | | Multiply(a, s0, rho, product); |
| | 6404 | 182 | | Multiply(b, s0, d, product); |
| | 6404 | 183 | | Multiply(s1, rho, bigK, product); |
| | 6404 | 184 | | Add(b, s1, -4); |
| | | 185 | | |
| | 6404 | 186 | | temporary0.Clear(); |
| | 6404 | 187 | | Multiply(s0, rho, e, product); |
| | 6404 | 188 | | Add(temporary0, s0, 4); |
| | 6404 | 189 | | Add(temporary0, bigK, 4, 1); |
| | 6404 | 190 | | Add(temporary0, bigL, -1); |
| | 6404 | 191 | | Multiply(s0, delta, tau, product); |
| | 6404 | 192 | | Multiply(temporary1, s0, temporary0, product); |
| | 6404 | 193 | | Add(c, temporary1, 4); |
| | | 194 | | |
| | | 195 | | // The completed-square form makes D linear, not quadratic. |
| | 6404 | 196 | | Multiply(s0, delta, delta, product); |
| | 6404 | 197 | | Multiply(s1, tau, tau, product); |
| | 6404 | 198 | | Multiply(s2, s0, s1, product); |
| | 6404 | 199 | | Multiply(s3, s2, e, product); |
| | 6404 | 200 | | Add(constant, s3, -64, 1); |
| | | 201 | | |
| | | 202 | | // P=B^2-3*A*C1; R=C1^2-3*B*D; N=9*A*D-B*C1. |
| | 6404 | 203 | | Multiply(p, b, b, product); |
| | 6404 | 204 | | Multiply(temporary0, a, c, product); |
| | 6404 | 205 | | Add(p, temporary0, -3); |
| | 6404 | 206 | | Multiply(r, c, c, product); |
| | 6404 | 207 | | Multiply(temporary0, b, constant, product); |
| | 6404 | 208 | | Add(r, temporary0, -3); |
| | 6404 | 209 | | Multiply(temporary0, a, constant, product); |
| | 6404 | 210 | | Add(n, temporary0, 9); |
| | 6404 | 211 | | Multiply(temporary0, b, c, product); |
| | 6404 | 212 | | Add(n, temporary0, -1); |
| | | 213 | | |
| | | 214 | | // Three times the discriminant; preserve this exact envelope scale. |
| | 6404 | 215 | | Multiply(temporary0, p, r, product); |
| | 6404 | 216 | | Multiply(temporary1, n, n, product); |
| | 6404 | 217 | | Polynomial output = new(coefficients.Slice(0, CoefficientCount * words), |
| | 6404 | 218 | | signs.Slice(0, CoefficientCount), |
| | 6404 | 219 | | differentiate ? derivative.Coefficients.Slice(0, CoefficientCount * words) : default, |
| | 6404 | 220 | | differentiate ? derivative.Signs.Slice(0, CoefficientCount) : default); |
| | 6404 | 221 | | output.Clear(); |
| | 6404 | 222 | | Add(output, temporary0, 4); |
| | 6404 | 223 | | Add(output, temporary1, -1); |
| | 6404 | 224 | | int degree = CoefficientCount - 1; |
| | 6424 | 225 | | while (degree >= 0 && output.Signs[degree] == 0) |
| | 20 | 226 | | degree--; |
| | 6404 | 227 | | return degree; |
| | | 228 | | } |
| | | 229 | | |
| | | 230 | | // These are construction views, not retained roots or a second arithmetic |
| | | 231 | | // implementation. Scalar products and signed sums use WideArithmetic; |
| | | 232 | | // every destination in the schedule is disjoint from its sources. |
| | | 233 | | private readonly ref struct Polynomial |
| | | 234 | | { |
| | | 235 | | internal readonly Span<ulong> Magnitudes; |
| | | 236 | | internal readonly Span<sbyte> Signs; |
| | | 237 | | internal readonly Span<ulong> DerivativeMagnitudes; |
| | | 238 | | internal readonly Span<sbyte> DerivativeSigns; |
| | | 239 | | |
| | | 240 | | internal Polynomial(Span<ulong> magnitudes, Span<sbyte> signs, |
| | | 241 | | Span<ulong> derivativeMagnitudes = default, Span<sbyte> derivativeSigns = default) |
| | | 242 | | { |
| | 211332 | 243 | | Magnitudes = magnitudes; |
| | 211332 | 244 | | Signs = signs; |
| | 211332 | 245 | | DerivativeMagnitudes = derivativeMagnitudes; |
| | 211332 | 246 | | DerivativeSigns = derivativeSigns; |
| | 211332 | 247 | | } |
| | | 248 | | |
| | | 249 | | internal void Clear() |
| | | 250 | | { |
| | 19212 | 251 | | Magnitudes.Clear(); |
| | 19212 | 252 | | Signs.Clear(); |
| | 19212 | 253 | | DerivativeMagnitudes.Clear(); |
| | 19212 | 254 | | DerivativeSigns.Clear(); |
| | 19212 | 255 | | } |
| | | 256 | | } |
| | | 257 | | |
| | | 258 | | private readonly ref struct ConstructionScratch |
| | | 259 | | { |
| | | 260 | | private readonly Span<ulong> _magnitudes; |
| | | 261 | | private readonly Span<sbyte> _signs; |
| | | 262 | | private readonly int _words; |
| | | 263 | | private readonly CylinderPairValueDerivative _derivative; |
| | | 264 | | |
| | | 265 | | internal ConstructionScratch(Span<ulong> magnitudes, Span<sbyte> signs, |
| | | 266 | | int words, CylinderPairValueDerivative derivative) |
| | | 267 | | { |
| | 6404 | 268 | | _magnitudes = magnitudes; |
| | 6404 | 269 | | _signs = signs; |
| | 6404 | 270 | | _words = words; |
| | 6404 | 271 | | _derivative = derivative; |
| | 6404 | 272 | | } |
| | | 273 | | |
| | | 274 | | internal Polynomial Slice(int first, int count) => |
| | 204928 | 275 | | new(_magnitudes.Slice(first * _words, count * _words), |
| | 204928 | 276 | | _signs.Slice(first, count), |
| | 204928 | 277 | | _derivative.Invariants.IsEmpty ? default : |
| | 204928 | 278 | | _derivative.Scratch.Slice(first * _words, count * _words), |
| | 204928 | 279 | | _derivative.Invariants.IsEmpty ? default : |
| | 204928 | 280 | | _derivative.ScratchSigns.Slice(first, count)); |
| | | 281 | | } |
| | | 282 | | |
| | | 283 | | private static void Multiply(Polynomial destination, Polynomial left, |
| | | 284 | | Polynomial right, Span<ulong> product) |
| | | 285 | | { |
| | 288180 | 286 | | WideFiniteAxisIntersection.MultiplyFiniteAxisPolynomials( |
| | 288180 | 287 | | left.Magnitudes, left.Signs, right.Magnitudes, right.Signs, |
| | 288180 | 288 | | destination.Magnitudes, destination.Signs, product); |
| | 288180 | 289 | | if (!destination.DerivativeSigns.IsEmpty) |
| | | 290 | | { |
| | 230760 | 291 | | WideFiniteAxisIntersection.MultiplyFiniteAxisPolynomials( |
| | 230760 | 292 | | left.DerivativeMagnitudes, left.DerivativeSigns, right.Magnitudes, right.Signs, |
| | 230760 | 293 | | destination.DerivativeMagnitudes, destination.DerivativeSigns, product); |
| | 230760 | 294 | | WideFiniteAxisIntersection.MultiplyFiniteAxisPolynomials( |
| | 230760 | 295 | | left.Magnitudes, left.Signs, right.DerivativeMagnitudes, right.DerivativeSigns, |
| | 230760 | 296 | | destination.DerivativeMagnitudes, destination.DerivativeSigns, product, accumulate: true); |
| | | 297 | | } |
| | 288180 | 298 | | } |
| | | 299 | | |
| | | 300 | | private static void Add(Polynomial destination, Polynomial source, |
| | | 301 | | int multiplier, int degreeOffset = 0) |
| | | 302 | | { |
| | 172908 | 303 | | WideFiniteAxisIntersection.AddFiniteAxisPolynomial( |
| | 172908 | 304 | | source.Magnitudes, source.Signs, destination.Magnitudes, destination.Signs, |
| | 172908 | 305 | | multiplier, degreeOffset); |
| | 172908 | 306 | | if (!destination.DerivativeSigns.IsEmpty) |
| | 138456 | 307 | | WideFiniteAxisIntersection.AddFiniteAxisPolynomial( |
| | 138456 | 308 | | source.DerivativeMagnitudes, source.DerivativeSigns, |
| | 138456 | 309 | | destination.DerivativeMagnitudes, destination.DerivativeSigns, |
| | 138456 | 310 | | multiplier, degreeOffset); |
| | 172908 | 311 | | } |
| | | 312 | | } |