| | | 1 | | //======================================================================= |
| | | 2 | | // WideFiniteAxisIntersection.RoundedCylinderRootSolver.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 | | /// <content> |
| | | 13 | | /// Implements a wide-precision root solver for rounded-cylinder (capsule-like) torus |
| | | 14 | | /// intersection tests, using a Sturm sequence built over arbitrary-precision integer |
| | | 15 | | /// coefficients to isolate and count roots within a bounded segment length. |
| | | 16 | | /// </content> |
| | | 17 | | internal static partial class WideFiniteAxisIntersection |
| | | 18 | | { |
| | | 19 | | private const int RoundedCylinderPolynomialCount = 5; |
| | | 20 | | private const int RoundedCylinderCoefficientCount = 5; |
| | | 21 | | // The largest factorized quartic subresultant transient is below 2,772 bits. |
| | | 22 | | private const int RoundedCylinderWideLimbCount = 44; |
| | | 23 | | |
| | | 24 | | /// <summary> |
| | | 25 | | /// Counts real roots for quartics whose coefficient growth is bounded by |
| | | 26 | | /// the finite-axis and planar-projection derivations in this assembly. |
| | | 27 | | /// </summary> |
| | | 28 | | internal static int CountProjectionDomainQuarticRoots( |
| | | 29 | | Signed832 constant, |
| | | 30 | | Signed832 linear, |
| | | 31 | | Signed832 quadratic, |
| | | 32 | | Signed832 cubic, |
| | | 33 | | Signed832 quartic) |
| | | 34 | | { |
| | 114 | 35 | | int leadingSign = quartic.Sign; |
| | 114 | 36 | | if (leadingSign == 0) |
| | 1 | 37 | | throw new ArgumentOutOfRangeException(nameof(quartic)); |
| | | 38 | | |
| | 113 | 39 | | if (leadingSign < 0) |
| | | 40 | | { |
| | 4 | 41 | | constant = WideArithmetic.SubtractSigned832(default, constant); |
| | 4 | 42 | | linear = WideArithmetic.SubtractSigned832(default, linear); |
| | 4 | 43 | | quadratic = WideArithmetic.SubtractSigned832(default, quadratic); |
| | 4 | 44 | | cubic = WideArithmetic.SubtractSigned832(default, cubic); |
| | 4 | 45 | | quartic = WideArithmetic.SubtractSigned832(default, quartic); |
| | | 46 | | } |
| | | 47 | | |
| | 113 | 48 | | Span<ulong> coefficients = stackalloc ulong[ |
| | 113 | 49 | | RoundedCylinderPolynomialCount |
| | 113 | 50 | | * RoundedCylinderCoefficientCount |
| | 113 | 51 | | * RoundedCylinderWideLimbCount]; |
| | 113 | 52 | | Span<sbyte> signs = stackalloc sbyte[ |
| | 113 | 53 | | RoundedCylinderPolynomialCount * RoundedCylinderCoefficientCount]; |
| | 113 | 54 | | Span<int> degrees = stackalloc int[RoundedCylinderPolynomialCount]; |
| | 113 | 55 | | coefficients.Clear(); |
| | 113 | 56 | | signs.Clear(); |
| | 113 | 57 | | degrees.Clear(); |
| | 113 | 58 | | ImportRoundedCylinderCoefficient(constant, coefficients, signs, 0, 0); |
| | 113 | 59 | | ImportRoundedCylinderCoefficient(linear, coefficients, signs, 0, 1); |
| | 113 | 60 | | ImportRoundedCylinderCoefficient(quadratic, coefficients, signs, 0, 2); |
| | 113 | 61 | | ImportRoundedCylinderCoefficient(cubic, coefficients, signs, 0, 3); |
| | 113 | 62 | | ImportRoundedCylinderCoefficient(quartic, coefficients, signs, 0, 4); |
| | 113 | 63 | | degrees[0] = 4; |
| | 1130 | 64 | | for (int coefficientIndex = 1; coefficientIndex <= 4; coefficientIndex++) |
| | | 65 | | { |
| | 452 | 66 | | int destinationIndex = coefficientIndex - 1; |
| | 452 | 67 | | MultiplyRoundedCylinderWideByWord( |
| | 452 | 68 | | GetRoundedCylinderCoefficient( |
| | 452 | 69 | | coefficients, |
| | 452 | 70 | | 0, |
| | 452 | 71 | | coefficientIndex), |
| | 452 | 72 | | (ulong)coefficientIndex, |
| | 452 | 73 | | GetRoundedCylinderCoefficient( |
| | 452 | 74 | | coefficients, |
| | 452 | 75 | | 1, |
| | 452 | 76 | | destinationIndex)); |
| | 452 | 77 | | signs[GetRoundedCylinderCoefficientIndex(1, destinationIndex)] = |
| | 452 | 78 | | signs[GetRoundedCylinderCoefficientIndex( |
| | 452 | 79 | | 0, |
| | 452 | 80 | | coefficientIndex)]; |
| | | 81 | | } |
| | 113 | 82 | | degrees[1] = 3; |
| | | 83 | | |
| | 113 | 84 | | int sequenceCount = BuildRoundedCylinderQuarticSturmSequence( |
| | 113 | 85 | | coefficients, |
| | 113 | 86 | | signs, |
| | 113 | 87 | | degrees); |
| | 113 | 88 | | return GetRoundedCylinderInfinityVariations( |
| | 113 | 89 | | signs, |
| | 113 | 90 | | degrees, |
| | 113 | 91 | | sequenceCount, |
| | 113 | 92 | | negative: true) |
| | 113 | 93 | | - GetRoundedCylinderInfinityVariations( |
| | 113 | 94 | | signs, |
| | 113 | 95 | | degrees, |
| | 113 | 96 | | sequenceCount, |
| | 113 | 97 | | negative: false); |
| | | 98 | | } |
| | | 99 | | |
| | | 100 | | private static int GetRoundedCylinderInfinityVariations( |
| | | 101 | | ReadOnlySpan<sbyte> signs, |
| | | 102 | | ReadOnlySpan<int> degrees, |
| | | 103 | | int sequenceCount, |
| | | 104 | | bool negative) |
| | | 105 | | { |
| | 226 | 106 | | int variations = 0; |
| | 226 | 107 | | int previousSign = 0; |
| | 226 | 108 | | for (int polynomialIndex = 0; |
| | 1320 | 109 | | polynomialIndex < sequenceCount; |
| | 1094 | 110 | | polynomialIndex++) |
| | | 111 | | { |
| | 1094 | 112 | | int degree = degrees[polynomialIndex]; |
| | 1094 | 113 | | int sign = signs[ |
| | 1094 | 114 | | GetRoundedCylinderCoefficientIndex( |
| | 1094 | 115 | | polynomialIndex, |
| | 1094 | 116 | | degree)]; |
| | 1094 | 117 | | if (negative && (degree & 1) != 0) |
| | 221 | 118 | | sign = -sign; |
| | 1094 | 119 | | if (previousSign != 0 && sign != previousSign) |
| | 434 | 120 | | variations++; |
| | 1094 | 121 | | previousSign = sign; |
| | | 122 | | } |
| | | 123 | | |
| | 226 | 124 | | return variations; |
| | | 125 | | } |
| | | 126 | | |
| | | 127 | | private static bool TryGetRoundedCylinderRimDistanceInterval( |
| | | 128 | | RoundedCylinderTorusPolynomial polynomial, |
| | | 129 | | Fixed64 segmentLength, |
| | | 130 | | bool startContained, |
| | | 131 | | bool endContained, |
| | | 132 | | bool calculateExit, |
| | | 133 | | out Fixed64 entry, |
| | | 134 | | out Fixed64 exit) |
| | | 135 | | { |
| | 48 | 136 | | entry = default; |
| | 48 | 137 | | exit = default; |
| | 48 | 138 | | if (segmentLength == Fixed64.Zero) |
| | 3 | 139 | | return startContained; |
| | | 140 | | |
| | 45 | 141 | | Span<ulong> coefficients = stackalloc ulong[ |
| | 45 | 142 | | RoundedCylinderPolynomialCount |
| | 45 | 143 | | * RoundedCylinderCoefficientCount |
| | 45 | 144 | | * RoundedCylinderWideLimbCount]; |
| | 45 | 145 | | Span<sbyte> signs = stackalloc sbyte[ |
| | 45 | 146 | | RoundedCylinderPolynomialCount * RoundedCylinderCoefficientCount]; |
| | 45 | 147 | | Span<int> degrees = stackalloc int[RoundedCylinderPolynomialCount]; |
| | 45 | 148 | | coefficients.Clear(); |
| | 45 | 149 | | signs.Clear(); |
| | 45 | 150 | | degrees.Clear(); |
| | | 151 | | |
| | 45 | 152 | | ulong maximum = unchecked((ulong)segmentLength.m_rawValue); |
| | 45 | 153 | | int sequenceCount = BuildRoundedCylinderSturmSequence( |
| | 45 | 154 | | polynomial, |
| | 45 | 155 | | coefficients, |
| | 45 | 156 | | signs, |
| | 45 | 157 | | degrees); |
| | 45 | 158 | | if (sequenceCount == 0) |
| | | 159 | | { |
| | 1 | 160 | | exit = segmentLength; |
| | 1 | 161 | | return true; |
| | | 162 | | } |
| | 44 | 163 | | HomogenizeRoundedCylinderSturmSequence( |
| | 44 | 164 | | coefficients, |
| | 44 | 165 | | degrees, |
| | 44 | 166 | | sequenceCount, |
| | 44 | 167 | | maximum); |
| | | 168 | | |
| | 44 | 169 | | int lowerVariations = GetRoundedCylinderSignVariations( |
| | 44 | 170 | | coefficients, |
| | 44 | 171 | | signs, |
| | 44 | 172 | | degrees, |
| | 44 | 173 | | sequenceCount, |
| | 44 | 174 | | 0UL, |
| | 44 | 175 | | 1UL); |
| | 44 | 176 | | int upperVariations = GetRoundedCylinderSignVariations( |
| | 44 | 177 | | coefficients, |
| | 44 | 178 | | signs, |
| | 44 | 179 | | degrees, |
| | 44 | 180 | | sequenceCount, |
| | 44 | 181 | | maximum, |
| | 44 | 182 | | 1UL); |
| | 44 | 183 | | int rootCount = lowerVariations - upperVariations; |
| | 44 | 184 | | bool rootAtStart = polynomial.H0.IsZero; |
| | 44 | 185 | | if (rootCount == 0) |
| | | 186 | | { |
| | 7 | 187 | | if (!startContained && !rootAtStart) |
| | 4 | 188 | | return false; |
| | | 189 | | |
| | 3 | 190 | | entry = Fixed64.Zero; |
| | 3 | 191 | | exit = endContained ? segmentLength : Fixed64.Zero; |
| | 3 | 192 | | return true; |
| | | 193 | | } |
| | | 194 | | |
| | 37 | 195 | | entry = startContained || rootAtStart |
| | 37 | 196 | | ? Fixed64.Zero |
| | 37 | 197 | | : FindRoundedCylinderRoot( |
| | 37 | 198 | | coefficients, |
| | 37 | 199 | | signs, |
| | 37 | 200 | | degrees, |
| | 37 | 201 | | sequenceCount, |
| | 37 | 202 | | maximum, |
| | 37 | 203 | | lowerVariations, |
| | 37 | 204 | | upperVariations, |
| | 37 | 205 | | rootCount, |
| | 37 | 206 | | findFirst: true); |
| | 37 | 207 | | if (!calculateExit) |
| | | 208 | | { |
| | 8 | 209 | | exit = entry; |
| | 8 | 210 | | return true; |
| | | 211 | | } |
| | 29 | 212 | | exit = endContained |
| | 29 | 213 | | ? segmentLength |
| | 29 | 214 | | : FindRoundedCylinderRoot( |
| | 29 | 215 | | coefficients, |
| | 29 | 216 | | signs, |
| | 29 | 217 | | degrees, |
| | 29 | 218 | | sequenceCount, |
| | 29 | 219 | | maximum, |
| | 29 | 220 | | lowerVariations, |
| | 29 | 221 | | upperVariations, |
| | 29 | 222 | | rootCount, |
| | 29 | 223 | | findFirst: false); |
| | 29 | 224 | | return true; |
| | | 225 | | } |
| | | 226 | | |
| | | 227 | | private static Fixed64 FindRoundedCylinderRoot( |
| | | 228 | | Span<ulong> coefficients, |
| | | 229 | | Span<sbyte> signs, |
| | | 230 | | Span<int> degrees, |
| | | 231 | | int sequenceCount, |
| | | 232 | | ulong maximum, |
| | | 233 | | int lowerVariations, |
| | | 234 | | int upperVariations, |
| | | 235 | | int totalRootCount, |
| | | 236 | | bool findFirst) |
| | | 237 | | { |
| | 52 | 238 | | ulong lower = 0UL; |
| | 52 | 239 | | ulong upper = maximum; |
| | 52 | 240 | | int intervalRootCount = totalRootCount; |
| | 742 | 241 | | while (upper - lower > 1UL) |
| | | 242 | | { |
| | 721 | 243 | | if (intervalRootCount == 1) |
| | | 244 | | { |
| | 593 | 245 | | int lowerSign = EvaluateRoundedCylinderPolynomialSign( |
| | 593 | 246 | | coefficients, |
| | 593 | 247 | | signs, |
| | 593 | 248 | | degrees[0], |
| | 593 | 249 | | 0, |
| | 593 | 250 | | lower, |
| | 593 | 251 | | 1UL); |
| | | 252 | | // The global start root is handled before isolation; midpoint |
| | | 253 | | // roots remain on the upper boundary of the retained interval. |
| | 593 | 254 | | int upperSign = EvaluateRoundedCylinderPolynomialSign( |
| | 593 | 255 | | coefficients, |
| | 593 | 256 | | signs, |
| | 593 | 257 | | degrees[0], |
| | 593 | 258 | | 0, |
| | 593 | 259 | | upper, |
| | 593 | 260 | | 1UL); |
| | 593 | 261 | | if (upperSign == 0) |
| | 4 | 262 | | return Fixed64.FromRaw((long)upper); |
| | | 263 | | |
| | 589 | 264 | | if (lowerSign != upperSign) |
| | | 265 | | { |
| | 27 | 266 | | return FindRoundedCylinderSignChangingRoot( |
| | 27 | 267 | | coefficients, |
| | 27 | 268 | | signs, |
| | 27 | 269 | | degrees[0], |
| | 27 | 270 | | maximum, |
| | 27 | 271 | | lower, |
| | 27 | 272 | | upper, |
| | 27 | 273 | | lowerSign); |
| | | 274 | | } |
| | | 275 | | } |
| | | 276 | | |
| | 690 | 277 | | ulong midpoint = lower + ((upper - lower) >> 1); |
| | 690 | 278 | | int midpointVariations = GetRoundedCylinderSignVariations( |
| | 690 | 279 | | coefficients, |
| | 690 | 280 | | signs, |
| | 690 | 281 | | degrees, |
| | 690 | 282 | | sequenceCount, |
| | 690 | 283 | | midpoint, |
| | 690 | 284 | | 1UL); |
| | 690 | 285 | | int leftRootCount = lowerVariations - midpointVariations; |
| | 690 | 286 | | if (findFirst ? leftRootCount > 0 : leftRootCount == intervalRootCount) |
| | | 287 | | { |
| | 341 | 288 | | upper = midpoint; |
| | 341 | 289 | | upperVariations = midpointVariations; |
| | 341 | 290 | | intervalRootCount = leftRootCount; |
| | | 291 | | } |
| | | 292 | | else |
| | | 293 | | { |
| | 349 | 294 | | lower = midpoint; |
| | 349 | 295 | | lowerVariations = midpointVariations; |
| | 349 | 296 | | intervalRootCount -= leftRootCount; |
| | | 297 | | } |
| | | 298 | | } |
| | | 299 | | |
| | 21 | 300 | | ulong midpointNumerator = (lower << 1) | 1UL; |
| | 21 | 301 | | int halfVariations = GetRoundedCylinderSignVariations( |
| | 21 | 302 | | coefficients, |
| | 21 | 303 | | signs, |
| | 21 | 304 | | degrees, |
| | 21 | 305 | | sequenceCount, |
| | 21 | 306 | | midpointNumerator, |
| | 21 | 307 | | 2UL); |
| | 21 | 308 | | int leftHalfRootCount = lowerVariations - halfVariations; |
| | 21 | 309 | | int midpointSign = EvaluateRoundedCylinderPolynomialSign( |
| | 21 | 310 | | coefficients, |
| | 21 | 311 | | signs, |
| | 21 | 312 | | degrees[0], |
| | 21 | 313 | | 0, |
| | 21 | 314 | | midpointNumerator, |
| | 21 | 315 | | 2UL); |
| | | 316 | | |
| | | 317 | | ulong rounded; |
| | 21 | 318 | | if (findFirst) |
| | | 319 | | { |
| | 14 | 320 | | rounded = leftHalfRootCount > 0 ? lower : upper; |
| | 14 | 321 | | if ((midpointSign == 0) & (leftHalfRootCount == 1)) |
| | 1 | 322 | | rounded = lower + (lower & 1UL); |
| | | 323 | | } |
| | | 324 | | else |
| | | 325 | | { |
| | 7 | 326 | | int rightHalfRootCount = intervalRootCount - leftHalfRootCount; |
| | 7 | 327 | | rounded = rightHalfRootCount > 0 ? upper : lower; |
| | 7 | 328 | | if ((midpointSign == 0) & (rightHalfRootCount == 0)) |
| | 1 | 329 | | rounded = lower + (lower & 1UL); |
| | | 330 | | } |
| | | 331 | | |
| | 21 | 332 | | return Fixed64.FromRaw((long)rounded); |
| | | 333 | | } |
| | | 334 | | |
| | | 335 | | private static Fixed64 FindRoundedCylinderSignChangingRoot( |
| | | 336 | | Span<ulong> coefficients, |
| | | 337 | | Span<sbyte> signs, |
| | | 338 | | int degree, |
| | | 339 | | ulong maximum, |
| | | 340 | | ulong lower, |
| | | 341 | | ulong upper, |
| | | 342 | | int lowerSign) |
| | | 343 | | { |
| | 869 | 344 | | while (upper - lower > 1UL) |
| | | 345 | | { |
| | 847 | 346 | | ulong midpoint = lower + ((upper - lower) >> 1); |
| | 847 | 347 | | int midpointSign = EvaluateRoundedCylinderPolynomialSign( |
| | 847 | 348 | | coefficients, |
| | 847 | 349 | | signs, |
| | 847 | 350 | | degree, |
| | 847 | 351 | | 0, |
| | 847 | 352 | | midpoint, |
| | 847 | 353 | | 1UL); |
| | 847 | 354 | | if (midpointSign == 0) |
| | 5 | 355 | | return Fixed64.FromRaw((long)midpoint); |
| | 842 | 356 | | if (midpointSign == lowerSign) |
| | 445 | 357 | | lower = midpoint; |
| | | 358 | | else |
| | 397 | 359 | | upper = midpoint; |
| | | 360 | | } |
| | | 361 | | |
| | 22 | 362 | | ulong midpointNumerator = (lower << 1) | 1UL; |
| | 22 | 363 | | int halfSign = EvaluateRoundedCylinderPolynomialSign( |
| | 22 | 364 | | coefficients, |
| | 22 | 365 | | signs, |
| | 22 | 366 | | degree, |
| | 22 | 367 | | 0, |
| | 22 | 368 | | midpointNumerator, |
| | 22 | 369 | | 2UL); |
| | 22 | 370 | | if (halfSign == 0) |
| | 2 | 371 | | return Fixed64.FromRaw((long)(lower + (lower & 1UL))); |
| | 20 | 372 | | return Fixed64.FromRaw((long)(halfSign == lowerSign ? upper : lower)); |
| | | 373 | | } |
| | | 374 | | |
| | | 375 | | private static void HomogenizeRoundedCylinderSturmSequence( |
| | | 376 | | Span<ulong> coefficients, |
| | | 377 | | Span<int> degrees, |
| | | 378 | | int sequenceCount, |
| | | 379 | | ulong segmentLengthRaw) |
| | | 380 | | { |
| | 44 | 381 | | Span<ulong> product = stackalloc ulong[RoundedCylinderWideLimbCount]; |
| | 420 | 382 | | for (int polynomialIndex = 0; polynomialIndex < sequenceCount; polynomialIndex++) |
| | | 383 | | { |
| | 166 | 384 | | int degree = degrees[polynomialIndex]; |
| | 1098 | 385 | | for (int coefficientIndex = 0; coefficientIndex < degree; coefficientIndex++) |
| | | 386 | | { |
| | 383 | 387 | | Span<ulong> coefficient = GetRoundedCylinderCoefficient( |
| | 383 | 388 | | coefficients, |
| | 383 | 389 | | polynomialIndex, |
| | 383 | 390 | | coefficientIndex); |
| | 383 | 391 | | int power = degree - coefficientIndex; |
| | 2350 | 392 | | for (int factor = 0; factor < power; factor++) |
| | | 393 | | { |
| | 792 | 394 | | MultiplyRoundedCylinderWideByWord(coefficient, segmentLengthRaw, product); |
| | 792 | 395 | | CopyRoundedCylinderWide(product, coefficient); |
| | | 396 | | } |
| | | 397 | | } |
| | | 398 | | } |
| | 44 | 399 | | } |
| | | 400 | | |
| | | 401 | | private static bool IsRoundedCylinderEntryBeforeOrEqualToExit( |
| | | 402 | | RoundedCylinderQuadraticBound entry, |
| | | 403 | | RoundedCylinderQuadraticBound exit) |
| | | 404 | | { |
| | 25 | 405 | | Signed576 cross = WideArithmetic.SubtractSigned576( |
| | 25 | 406 | | WideArithmetic.MultiplySigned320(entry.Numerator, exit.Denominator), |
| | 25 | 407 | | WideArithmetic.MultiplySigned320(exit.Numerator, entry.Denominator)); |
| | 25 | 408 | | bool hasEntryRadical = entry.RadicalSign < 0; |
| | 25 | 409 | | bool hasExitRadical = exit.RadicalSign > 0; |
| | 25 | 410 | | if (!hasEntryRadical && !hasExitRadical) |
| | 12 | 411 | | return cross.Sign <= 0; |
| | 13 | 412 | | if (cross.Sign <= 0) |
| | 5 | 413 | | return true; |
| | | 414 | | |
| | 8 | 415 | | Span<ulong> crossMagnitude = stackalloc ulong[RoundedCylinderWideLimbCount]; |
| | 8 | 416 | | Span<ulong> entryTerm = stackalloc ulong[RoundedCylinderWideLimbCount]; |
| | 8 | 417 | | Span<ulong> exitTerm = stackalloc ulong[RoundedCylinderWideLimbCount]; |
| | 8 | 418 | | Span<ulong> crossSquared = stackalloc ulong[RoundedCylinderWideLimbCount]; |
| | 8 | 419 | | Span<ulong> termSum = stackalloc ulong[RoundedCylinderWideLimbCount]; |
| | 8 | 420 | | Span<ulong> difference = stackalloc ulong[RoundedCylinderWideLimbCount]; |
| | 8 | 421 | | Span<ulong> differenceSquared = stackalloc ulong[RoundedCylinderWideLimbCount]; |
| | 8 | 422 | | Span<ulong> termProduct = stackalloc ulong[RoundedCylinderWideLimbCount]; |
| | 8 | 423 | | Span<ulong> right = stackalloc ulong[RoundedCylinderWideLimbCount]; |
| | 8 | 424 | | Span<ulong> first = stackalloc ulong[RoundedCylinderWideLimbCount]; |
| | 8 | 425 | | Span<ulong> second = stackalloc ulong[RoundedCylinderWideLimbCount]; |
| | 8 | 426 | | crossMagnitude.Clear(); |
| | 8 | 427 | | entryTerm.Clear(); |
| | 8 | 428 | | exitTerm.Clear(); |
| | 8 | 429 | | ImportRoundedCylinderWide(cross, crossMagnitude); |
| | 8 | 430 | | if (hasEntryRadical) |
| | | 431 | | { |
| | 5 | 432 | | GetRoundedCylinderRadicalSquare( |
| | 5 | 433 | | exit.Denominator, |
| | 5 | 434 | | entry.Discriminant, |
| | 5 | 435 | | entryTerm, |
| | 5 | 436 | | first, |
| | 5 | 437 | | second); |
| | | 438 | | } |
| | 8 | 439 | | if (hasExitRadical) |
| | | 440 | | { |
| | 6 | 441 | | GetRoundedCylinderRadicalSquare( |
| | 6 | 442 | | entry.Denominator, |
| | 6 | 443 | | exit.Discriminant, |
| | 6 | 444 | | exitTerm, |
| | 6 | 445 | | first, |
| | 6 | 446 | | second); |
| | | 447 | | } |
| | | 448 | | |
| | 8 | 449 | | MultiplyRoundedCylinderWide(crossMagnitude, crossMagnitude, crossSquared); |
| | 8 | 450 | | AddRoundedCylinderSigned( |
| | 8 | 451 | | entryTerm, |
| | 8 | 452 | | GetRoundedCylinderWideLength(entryTerm) == 0 ? (sbyte)0 : (sbyte)1, |
| | 8 | 453 | | exitTerm, |
| | 8 | 454 | | GetRoundedCylinderWideLength(exitTerm) == 0 ? (sbyte)0 : (sbyte)1, |
| | 8 | 455 | | termSum, |
| | 8 | 456 | | out _); |
| | 8 | 457 | | SubtractRoundedCylinderSigned( |
| | 8 | 458 | | crossSquared, |
| | 8 | 459 | | 1, |
| | 8 | 460 | | termSum, |
| | 8 | 461 | | 1, |
| | 8 | 462 | | difference, |
| | 8 | 463 | | out sbyte differenceSign); |
| | 8 | 464 | | if (differenceSign <= 0) |
| | 7 | 465 | | return true; |
| | | 466 | | |
| | 1 | 467 | | MultiplyRoundedCylinderWide(difference, difference, differenceSquared); |
| | 1 | 468 | | MultiplyRoundedCylinderWide(entryTerm, exitTerm, termProduct); |
| | 1 | 469 | | MultiplyRoundedCylinderWideByWord(termProduct, 4UL, right); |
| | 1 | 470 | | return WideArithmetic.CompareMagnitudeEqualLength(differenceSquared, right) <= 0; |
| | | 471 | | } |
| | | 472 | | |
| | | 473 | | private static void GetRoundedCylinderRadicalSquare( |
| | | 474 | | Signed320 coefficient, |
| | | 475 | | Signed576 discriminant, |
| | | 476 | | Span<ulong> destination, |
| | | 477 | | Span<ulong> first, |
| | | 478 | | Span<ulong> second) |
| | | 479 | | { |
| | 11 | 480 | | ImportRoundedCylinderWide(coefficient, first); |
| | 11 | 481 | | MultiplyRoundedCylinderWide(first, first, second); |
| | 11 | 482 | | ImportRoundedCylinderWide(discriminant, first); |
| | 11 | 483 | | MultiplyRoundedCylinderWide(second, first, destination); |
| | 11 | 484 | | } |
| | | 485 | | |
| | | 486 | | private static void ImportRoundedCylinderWide(Signed320 value, Span<ulong> destination) |
| | | 487 | | { |
| | 11 | 488 | | destination.Clear(); |
| | 11 | 489 | | WideArithmetic.GetMagnitude( |
| | 11 | 490 | | value, |
| | 11 | 491 | | out destination[4], |
| | 11 | 492 | | out destination[3], |
| | 11 | 493 | | out destination[2], |
| | 11 | 494 | | out destination[1], |
| | 11 | 495 | | out destination[0]); |
| | 11 | 496 | | } |
| | | 497 | | |
| | | 498 | | private static void ImportRoundedCylinderWide(Signed576 value, Span<ulong> destination) |
| | | 499 | | { |
| | 19 | 500 | | destination.Clear(); |
| | 19 | 501 | | WideArithmetic.GetMagnitude(value, destination[..9]); |
| | 19 | 502 | | } |
| | | 503 | | |
| | | 504 | | private static int GetRoundedCylinderSignVariations( |
| | | 505 | | Span<ulong> coefficients, |
| | | 506 | | Span<sbyte> signs, |
| | | 507 | | Span<int> degrees, |
| | | 508 | | int sequenceCount, |
| | | 509 | | ulong numerator, |
| | | 510 | | ulong denominator) |
| | | 511 | | { |
| | 799 | 512 | | int variations = 0; |
| | 799 | 513 | | int previousSign = 0; |
| | 6344 | 514 | | for (int polynomialIndex = 0; polynomialIndex < sequenceCount; polynomialIndex++) |
| | | 515 | | { |
| | 2373 | 516 | | int sign = EvaluateRoundedCylinderPolynomialSign( |
| | 2373 | 517 | | coefficients, |
| | 2373 | 518 | | signs, |
| | 2373 | 519 | | degrees[polynomialIndex], |
| | 2373 | 520 | | polynomialIndex, |
| | 2373 | 521 | | numerator, |
| | 2373 | 522 | | denominator); |
| | 2373 | 523 | | if (sign == 0) |
| | | 524 | | continue; |
| | 2289 | 525 | | if (previousSign != 0 && sign != previousSign) |
| | 930 | 526 | | variations++; |
| | 2289 | 527 | | previousSign = sign; |
| | | 528 | | } |
| | | 529 | | |
| | 799 | 530 | | return variations; |
| | | 531 | | } |
| | | 532 | | |
| | | 533 | | private static int EvaluateRoundedCylinderPolynomialSign( |
| | | 534 | | Span<ulong> coefficients, |
| | | 535 | | Span<sbyte> signs, |
| | | 536 | | int degree, |
| | | 537 | | int polynomialIndex, |
| | | 538 | | ulong numerator, |
| | | 539 | | ulong denominator) |
| | | 540 | | { |
| | 4449 | 541 | | if (denominator == 1UL) |
| | | 542 | | { |
| | 4341 | 543 | | return EvaluateRoundedCylinderPolynomialSignAtInteger( |
| | 4341 | 544 | | coefficients, |
| | 4341 | 545 | | signs, |
| | 4341 | 546 | | degree, |
| | 4341 | 547 | | polynomialIndex, |
| | 4341 | 548 | | numerator); |
| | | 549 | | } |
| | | 550 | | |
| | 108 | 551 | | Span<ulong> result = stackalloc ulong[RoundedCylinderWideLimbCount]; |
| | 108 | 552 | | Span<ulong> scale = stackalloc ulong[RoundedCylinderWideLimbCount]; |
| | 108 | 553 | | Span<ulong> first = stackalloc ulong[RoundedCylinderWideLimbCount]; |
| | 108 | 554 | | Span<ulong> second = stackalloc ulong[RoundedCylinderWideLimbCount]; |
| | 108 | 555 | | Span<ulong> sum = stackalloc ulong[RoundedCylinderWideLimbCount]; |
| | 108 | 556 | | result.Clear(); |
| | 108 | 557 | | scale.Clear(); |
| | 108 | 558 | | CopyRoundedCylinderWide( |
| | 108 | 559 | | GetRoundedCylinderCoefficient(coefficients, polynomialIndex, degree), |
| | 108 | 560 | | result); |
| | 108 | 561 | | sbyte resultSign = signs[GetRoundedCylinderCoefficientIndex(polynomialIndex, degree)]; |
| | 108 | 562 | | scale[0] = denominator; |
| | | 563 | | |
| | 926 | 564 | | for (int coefficientIndex = degree - 1; coefficientIndex >= 0; coefficientIndex--) |
| | | 565 | | { |
| | 355 | 566 | | MultiplyRoundedCylinderWideByWord(result, numerator, first); |
| | 355 | 567 | | sbyte firstSign = resultSign; |
| | 355 | 568 | | MultiplyRoundedCylinderWide( |
| | 355 | 569 | | GetRoundedCylinderCoefficient(coefficients, polynomialIndex, coefficientIndex), |
| | 355 | 570 | | scale, |
| | 355 | 571 | | second); |
| | 355 | 572 | | sbyte secondSign = signs[GetRoundedCylinderCoefficientIndex(polynomialIndex, coefficientIndex)]; |
| | 355 | 573 | | AddRoundedCylinderSigned( |
| | 355 | 574 | | first, |
| | 355 | 575 | | firstSign, |
| | 355 | 576 | | second, |
| | 355 | 577 | | secondSign, |
| | 355 | 578 | | sum, |
| | 355 | 579 | | out resultSign); |
| | 355 | 580 | | CopyRoundedCylinderWide(sum, result); |
| | 355 | 581 | | MultiplyRoundedCylinderWideByWord(scale, denominator, first); |
| | 355 | 582 | | CopyRoundedCylinderWide(first, scale); |
| | | 583 | | } |
| | | 584 | | |
| | 108 | 585 | | return resultSign; |
| | | 586 | | } |
| | | 587 | | |
| | | 588 | | private static int EvaluateRoundedCylinderPolynomialSignAtInteger( |
| | | 589 | | Span<ulong> coefficients, |
| | | 590 | | Span<sbyte> signs, |
| | | 591 | | int degree, |
| | | 592 | | int polynomialIndex, |
| | | 593 | | ulong parameter) |
| | | 594 | | { |
| | 4341 | 595 | | Span<ulong> result = stackalloc ulong[RoundedCylinderWideLimbCount]; |
| | 4341 | 596 | | Span<ulong> product = stackalloc ulong[RoundedCylinderWideLimbCount]; |
| | 4341 | 597 | | Span<ulong> sum = stackalloc ulong[RoundedCylinderWideLimbCount]; |
| | 4341 | 598 | | result.Clear(); |
| | 4341 | 599 | | CopyRoundedCylinderWide( |
| | 4341 | 600 | | GetRoundedCylinderCoefficient(coefficients, polynomialIndex, degree), |
| | 4341 | 601 | | result); |
| | 4341 | 602 | | sbyte resultSign = signs[GetRoundedCylinderCoefficientIndex(polynomialIndex, degree)]; |
| | | 603 | | |
| | 37894 | 604 | | for (int coefficientIndex = degree - 1; coefficientIndex >= 0; coefficientIndex--) |
| | | 605 | | { |
| | 14606 | 606 | | MultiplyRoundedCylinderWideByWord(result, parameter, product); |
| | 14606 | 607 | | AddRoundedCylinderSigned( |
| | 14606 | 608 | | product, |
| | 14606 | 609 | | resultSign, |
| | 14606 | 610 | | GetRoundedCylinderCoefficient(coefficients, polynomialIndex, coefficientIndex), |
| | 14606 | 611 | | signs[GetRoundedCylinderCoefficientIndex(polynomialIndex, coefficientIndex)], |
| | 14606 | 612 | | sum, |
| | 14606 | 613 | | out resultSign); |
| | 14606 | 614 | | CopyRoundedCylinderWide(sum, result); |
| | | 615 | | } |
| | | 616 | | |
| | 4341 | 617 | | return resultSign; |
| | | 618 | | } |
| | | 619 | | |
| | | 620 | | private static int BuildRoundedCylinderSturmSequence( |
| | | 621 | | RoundedCylinderTorusPolynomial polynomial, |
| | | 622 | | Span<ulong> coefficients, |
| | | 623 | | Span<sbyte> signs, |
| | | 624 | | Span<int> degrees) |
| | | 625 | | { |
| | 45 | 626 | | ImportRoundedCylinderCoefficient(polynomial.H0, coefficients, signs, 0, 0); |
| | 45 | 627 | | ImportRoundedCylinderCoefficient(polynomial.H1, coefficients, signs, 0, 1); |
| | 45 | 628 | | ImportRoundedCylinderCoefficient(polynomial.H2, coefficients, signs, 0, 2); |
| | 45 | 629 | | ImportRoundedCylinderCoefficient(polynomial.H3, coefficients, signs, 0, 3); |
| | 45 | 630 | | ImportRoundedCylinderCoefficient(polynomial.H4, coefficients, signs, 0, 4); |
| | 45 | 631 | | degrees[0] = TrimRoundedCylinderPolynomial(signs, 0, 4); |
| | 45 | 632 | | if (degrees[0] < 0) |
| | 1 | 633 | | return 0; |
| | 44 | 634 | | if (degrees[0] == 0) |
| | 2 | 635 | | return 1; |
| | | 636 | | |
| | 420 | 637 | | for (int coefficientIndex = 1; coefficientIndex <= degrees[0]; coefficientIndex++) |
| | | 638 | | { |
| | 168 | 639 | | int destinationIndex = coefficientIndex - 1; |
| | 168 | 640 | | MultiplyRoundedCylinderWideByWord( |
| | 168 | 641 | | GetRoundedCylinderCoefficient(coefficients, 0, coefficientIndex), |
| | 168 | 642 | | (ulong)coefficientIndex, |
| | 168 | 643 | | GetRoundedCylinderCoefficient(coefficients, 1, destinationIndex)); |
| | 168 | 644 | | signs[GetRoundedCylinderCoefficientIndex(1, destinationIndex)] = |
| | 168 | 645 | | signs[GetRoundedCylinderCoefficientIndex(0, coefficientIndex)]; |
| | | 646 | | } |
| | 42 | 647 | | degrees[1] = degrees[0] - 1; |
| | | 648 | | |
| | | 649 | | // A nonstationary authored chord makes H4 strictly positive. A |
| | | 650 | | // stationary chord was already reduced to the constant cases above. |
| | 42 | 651 | | return BuildRoundedCylinderQuarticSturmSequence(coefficients, signs, degrees); |
| | | 652 | | } |
| | | 653 | | |
| | | 654 | | private static int BuildRoundedCylinderQuarticSturmSequence( |
| | | 655 | | Span<ulong> coefficients, |
| | | 656 | | Span<sbyte> signs, |
| | | 657 | | Span<int> degrees) |
| | | 658 | | { |
| | 155 | 659 | | Span<ulong> first = stackalloc ulong[RoundedCylinderWideLimbCount]; |
| | 155 | 660 | | Span<ulong> second = stackalloc ulong[RoundedCylinderWideLimbCount]; |
| | 155 | 661 | | Span<ulong> scratch = stackalloc ulong[RoundedCylinderWideLimbCount]; |
| | 155 | 662 | | Span<ulong> result = stackalloc ulong[RoundedCylinderWideLimbCount]; |
| | | 663 | | |
| | 155 | 664 | | ReadOnlySpan<ulong> a = GetRoundedCylinderCoefficient(coefficients, 0, 4); |
| | 155 | 665 | | ReadOnlySpan<ulong> b = GetRoundedCylinderCoefficient(coefficients, 0, 3); |
| | 155 | 666 | | ReadOnlySpan<ulong> c = GetRoundedCylinderCoefficient(coefficients, 0, 2); |
| | 155 | 667 | | ReadOnlySpan<ulong> d = GetRoundedCylinderCoefficient(coefficients, 0, 1); |
| | 155 | 668 | | ReadOnlySpan<ulong> e = GetRoundedCylinderCoefficient(coefficients, 0, 0); |
| | 155 | 669 | | sbyte aSign = signs[GetRoundedCylinderCoefficientIndex(0, 4)]; |
| | 155 | 670 | | sbyte bSign = signs[GetRoundedCylinderCoefficientIndex(0, 3)]; |
| | 155 | 671 | | sbyte cSign = signs[GetRoundedCylinderCoefficientIndex(0, 2)]; |
| | 155 | 672 | | sbyte dSign = signs[GetRoundedCylinderCoefficientIndex(0, 1)]; |
| | 155 | 673 | | sbyte eSign = signs[GetRoundedCylinderCoefficientIndex(0, 0)]; |
| | | 674 | | |
| | | 675 | | // Positive-scaled negative remainder of the quartic and its derivative. |
| | 155 | 676 | | SetRoundedCylinderScaledProduct(b, bSign, b, bSign, 3UL, first, out sbyte firstSign, scratch); |
| | 155 | 677 | | SetRoundedCylinderScaledProduct(a, aSign, c, cSign, 8UL, second, out sbyte secondSign, scratch); |
| | 155 | 678 | | Span<ulong> u = GetRoundedCylinderCoefficient(coefficients, 2, 2); |
| | 155 | 679 | | SubtractRoundedCylinderSigned(first, firstSign, second, secondSign, u, out sbyte uSign); |
| | 155 | 680 | | signs[GetRoundedCylinderCoefficientIndex(2, 2)] = uSign; |
| | | 681 | | |
| | 155 | 682 | | SetRoundedCylinderProduct(b, bSign, c, cSign, first, out firstSign); |
| | 155 | 683 | | SetRoundedCylinderScaledProduct(a, aSign, d, dSign, 6UL, second, out secondSign, scratch); |
| | 155 | 684 | | SubtractRoundedCylinderSigned(first, firstSign, second, secondSign, result, out sbyte resultSign); |
| | 155 | 685 | | Span<ulong> v = GetRoundedCylinderCoefficient(coefficients, 2, 1); |
| | 155 | 686 | | MultiplyRoundedCylinderWideByWord(result, 2UL, v); |
| | 155 | 687 | | sbyte vSign = resultSign; |
| | 155 | 688 | | signs[GetRoundedCylinderCoefficientIndex(2, 1)] = vSign; |
| | | 689 | | |
| | 155 | 690 | | SetRoundedCylinderProduct(b, bSign, d, dSign, first, out firstSign); |
| | 155 | 691 | | SetRoundedCylinderScaledProduct(a, aSign, e, eSign, 16UL, second, out secondSign, scratch); |
| | 155 | 692 | | Span<ulong> w = GetRoundedCylinderCoefficient(coefficients, 2, 0); |
| | 155 | 693 | | SubtractRoundedCylinderSigned(first, firstSign, second, secondSign, w, out sbyte wSign); |
| | 155 | 694 | | signs[GetRoundedCylinderCoefficientIndex(2, 0)] = wSign; |
| | 155 | 695 | | degrees[2] = TrimRoundedCylinderPolynomial(signs, 2, 2); |
| | 155 | 696 | | if (degrees[2] < 0) |
| | 6 | 697 | | return 2; |
| | 149 | 698 | | if (degrees[2] == 0) |
| | | 699 | | { |
| | 5 | 700 | | RemoveRoundedCylinderCommonPowerOfTwo(coefficients, signs, 2, 0); |
| | 5 | 701 | | return 3; |
| | | 702 | | } |
| | 144 | 703 | | if (degrees[2] == 1) |
| | | 704 | | { |
| | 3 | 705 | | RemoveRoundedCylinderCommonPowerOfTwo(coefficients, signs, 2, 1); |
| | 3 | 706 | | int finalDegree = BuildRoundedCylinderNegativeRemainder( |
| | 3 | 707 | | coefficients, |
| | 3 | 708 | | signs, |
| | 3 | 709 | | 1, |
| | 3 | 710 | | 3, |
| | 3 | 711 | | 2, |
| | 3 | 712 | | 1, |
| | 3 | 713 | | 3); |
| | 3 | 714 | | if (finalDegree < 0) |
| | 1 | 715 | | return 3; |
| | | 716 | | |
| | 2 | 717 | | degrees[3] = finalDegree; |
| | 2 | 718 | | return 4; |
| | | 719 | | } |
| | 141 | 720 | | Span<ulong> l = GetRoundedCylinderCoefficient(coefficients, 3, 1); |
| | 141 | 721 | | Span<ulong> m = GetRoundedCylinderCoefficient(coefficients, 3, 0); |
| | | 722 | | |
| | | 723 | | // This factorization divides the ordinary third Sturm polynomial by |
| | | 724 | | // positive 16a. Both two-bit shifts are algebraically exact. |
| | 141 | 725 | | SetRoundedCylinderScaledProduct(c, cSign, c, cSign, 4UL, first, out firstSign, scratch); |
| | 141 | 726 | | AddRoundedCylinderSigned(w, wSign, first, firstSign, result, out resultSign); |
| | 141 | 727 | | SetRoundedCylinderScaledProduct(b, bSign, d, dSign, 9UL, first, out firstSign, scratch); |
| | 141 | 728 | | SubtractRoundedCylinderSigned(result, resultSign, first, firstSign, l, out sbyte xSign); |
| | | 729 | | |
| | 141 | 730 | | SetRoundedCylinderProduct(c, cSign, d, dSign, first, out firstSign); |
| | 141 | 731 | | SetRoundedCylinderScaledProduct(b, bSign, e, eSign, 6UL, second, out secondSign, scratch); |
| | 141 | 732 | | SubtractRoundedCylinderSigned(first, firstSign, second, secondSign, m, out sbyte ySign); |
| | | 733 | | |
| | 141 | 734 | | SetRoundedCylinderProduct(u, uSign, l, xSign, first, out firstSign); |
| | 141 | 735 | | SetRoundedCylinderProduct(v, vSign, v, vSign, second, out secondSign); |
| | 141 | 736 | | SubtractRoundedCylinderSigned(first, firstSign, second, secondSign, result, out resultSign); |
| | 141 | 737 | | ShiftRoundedCylinderWideRight(result, 2); |
| | 141 | 738 | | CopyRoundedCylinderWide(result, l); |
| | 141 | 739 | | sbyte lSign = resultSign; |
| | 141 | 740 | | signs[GetRoundedCylinderCoefficientIndex(3, 1)] = lSign; |
| | | 741 | | |
| | 141 | 742 | | SetRoundedCylinderProduct(u, uSign, m, ySign, first, out firstSign); |
| | 141 | 743 | | MultiplyRoundedCylinderWideByWord(first, 2UL, second); |
| | 141 | 744 | | secondSign = firstSign; |
| | 141 | 745 | | SetRoundedCylinderProduct(v, vSign, w, wSign, first, out firstSign); |
| | 141 | 746 | | SubtractRoundedCylinderSigned(second, secondSign, first, firstSign, result, out resultSign); |
| | 141 | 747 | | ShiftRoundedCylinderWideRight(result, 2); |
| | 141 | 748 | | CopyRoundedCylinderWide(result, m); |
| | 141 | 749 | | sbyte mSign = resultSign; |
| | 141 | 750 | | signs[GetRoundedCylinderCoefficientIndex(3, 0)] = mSign; |
| | | 751 | | |
| | 141 | 752 | | RemoveRoundedCylinderCommonPowerOfTwo(coefficients, signs, 2, 2); |
| | 141 | 753 | | uSign = signs[GetRoundedCylinderCoefficientIndex(2, 2)]; |
| | 141 | 754 | | vSign = signs[GetRoundedCylinderCoefficientIndex(2, 1)]; |
| | 141 | 755 | | wSign = signs[GetRoundedCylinderCoefficientIndex(2, 0)]; |
| | 141 | 756 | | degrees[3] = TrimRoundedCylinderPolynomial(signs, 3, 1); |
| | 141 | 757 | | if (degrees[3] < 0) |
| | 7 | 758 | | return 3; |
| | 134 | 759 | | RemoveRoundedCylinderCommonPowerOfTwo(coefficients, signs, 3, degrees[3]); |
| | 134 | 760 | | if (degrees[3] == 0) |
| | 2 | 761 | | return 4; |
| | | 762 | | |
| | 132 | 763 | | Span<ulong> k = GetRoundedCylinderCoefficient(coefficients, 4, 0); |
| | | 764 | | // The last remainder has the quartic discriminant's sign when u != 0. |
| | | 765 | | // Factorized invariants avoid materializing its much wider magnitude. |
| | 132 | 766 | | SetRoundedCylinderProduct(c, cSign, c, cSign, first, out firstSign); |
| | 132 | 767 | | SetRoundedCylinderScaledProduct(b, bSign, d, dSign, 3UL, second, out secondSign, scratch); |
| | 132 | 768 | | SubtractRoundedCylinderSigned(first, firstSign, second, secondSign, result, out resultSign); |
| | 132 | 769 | | SetRoundedCylinderScaledProduct(a, aSign, e, eSign, 12UL, first, out firstSign, scratch); |
| | 132 | 770 | | AddRoundedCylinderSigned(result, resultSign, first, firstSign, k, out sbyte delta0Sign); |
| | | 771 | | |
| | 132 | 772 | | SetRoundedCylinderTripleProduct(c, cSign, c, cSign, c, cSign, first, out firstSign, scratch); |
| | 132 | 773 | | MultiplyRoundedCylinderWideByWord(first, 2UL, second); |
| | 132 | 774 | | secondSign = firstSign; |
| | 132 | 775 | | SetRoundedCylinderTripleProduct(b, bSign, c, cSign, d, dSign, first, out firstSign, scratch); |
| | 132 | 776 | | MultiplyRoundedCylinderWideByWord(first, 9UL, result); |
| | 132 | 777 | | SubtractRoundedCylinderSigned(second, secondSign, result, firstSign, scratch, out secondSign); |
| | 132 | 778 | | CopyRoundedCylinderWide(scratch, second); |
| | 132 | 779 | | SetRoundedCylinderTripleProduct(b, bSign, b, bSign, e, eSign, first, out firstSign, scratch); |
| | 132 | 780 | | MultiplyRoundedCylinderWideByWord(first, 27UL, result); |
| | 132 | 781 | | AddRoundedCylinderSigned(second, secondSign, result, firstSign, scratch, out secondSign); |
| | 132 | 782 | | CopyRoundedCylinderWide(scratch, second); |
| | 132 | 783 | | SetRoundedCylinderTripleProduct(a, aSign, d, dSign, d, dSign, first, out firstSign, scratch); |
| | 132 | 784 | | MultiplyRoundedCylinderWideByWord(first, 27UL, result); |
| | 132 | 785 | | AddRoundedCylinderSigned(second, secondSign, result, firstSign, scratch, out secondSign); |
| | 132 | 786 | | CopyRoundedCylinderWide(scratch, second); |
| | 132 | 787 | | SetRoundedCylinderTripleProduct(a, aSign, c, cSign, e, eSign, first, out firstSign, scratch); |
| | 132 | 788 | | MultiplyRoundedCylinderWideByWord(first, 72UL, result); |
| | 132 | 789 | | SubtractRoundedCylinderSigned(second, secondSign, result, firstSign, scratch, out secondSign); |
| | 132 | 790 | | CopyRoundedCylinderWide(scratch, second); |
| | | 791 | | |
| | 132 | 792 | | SetRoundedCylinderTripleProduct(k, delta0Sign, k, delta0Sign, k, delta0Sign, first, out firstSign, scratch); |
| | 132 | 793 | | MultiplyRoundedCylinderWideByWord(first, 4UL, result); |
| | 132 | 794 | | SetRoundedCylinderProduct(second, secondSign, second, secondSign, scratch, out sbyte delta1SquaredSign); |
| | 132 | 795 | | SubtractRoundedCylinderSigned(result, firstSign, scratch, delta1SquaredSign, first, out sbyte kSign); |
| | 132 | 796 | | k.Clear(); |
| | 132 | 797 | | if (kSign != 0) |
| | 116 | 798 | | k[0] = 1UL; |
| | 132 | 799 | | signs[GetRoundedCylinderCoefficientIndex(4, 0)] = kSign; |
| | 132 | 800 | | if (kSign == 0) |
| | 16 | 801 | | return 4; |
| | 116 | 802 | | degrees[4] = 0; |
| | 116 | 803 | | return 5; |
| | | 804 | | } |
| | | 805 | | |
| | | 806 | | private static void SetRoundedCylinderProduct( |
| | | 807 | | ReadOnlySpan<ulong> left, |
| | | 808 | | sbyte leftSign, |
| | | 809 | | ReadOnlySpan<ulong> right, |
| | | 810 | | sbyte rightSign, |
| | | 811 | | Span<ulong> destination, |
| | | 812 | | out sbyte destinationSign) |
| | | 813 | | { |
| | 1279 | 814 | | MultiplyRoundedCylinderWide(left, right, destination); |
| | 1279 | 815 | | destinationSign = MultiplyRoundedCylinderSigns(leftSign, rightSign); |
| | 1279 | 816 | | } |
| | | 817 | | |
| | | 818 | | private static void SetRoundedCylinderScaledProduct( |
| | | 819 | | ReadOnlySpan<ulong> left, |
| | | 820 | | sbyte leftSign, |
| | | 821 | | ReadOnlySpan<ulong> right, |
| | | 822 | | sbyte rightSign, |
| | | 823 | | ulong scale, |
| | | 824 | | Span<ulong> destination, |
| | | 825 | | out sbyte destinationSign, |
| | | 826 | | Span<ulong> scratch) |
| | | 827 | | { |
| | 1307 | 828 | | MultiplyRoundedCylinderWide(left, right, scratch); |
| | 1307 | 829 | | MultiplyRoundedCylinderWideByWord(scratch, scale, destination); |
| | 1307 | 830 | | destinationSign = MultiplyRoundedCylinderSigns(leftSign, rightSign); |
| | 1307 | 831 | | } |
| | | 832 | | |
| | | 833 | | private static void SetRoundedCylinderTripleProduct( |
| | | 834 | | ReadOnlySpan<ulong> first, |
| | | 835 | | sbyte firstSign, |
| | | 836 | | ReadOnlySpan<ulong> second, |
| | | 837 | | sbyte secondSign, |
| | | 838 | | ReadOnlySpan<ulong> third, |
| | | 839 | | sbyte thirdSign, |
| | | 840 | | Span<ulong> destination, |
| | | 841 | | out sbyte destinationSign, |
| | | 842 | | Span<ulong> scratch) |
| | | 843 | | { |
| | 792 | 844 | | MultiplyRoundedCylinderWide(first, second, scratch); |
| | 792 | 845 | | MultiplyRoundedCylinderWide(scratch, third, destination); |
| | 792 | 846 | | destinationSign = MultiplyRoundedCylinderSigns( |
| | 792 | 847 | | MultiplyRoundedCylinderSigns(firstSign, secondSign), |
| | 792 | 848 | | thirdSign); |
| | 792 | 849 | | } |
| | | 850 | | |
| | | 851 | | private static int BuildRoundedCylinderNegativeRemainder( |
| | | 852 | | Span<ulong> coefficients, |
| | | 853 | | Span<sbyte> signs, |
| | | 854 | | int dividendIndex, |
| | | 855 | | int dividendDegree, |
| | | 856 | | int divisorIndex, |
| | | 857 | | int divisorDegree, |
| | | 858 | | int destinationIndex) |
| | | 859 | | { |
| | 30 | 860 | | for (int coefficientIndex = 0; coefficientIndex <= dividendDegree; coefficientIndex++) |
| | | 861 | | { |
| | 12 | 862 | | CopyRoundedCylinderWide( |
| | 12 | 863 | | GetRoundedCylinderCoefficient(coefficients, dividendIndex, coefficientIndex), |
| | 12 | 864 | | GetRoundedCylinderCoefficient(coefficients, destinationIndex, coefficientIndex)); |
| | 12 | 865 | | signs[GetRoundedCylinderCoefficientIndex(destinationIndex, coefficientIndex)] = |
| | 12 | 866 | | signs[GetRoundedCylinderCoefficientIndex(dividendIndex, coefficientIndex)]; |
| | | 867 | | } |
| | | 868 | | |
| | 3 | 869 | | Span<ulong> leadingRemainder = stackalloc ulong[RoundedCylinderWideLimbCount]; |
| | 3 | 870 | | Span<ulong> leadingDivisor = stackalloc ulong[RoundedCylinderWideLimbCount]; |
| | 3 | 871 | | Span<ulong> first = stackalloc ulong[RoundedCylinderWideLimbCount]; |
| | 3 | 872 | | Span<ulong> second = stackalloc ulong[RoundedCylinderWideLimbCount]; |
| | 3 | 873 | | Span<ulong> difference = stackalloc ulong[RoundedCylinderWideLimbCount]; |
| | 3 | 874 | | CopyRoundedCylinderWide( |
| | 3 | 875 | | GetRoundedCylinderCoefficient(coefficients, divisorIndex, divisorDegree), |
| | 3 | 876 | | leadingDivisor); |
| | 3 | 877 | | sbyte leadingDivisorSign = signs[GetRoundedCylinderCoefficientIndex(divisorIndex, divisorDegree)]; |
| | 3 | 878 | | int remainderDegree = dividendDegree; |
| | 3 | 879 | | int multiplicationCount = 0; |
| | 11 | 880 | | while (remainderDegree >= divisorDegree) |
| | | 881 | | { |
| | 9 | 882 | | CopyRoundedCylinderWide( |
| | 9 | 883 | | GetRoundedCylinderCoefficient(coefficients, destinationIndex, remainderDegree), |
| | 9 | 884 | | leadingRemainder); |
| | 9 | 885 | | sbyte leadingRemainderSign = signs[ |
| | 9 | 886 | | GetRoundedCylinderCoefficientIndex(destinationIndex, remainderDegree)]; |
| | 9 | 887 | | int offset = remainderDegree - divisorDegree; |
| | | 888 | | |
| | 72 | 889 | | for (int coefficientIndex = 0; coefficientIndex <= remainderDegree; coefficientIndex++) |
| | | 890 | | { |
| | 27 | 891 | | Span<ulong> remainderCoefficient = GetRoundedCylinderCoefficient( |
| | 27 | 892 | | coefficients, |
| | 27 | 893 | | destinationIndex, |
| | 27 | 894 | | coefficientIndex); |
| | 27 | 895 | | MultiplyRoundedCylinderWide( |
| | 27 | 896 | | leadingDivisor, |
| | 27 | 897 | | remainderCoefficient, |
| | 27 | 898 | | first); |
| | 27 | 899 | | sbyte firstSign = MultiplyRoundedCylinderSigns( |
| | 27 | 900 | | leadingDivisorSign, |
| | 27 | 901 | | signs[GetRoundedCylinderCoefficientIndex(destinationIndex, coefficientIndex)]); |
| | | 902 | | |
| | 27 | 903 | | int divisorCoefficientIndex = coefficientIndex - offset; |
| | | 904 | | sbyte secondSign; |
| | 27 | 905 | | if (divisorCoefficientIndex >= 0 && divisorCoefficientIndex <= divisorDegree) |
| | | 906 | | { |
| | 18 | 907 | | MultiplyRoundedCylinderWide( |
| | 18 | 908 | | leadingRemainder, |
| | 18 | 909 | | GetRoundedCylinderCoefficient(coefficients, divisorIndex, divisorCoefficientIndex), |
| | 18 | 910 | | second); |
| | 18 | 911 | | secondSign = MultiplyRoundedCylinderSigns( |
| | 18 | 912 | | leadingRemainderSign, |
| | 18 | 913 | | signs[GetRoundedCylinderCoefficientIndex(divisorIndex, divisorCoefficientIndex)]); |
| | | 914 | | } |
| | | 915 | | else |
| | | 916 | | { |
| | 9 | 917 | | second.Clear(); |
| | 9 | 918 | | secondSign = 0; |
| | | 919 | | } |
| | | 920 | | |
| | 27 | 921 | | SubtractRoundedCylinderSigned( |
| | 27 | 922 | | first, |
| | 27 | 923 | | firstSign, |
| | 27 | 924 | | second, |
| | 27 | 925 | | secondSign, |
| | 27 | 926 | | difference, |
| | 27 | 927 | | out sbyte differenceSign); |
| | 27 | 928 | | CopyRoundedCylinderWide(difference, remainderCoefficient); |
| | 27 | 929 | | signs[GetRoundedCylinderCoefficientIndex(destinationIndex, coefficientIndex)] = differenceSign; |
| | | 930 | | } |
| | | 931 | | |
| | 9 | 932 | | multiplicationCount++; |
| | 9 | 933 | | remainderDegree = TrimRoundedCylinderPolynomial( |
| | 9 | 934 | | signs, |
| | 9 | 935 | | destinationIndex, |
| | 9 | 936 | | remainderDegree - 1); |
| | 9 | 937 | | if (remainderDegree < 0) |
| | 1 | 938 | | return -1; |
| | | 939 | | } |
| | | 940 | | |
| | 2 | 941 | | bool pseudoScaleIsPositive = leadingDivisorSign > 0 || (multiplicationCount & 1) == 0; |
| | 2 | 942 | | if (pseudoScaleIsPositive) |
| | | 943 | | { |
| | 4 | 944 | | for (int coefficientIndex = 0; coefficientIndex <= remainderDegree; coefficientIndex++) |
| | | 945 | | { |
| | 1 | 946 | | int signIndex = GetRoundedCylinderCoefficientIndex(destinationIndex, coefficientIndex); |
| | 1 | 947 | | signs[signIndex] = (sbyte)-signs[signIndex]; |
| | | 948 | | } |
| | | 949 | | } |
| | | 950 | | |
| | 2 | 951 | | RemoveRoundedCylinderCommonPowerOfTwo( |
| | 2 | 952 | | coefficients, |
| | 2 | 953 | | signs, |
| | 2 | 954 | | destinationIndex, |
| | 2 | 955 | | remainderDegree); |
| | 2 | 956 | | return remainderDegree; |
| | | 957 | | } |
| | | 958 | | |
| | | 959 | | private static void RemoveRoundedCylinderCommonPowerOfTwo( |
| | | 960 | | Span<ulong> coefficients, |
| | | 961 | | Span<sbyte> signs, |
| | | 962 | | int polynomialIndex, |
| | | 963 | | int degree) |
| | | 964 | | { |
| | 285 | 965 | | int commonShift = int.MaxValue; |
| | 1974 | 966 | | for (int coefficientIndex = 0; coefficientIndex <= degree; coefficientIndex++) |
| | | 967 | | { |
| | 702 | 968 | | if (signs[GetRoundedCylinderCoefficientIndex(polynomialIndex, coefficientIndex)] == 0) |
| | | 969 | | continue; |
| | 605 | 970 | | commonShift = Math.Min( |
| | 605 | 971 | | commonShift, |
| | 605 | 972 | | CountRoundedCylinderTrailingZeroes( |
| | 605 | 973 | | GetRoundedCylinderCoefficient(coefficients, polynomialIndex, coefficientIndex))); |
| | | 974 | | } |
| | | 975 | | |
| | 1974 | 976 | | for (int coefficientIndex = 0; coefficientIndex <= degree; coefficientIndex++) |
| | | 977 | | { |
| | 702 | 978 | | if (signs[GetRoundedCylinderCoefficientIndex(polynomialIndex, coefficientIndex)] != 0) |
| | | 979 | | { |
| | 605 | 980 | | ShiftRoundedCylinderWideRight( |
| | 605 | 981 | | GetRoundedCylinderCoefficient(coefficients, polynomialIndex, coefficientIndex), |
| | 605 | 982 | | commonShift); |
| | | 983 | | } |
| | | 984 | | } |
| | 285 | 985 | | } |
| | | 986 | | |
| | | 987 | | private static int TrimRoundedCylinderPolynomial( |
| | | 988 | | Span<sbyte> signs, |
| | | 989 | | int polynomialIndex, |
| | | 990 | | int degree) |
| | | 991 | | { |
| | 411 | 992 | | while (degree >= 0 |
| | 411 | 993 | | && signs[GetRoundedCylinderCoefficientIndex(polynomialIndex, degree)] == 0) |
| | | 994 | | { |
| | 61 | 995 | | degree--; |
| | | 996 | | } |
| | 350 | 997 | | return degree; |
| | | 998 | | } |
| | | 999 | | |
| | | 1000 | | } |