| | | 1 | | //======================================================================= |
| | | 2 | | // WideFiniteAxisIntersection.ValueRoot.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 | | /// All-root value isolation, separate from the optimized largest-root quartic |
| | | 14 | | /// contract. Signed raw subresultants preserve exact Sturm variations. |
| | | 15 | | /// </content> |
| | | 16 | | internal static partial class WideFiniteAxisIntersection |
| | | 17 | | { |
| | | 18 | | private static int GetFiniteValueSturmArenaBits(int coefficientBits, int degree) |
| | | 19 | | { |
| | | 20 | | // C=B+64. For degree <=4, Brown's retained minors use <=39C; |
| | | 21 | | // graded pseudo-remainders plus two products use <=50C and exact |
| | | 22 | | // division <=52C. The existing eight scalar slots each use <=28C, |
| | | 23 | | // hence construction fits (39+52+8*28)C=315C. Under the unchanged |
| | | 24 | | // isolation shift bound, retained rows plus dyadic evaluation need |
| | | 25 | | // <=170C+318 bits. Callers add the existing 512-word carry padding. |
| | | 26 | | // Deeper refinement separately takes its target evaluation maximum. |
| | 9895 | 27 | | return (degree <= 4 ? 315 : 880) * (coefficientBits + 64); |
| | | 28 | | } |
| | | 29 | | |
| | | 30 | | internal static int GetFiniteValueRootCellWords(ReadOnlySpan<ulong> coefficients, |
| | | 31 | | ReadOnlySpan<sbyte> signs) => |
| | 30479 | 32 | | (16 * (GetFiniteRootCoefficientBits(coefficients, signs.Length) + 11) + 255) / 64; |
| | | 33 | | |
| | | 34 | | internal static bool TryGetFiniteValueRoot(ReadOnlySpan<ulong> coefficients, |
| | | 35 | | ReadOnlySpan<sbyte> signs, int ordinal, Span<ulong> cell, out FiniteAxisValueRoot root) |
| | | 36 | | { |
| | 3859 | 37 | | root = default; |
| | 3859 | 38 | | int degree = signs.Length - 1; |
| | 9259 | 39 | | while (degree >= 0 && signs[degree] == 0) |
| | 5400 | 40 | | degree--; |
| | 3859 | 41 | | if (degree <= 0) |
| | 19 | 42 | | return false; |
| | 3840 | 43 | | int words = coefficients.Length / signs.Length; |
| | 3840 | 44 | | coefficients = coefficients[..((degree + 1) * words)]; |
| | 3840 | 45 | | signs = signs[..(degree + 1)]; |
| | 3840 | 46 | | int bits = GetFiniteRootCoefficientBits(coefficients, signs.Length); |
| | | 47 | | // Graded surviving-coefficient slots avoid the uniform-width |
| | | 48 | | // pseudo-remainder blow-up. Construction scratch is released before |
| | | 49 | | // endpoint evaluation; neither phase holds another chain. |
| | 3840 | 50 | | Span<ulong> arena = stackalloc ulong[(GetFiniteValueSturmArenaBits(bits, degree) + 63) / 64 + 512]; |
| | 3840 | 51 | | Span<int> offsets = stackalloc int[9]; |
| | 3840 | 52 | | Span<int> widths = stackalloc int[9]; |
| | 3840 | 53 | | Span<int> degrees = stackalloc int[9]; |
| | 3840 | 54 | | Span<sbyte> sturmSigns = stackalloc sbyte[81]; |
| | 3840 | 55 | | int count = BuildFiniteValueSturm(coefficients, signs, arena, |
| | 3840 | 56 | | offsets, widths, degrees, sturmSigns, out int retainedWords); |
| | 3840 | 57 | | int rootCount = GetFiniteValueRootCount(cell, arena[..retainedWords], offsets, widths, |
| | 3840 | 58 | | degrees, sturmSigns, count, arena[retainedWords..], out int lowerVariations); |
| | 3840 | 59 | | if (ordinal < 0 || ordinal >= rootCount) |
| | 914 | 60 | | return false; |
| | 2926 | 61 | | root = IsolateFiniteValueRoot(coefficients, signs, ordinal, rootCount, lowerVariations, cell, int.MaxValue, |
| | 2926 | 62 | | arena[..retainedWords], offsets, widths, degrees, sturmSigns, count, |
| | 2926 | 63 | | arena[retainedWords..], out _); |
| | 2926 | 64 | | return true; |
| | | 65 | | } |
| | | 66 | | |
| | | 67 | | internal static FiniteAxisValueRoots GetFiniteValueRoots(ReadOnlySpan<ulong> coefficients, |
| | | 68 | | ReadOnlySpan<sbyte> signs, Span<ulong> cells, Span<int> shifts) |
| | | 69 | | { |
| | 28817 | 70 | | bool compact = TryGetFiniteValueRoots(coefficients, signs, cells, shifts, |
| | 28817 | 71 | | out int count, out int rationalMask, out bool repeated); |
| | 28817 | 72 | | return new FiniteAxisValueRoots(cells, shifts, count, rationalMask, compact, repeated); |
| | | 73 | | } |
| | | 74 | | |
| | | 75 | | /// <summary> |
| | | 76 | | /// Isolates all distinct roots in (0,1] with one Sturm construction. |
| | | 77 | | /// Cells supplies eight equal, nonempty numerator slots and shifts supplies |
| | | 78 | | /// eight entries. Bit ordinal of rationalMask identifies singleton cells. |
| | | 79 | | /// Count is the total number of distinct roots on both outcomes. Returns |
| | | 80 | | /// false when compact storage is insufficient; discard all cell metadata |
| | | 81 | | /// on false and materialize valid ordinals with full-size storage. |
| | | 82 | | /// The hasRepeatedRoots flag concerns the whole polynomial, including roots |
| | | 83 | | /// outside (0,1], and stays valid on capacity failure. Zero and constant |
| | | 84 | | /// polynomials report no repeated roots. |
| | | 85 | | /// Inputs must not overlap outputs. Storage affects reuse, not semantics. |
| | | 86 | | /// </summary> |
| | | 87 | | internal static bool TryGetFiniteValueRoots(ReadOnlySpan<ulong> coefficients, |
| | | 88 | | ReadOnlySpan<sbyte> signs, Span<ulong> cells, Span<int> shifts, |
| | | 89 | | out int count, out int rationalMask, out bool hasRepeatedRoots) |
| | | 90 | | { |
| | 28868 | 91 | | count = rationalMask = 0; |
| | 28868 | 92 | | hasRepeatedRoots = false; |
| | 28868 | 93 | | cells.Clear(); |
| | 28868 | 94 | | shifts.Clear(); |
| | 28868 | 95 | | int degree = signs.Length - 1; |
| | 31530 | 96 | | while (degree >= 0 && signs[degree] == 0) |
| | 2662 | 97 | | degree--; |
| | 28868 | 98 | | if (degree <= 0) |
| | 629 | 99 | | return true; |
| | 28239 | 100 | | int inputWords = coefficients.Length / signs.Length; |
| | 28239 | 101 | | coefficients = coefficients[..((degree + 1) * inputWords)]; |
| | 28239 | 102 | | signs = signs[..(degree + 1)]; |
| | 28239 | 103 | | int cellWords = cells.Length / 8; |
| | | 104 | | System.Diagnostics.Debug.Assert(cellWords > 0 && cells.Length == 8 * cellWords && shifts.Length == 8); |
| | 28239 | 105 | | int bits = GetFiniteRootCoefficientBits(coefficients, signs.Length); |
| | 28239 | 106 | | if (TryGetFiniteValueRootsBernstein(coefficients, signs, cells, shifts, out count)) |
| | 24252 | 107 | | return true; |
| | | 108 | | // The same arena as single-root isolation; the only additional live |
| | | 109 | | // caller state is the compact output (544 bytes for eight-word cells). |
| | 3987 | 110 | | Span<ulong> arena = stackalloc ulong[(GetFiniteValueSturmArenaBits(bits, degree) + 63) / 64 + 512]; |
| | 3987 | 111 | | Span<int> offsets = stackalloc int[9]; |
| | 3987 | 112 | | Span<int> widths = stackalloc int[9]; |
| | 3987 | 113 | | Span<int> degrees = stackalloc int[9]; |
| | 3987 | 114 | | Span<sbyte> sturmSigns = stackalloc sbyte[81]; |
| | 3987 | 115 | | int chainCount = BuildFiniteValueSturm(coefficients, signs, arena, |
| | 3987 | 116 | | offsets, widths, degrees, sturmSigns, out int retainedWords); |
| | | 117 | | // The last nonzero row is gcd(F,F') up to a nonzero scalar. Its |
| | | 118 | | // degree answers this globally without normalizing or exporting it. |
| | 3987 | 119 | | hasRepeatedRoots = degrees[chainCount - 1] > 0; |
| | 3987 | 120 | | count = GetFiniteValueRootCount(cells[..cellWords], arena[..retainedWords], offsets, |
| | 3987 | 121 | | widths, degrees, sturmSigns, chainCount, arena[retainedWords..], out int lowerVariations); |
| | 3987 | 122 | | return TryIsolateFiniteValueRootCells(coefficients, signs, cells, shifts, count, lowerVariations, |
| | 3987 | 123 | | arena[..retainedWords], offsets, widths, degrees, sturmSigns, chainCount, |
| | 3987 | 124 | | arena[retainedWords..], out rationalMask); |
| | | 125 | | } |
| | | 126 | | |
| | | 127 | | private static bool TryIsolateFiniteValueRootCells(ReadOnlySpan<ulong> coefficients, |
| | | 128 | | ReadOnlySpan<sbyte> signs, Span<ulong> cells, Span<int> shifts, int rootCount, int lowerVariations, |
| | | 129 | | scoped ReadOnlySpan<ulong> chain, scoped ReadOnlySpan<int> offsets, scoped ReadOnlySpan<int> widths, |
| | | 130 | | scoped ReadOnlySpan<int> degrees, scoped ReadOnlySpan<sbyte> sturmSigns, int chainCount, |
| | | 131 | | scoped Span<ulong> evaluation, out int rationalMask) |
| | | 132 | | { |
| | 4295 | 133 | | rationalMask = 0; |
| | 4295 | 134 | | if (rootCount == 0) |
| | 1166 | 135 | | return true; |
| | 3129 | 136 | | int cellWords = cells.Length / 8; |
| | | 137 | | // A pending node owns the output slot of its first root. Continue left |
| | | 138 | | // and save only a genuine right sibling. Finalizing roots in ascending |
| | | 139 | | // order makes the next ordinal its next pending node, regardless of |
| | | 140 | | // subdivision depth. No second set of numerator cells is live. |
| | 3129 | 141 | | Span<int> remainingCounts = stackalloc int[8]; |
| | 3129 | 142 | | Span<int> lowerCounts = stackalloc int[8]; |
| | 3129 | 143 | | remainingCounts[0] = rootCount; |
| | 3129 | 144 | | lowerCounts[0] = lowerVariations; |
| | 3129 | 145 | | cells[..cellWords].Clear(); |
| | 3129 | 146 | | shifts[0] = 0; |
| | 23224 | 147 | | for (int ordinal = 0; ordinal < rootCount; ordinal++) |
| | | 148 | | { |
| | 8502 | 149 | | int remaining = remainingCounts[ordinal]; |
| | 8502 | 150 | | int lowerCount = lowerCounts[ordinal]; |
| | 8502 | 151 | | var root = new FiniteAxisValueRoot |
| | 8502 | 152 | | { |
| | 8502 | 153 | | Coefficients = coefficients, Signs = signs, Ordinal = ordinal, |
| | 8502 | 154 | | LowerNumerator = cells.Slice(ordinal * cellWords, cellWords), |
| | 8502 | 155 | | DenominatorShift = shifts[ordinal] |
| | 8502 | 156 | | }; |
| | 51923 | 157 | | while (remaining > 1) |
| | | 158 | | { |
| | 43432 | 159 | | if (root.DenominatorShift >= cellWords * 64 - 1) |
| | | 160 | | { |
| | 11 | 161 | | rationalMask = 0; |
| | 11 | 162 | | return false; |
| | | 163 | | } |
| | 43421 | 164 | | ShiftFiniteRootLeft(root.LowerNumerator, 1); |
| | 43421 | 165 | | AddRoundedCylinderWord(root.LowerNumerator, 0, 1); |
| | 43421 | 166 | | root.DenominatorShift++; |
| | 43421 | 167 | | int midpointCount = GetFiniteValueVariations(chain, offsets, widths, degrees, |
| | 43421 | 168 | | sturmSigns, chainCount, root.LowerNumerator, root.DenominatorShift, evaluation); |
| | 43421 | 169 | | int leftCount = lowerCount - midpointCount; |
| | 43421 | 170 | | if (leftCount == 0) |
| | 10711 | 171 | | lowerCount = midpointCount; |
| | | 172 | | else |
| | | 173 | | { |
| | 32710 | 174 | | if (leftCount < remaining) |
| | | 175 | | { |
| | 5373 | 176 | | int rightOrdinal = ordinal + leftCount; |
| | 5373 | 177 | | root.LowerNumerator.CopyTo(cells.Slice(rightOrdinal * cellWords, cellWords)); |
| | 5373 | 178 | | shifts[rightOrdinal] = root.DenominatorShift; |
| | 5373 | 179 | | remainingCounts[rightOrdinal] = remaining - leftCount; |
| | 5373 | 180 | | lowerCounts[rightOrdinal] = midpointCount; |
| | | 181 | | } |
| | 32710 | 182 | | root.LowerNumerator[0]--; // Odd midpoint; left interval includes it, right excludes it. |
| | 32710 | 183 | | remaining = leftCount; |
| | | 184 | | } |
| | | 185 | | } |
| | | 186 | | // A dyadic upper endpoint may have been recognized later than in |
| | | 187 | | // ordinal isolation. Restore its minimal denominator, including one. |
| | 8491 | 188 | | if (GetFiniteValueCellUpperSign(root, evaluation) == 0) |
| | | 189 | | { |
| | 1408 | 190 | | AddRoundedCylinderWord(root.LowerNumerator, 0, 1); |
| | 1408 | 191 | | int power = CountRoundedCylinderTrailingZeroes(root.LowerNumerator); |
| | 1408 | 192 | | ShiftRoundedCylinderWideRight(root.LowerNumerator, power); |
| | 1408 | 193 | | root.DenominatorShift -= power; |
| | 1408 | 194 | | root.IsRational = true; |
| | | 195 | | } |
| | | 196 | | else |
| | | 197 | | { |
| | 7083 | 198 | | IsolateFiniteValueCell(ref root, 0, 1, lowerCount, cellWords * 64 - 1, |
| | 7083 | 199 | | chain, offsets, widths, degrees, sturmSigns, chainCount, evaluation, out bool fits); |
| | 7083 | 200 | | if (!fits) |
| | | 201 | | { |
| | 8 | 202 | | rationalMask = 0; |
| | 8 | 203 | | return false; |
| | | 204 | | } |
| | | 205 | | } |
| | 8483 | 206 | | shifts[ordinal] = root.DenominatorShift; |
| | 8483 | 207 | | if (root.IsRational) |
| | 1840 | 208 | | rationalMask |= 1 << ordinal; |
| | | 209 | | } |
| | 3110 | 210 | | return true; |
| | | 211 | | } |
| | | 212 | | |
| | | 213 | | // The caller obtained rootCount for this exact polynomial and supplies a |
| | | 214 | | // valid ordinal and GetFiniteValueRootCellWords-sized storage. The shared |
| | | 215 | | // isolation bound therefore completes without a capacity-failure route. |
| | | 216 | | internal static FiniteAxisValueRoot GetKnownFiniteValueRoot(ReadOnlySpan<ulong> coefficients, |
| | | 217 | | ReadOnlySpan<sbyte> signs, int ordinal, int rootCount, Span<ulong> cell) |
| | | 218 | | { |
| | 36 | 219 | | int bits = GetFiniteRootCoefficientBits(coefficients, signs.Length); |
| | 36 | 220 | | Span<ulong> arena = stackalloc ulong[(GetFiniteValueSturmArenaBits(bits, signs.Length - 1) + 63) / 64 + 512]; |
| | 36 | 221 | | Span<int> offsets = stackalloc int[9]; |
| | 36 | 222 | | Span<int> widths = stackalloc int[9]; |
| | 36 | 223 | | Span<int> degrees = stackalloc int[9]; |
| | 36 | 224 | | Span<sbyte> sturmSigns = stackalloc sbyte[81]; |
| | 36 | 225 | | int count = BuildFiniteValueSturm(coefficients, signs, arena, |
| | 36 | 226 | | offsets, widths, degrees, sturmSigns, out int retainedWords); |
| | 36 | 227 | | cell.Clear(); |
| | 36 | 228 | | int lowerVariations = GetFiniteValueVariations(arena[..retainedWords], offsets, widths, |
| | 36 | 229 | | degrees, sturmSigns, count, cell, 0, arena[retainedWords..]); |
| | 36 | 230 | | return IsolateFiniteValueRoot(coefficients, signs, ordinal, rootCount, lowerVariations, cell, int.MaxValue, |
| | 36 | 231 | | arena[..retainedWords], offsets, widths, degrees, sturmSigns, count, |
| | 36 | 232 | | arena[retainedWords..], out _); |
| | | 233 | | } |
| | | 234 | | |
| | | 235 | | private static int GetFiniteValueRootCount(Span<ulong> cell, ReadOnlySpan<ulong> chain, |
| | | 236 | | ReadOnlySpan<int> offsets, ReadOnlySpan<int> widths, ReadOnlySpan<int> degrees, |
| | | 237 | | ReadOnlySpan<sbyte> sturmSigns, int count, Span<ulong> evaluation, out int lowerVariations) |
| | | 238 | | { |
| | 8135 | 239 | | cell.Clear(); |
| | 8135 | 240 | | lowerVariations = GetFiniteValueVariations(chain, offsets, widths, |
| | 8135 | 241 | | degrees, sturmSigns, count, cell, 0, evaluation); |
| | 8135 | 242 | | cell[0] = 1; |
| | 8135 | 243 | | return lowerVariations - GetFiniteValueVariations(chain, offsets, widths, |
| | 8135 | 244 | | degrees, sturmSigns, count, cell, 0, evaluation); |
| | | 245 | | } |
| | | 246 | | |
| | | 247 | | private static FiniteAxisValueRoot IsolateFiniteValueRoot(ReadOnlySpan<ulong> coefficients, |
| | | 248 | | ReadOnlySpan<sbyte> signs, int ordinal, int rootCount, int lowerVariations, Span<ulong> cell, int maximumShift, |
| | | 249 | | scoped ReadOnlySpan<ulong> chain, scoped ReadOnlySpan<int> offsets, |
| | | 250 | | scoped ReadOnlySpan<int> widths, scoped ReadOnlySpan<int> degrees, |
| | | 251 | | scoped ReadOnlySpan<sbyte> sturmSigns, int count, scoped Span<ulong> evaluation, |
| | | 252 | | out bool fits) |
| | | 253 | | { |
| | 2962 | 254 | | fits = true; |
| | 2962 | 255 | | cell.Clear(); |
| | 2962 | 256 | | cell[0] = 1; |
| | 2962 | 257 | | int remaining = rootCount; |
| | 2962 | 258 | | var root = new FiniteAxisValueRoot |
| | 2962 | 259 | | { |
| | 2962 | 260 | | Coefficients = coefficients, |
| | 2962 | 261 | | Signs = signs, |
| | 2962 | 262 | | LowerNumerator = cell, |
| | 2962 | 263 | | Ordinal = ordinal |
| | 2962 | 264 | | }; |
| | 2962 | 265 | | if (ordinal == remaining - 1 |
| | 2962 | 266 | | && EvaluateFiniteRootPolynomial(coefficients, signs, cell, 0, 0, evaluation) == 0) |
| | | 267 | | { |
| | 60 | 268 | | root.IsRational = true; |
| | 60 | 269 | | return root; |
| | | 270 | | } |
| | 2902 | 271 | | cell.Clear(); |
| | | 272 | | |
| | 2902 | 273 | | IsolateFiniteValueCell(ref root, ordinal, remaining, lowerVariations, maximumShift, |
| | 2902 | 274 | | chain, offsets, widths, degrees, sturmSigns, count, evaluation, out fits); |
| | 2902 | 275 | | return root; |
| | | 276 | | } |
| | | 277 | | |
| | | 278 | | private static void IsolateFiniteValueCell(ref FiniteAxisValueRoot root, |
| | | 279 | | int ordinal, int remaining, int lowerVariations, int maximumShift, |
| | | 280 | | scoped ReadOnlySpan<ulong> chain, scoped ReadOnlySpan<int> offsets, |
| | | 281 | | scoped ReadOnlySpan<int> widths, scoped ReadOnlySpan<int> degrees, |
| | | 282 | | scoped ReadOnlySpan<sbyte> sturmSigns, int count, scoped Span<ulong> evaluation, out bool fits) |
| | | 283 | | { |
| | 9985 | 284 | | fits = true; |
| | | 285 | | |
| | | 286 | | // Distinct roots of degree-eight integer polynomials of height B, |
| | | 287 | | // including their factors, have separation >2^(-16*(B+11)-89). |
| | | 288 | | // The physical value interval is (0,1], so the retained numerator |
| | | 289 | | // needs no extra Cauchy-scale integer part. Zero is excluded. |
| | 48697 | 290 | | while (remaining != 1 || GetRoundedCylinderWideLength(root.LowerNumerator) == 0 |
| | 48697 | 291 | | || EvaluateFiniteRootPolynomial(root.Coefficients, root.Signs, root.LowerNumerator, root.DenominatorShift, 0 |
| | 48697 | 292 | | || GetFiniteValueCellUpperSign(root, evaluation) == 0) |
| | | 293 | | { |
| | | 294 | | // For x in (0,1), a shift q needs at most q+1 numerator bits |
| | | 295 | | // including the upper endpoint. Stop before shifting beyond the |
| | | 296 | | // compact slot, even when a tiny numerator could fit farther. |
| | | 297 | | // Single-root callers supply the proven full isolation capacity. |
| | 39525 | 298 | | if (root.DenominatorShift >= maximumShift) |
| | | 299 | | { |
| | 8 | 300 | | fits = false; |
| | 8 | 301 | | return; |
| | | 302 | | } |
| | 39517 | 303 | | RefineFiniteValueCell(ref root, chain, offsets, widths, degrees, sturmSigns, |
| | 39517 | 304 | | count, evaluation, ref ordinal, ref remaining, ref lowerVariations); |
| | 39517 | 305 | | if (root.IsRational) |
| | 805 | 306 | | return; |
| | | 307 | | } |
| | 9172 | 308 | | } |
| | | 309 | | |
| | | 310 | | /// <summary> |
| | | 311 | | /// Refines an already isolated root without changing its ordinal or input |
| | | 312 | | /// polynomial. Caller-owned numerator storage must preserve the existing |
| | | 313 | | /// cell and provide at least targetShift+1 bits for the refined cell and |
| | | 314 | | /// its upper endpoint. Rational roots and already finer cells are unchanged. |
| | | 315 | | /// An optional knownLowerSign of -1 or +1 certifies that this same cell has |
| | | 316 | | /// opposite nonzero endpoint signs; zero requests exact endpoint evaluation. |
| | | 317 | | /// </summary> |
| | | 318 | | internal static void RefineFiniteValueRoot(ref FiniteAxisValueRoot root, int targetShift, int knownLowerSign = 0) |
| | | 319 | | { |
| | 42121 | 320 | | if (root.IsRational || root.DenominatorShift >= targetShift) |
| | 14556 | 321 | | return; |
| | 27565 | 322 | | int coefficientBits = GetFiniteRootCoefficientBits(root.Coefficients, root.Signs.Length); |
| | 27565 | 323 | | if (TryRefineFiniteValueCrossingRoot(ref root, targetShift, knownLowerSign, coefficientBits)) |
| | 27261 | 324 | | return; |
| | 304 | 325 | | int degree = root.Signs.Length - 1; |
| | 304 | 326 | | int bits = coefficientBits; |
| | 304 | 327 | | int unitBits = bits + 64; |
| | 304 | 328 | | int retainedBits = (3 * degree + 1) * unitBits; |
| | 304 | 329 | | int evaluationBits = 0; |
| | 4748 | 330 | | for (int index = 0; index <= degree; index++) |
| | | 331 | | { |
| | 2070 | 332 | | int height = index >= degree - 1 ? unitBits : (2 * degree - 2 * index - 1) * unitBits; |
| | 2070 | 333 | | if (index < degree - 1) |
| | 1462 | 334 | | retainedBits += (index + 1) * height; |
| | 2070 | 335 | | evaluationBits = Math.Max(evaluationBits, |
| | 2070 | 336 | | 2 * (height + index * (targetShift + 1) + 95) + height + 128); |
| | | 337 | | } |
| | | 338 | | // Raw subresultant minor bounds control the retained rows. A |
| | | 339 | | // refinement may need more precision than isolation, so separately |
| | | 340 | | // budget exact dyadic evaluation and reuse the construction arena. |
| | 304 | 341 | | int arenaBits = Math.Max(GetFiniteValueSturmArenaBits(bits, degree), retainedBits + evaluationBits); |
| | 304 | 342 | | Span<ulong> arena = stackalloc ulong[(arenaBits + 63) / 64 + 512]; |
| | 304 | 343 | | Span<int> offsets = stackalloc int[9]; |
| | 304 | 344 | | Span<int> widths = stackalloc int[9]; |
| | 304 | 345 | | Span<int> degrees = stackalloc int[9]; |
| | 304 | 346 | | Span<sbyte> signs = stackalloc sbyte[81]; |
| | 304 | 347 | | int count = BuildFiniteValueSturm(root.Coefficients, root.Signs, arena, |
| | 304 | 348 | | offsets, widths, degrees, signs, out int used); |
| | 304 | 349 | | Span<ulong> chain = arena[..used]; |
| | 304 | 350 | | Span<ulong> evaluation = arena[used..]; |
| | 304 | 351 | | int lowerVariations = GetFiniteValueVariations(chain, offsets, widths, |
| | 304 | 352 | | degrees, signs, count, root.LowerNumerator, root.DenominatorShift, evaluation); |
| | 304 | 353 | | int ordinal = 0; |
| | 304 | 354 | | int remaining = 1; |
| | 13293 | 355 | | while (!root.IsRational && root.DenominatorShift < targetShift) |
| | 12989 | 356 | | RefineFiniteValueCell(ref root, chain, offsets, widths, degrees, signs, |
| | 12989 | 357 | | count, evaluation, ref ordinal, ref remaining, ref lowerVariations); |
| | 304 | 358 | | } |
| | | 359 | | |
| | | 360 | | private static void RefineFiniteValueCell(ref FiniteAxisValueRoot root, |
| | | 361 | | scoped ReadOnlySpan<ulong> chain, scoped ReadOnlySpan<int> offsets, scoped ReadOnlySpan<int> widths, |
| | | 362 | | scoped ReadOnlySpan<int> degrees, scoped ReadOnlySpan<sbyte> signs, int count, |
| | | 363 | | scoped Span<ulong> evaluation, ref int ordinal, ref int remaining, ref int lowerVariations) |
| | | 364 | | { |
| | 52506 | 365 | | ShiftFiniteRootLeft(root.LowerNumerator, 1); |
| | 52506 | 366 | | AddRoundedCylinderWord(root.LowerNumerator, 0, 1); |
| | 52506 | 367 | | root.DenominatorShift++; |
| | 52506 | 368 | | int midpointVariations = GetFiniteValueVariations(chain, offsets, widths, |
| | 52506 | 369 | | degrees, signs, count, root.LowerNumerator, root.DenominatorShift, evaluation); |
| | 52506 | 370 | | int leftCount = lowerVariations - midpointVariations; |
| | 52506 | 371 | | if (ordinal < leftCount) |
| | | 372 | | { |
| | 39775 | 373 | | if (ordinal == leftCount - 1 |
| | 39775 | 374 | | && EvaluateFiniteRootPolynomial(root.Coefficients, root.Signs, |
| | 39775 | 375 | | root.LowerNumerator, root.DenominatorShift, 0, evaluation) == 0) |
| | 865 | 376 | | root.IsRational = true; |
| | | 377 | | else |
| | 38910 | 378 | | root.LowerNumerator[0]--; // Odd midpoint, no borrow. |
| | 39775 | 379 | | remaining = leftCount; |
| | | 380 | | } |
| | | 381 | | else |
| | | 382 | | { |
| | 12731 | 383 | | ordinal -= leftCount; |
| | 12731 | 384 | | remaining -= leftCount; |
| | 12731 | 385 | | lowerVariations = midpointVariations; |
| | | 386 | | } |
| | 12731 | 387 | | } |
| | | 388 | | |
| | | 389 | | private static int GetFiniteValueCellUpperSign(FiniteAxisValueRoot root, Span<ulong> evaluation = default) |
| | | 390 | | { |
| | 17748 | 391 | | AddRoundedCylinderWord(root.LowerNumerator, 0, 1); |
| | 17748 | 392 | | int sign = evaluation.IsEmpty |
| | 17748 | 393 | | ? EvaluateFiniteRootPolynomial(root.Coefficients, root.Signs, |
| | 17748 | 394 | | root.LowerNumerator, root.DenominatorShift, 0) |
| | 17748 | 395 | | : EvaluateFiniteRootPolynomial(root.Coefficients, root.Signs, |
| | 17748 | 396 | | root.LowerNumerator, root.DenominatorShift, 0, evaluation); |
| | 35612 | 397 | | for (int index = 0; index < root.LowerNumerator.Length; index++) |
| | | 398 | | { |
| | 17806 | 399 | | ulong previous = root.LowerNumerator[index]; |
| | 17806 | 400 | | root.LowerNumerator[index] = unchecked(previous - 1); |
| | 17806 | 401 | | if (previous != 0) |
| | | 402 | | break; |
| | | 403 | | } |
| | 17748 | 404 | | return sign; |
| | | 405 | | } |
| | | 406 | | |
| | | 407 | | private static int BuildFiniteValueSturm(ReadOnlySpan<ulong> coefficients, |
| | | 408 | | ReadOnlySpan<sbyte> signs, Span<ulong> arena, Span<int> offsets, |
| | | 409 | | Span<int> widths, Span<int> degrees, Span<sbyte> sturmSigns, out int used) |
| | | 410 | | { |
| | 9895 | 411 | | sturmSigns.Clear(); |
| | 9895 | 412 | | int degree = signs.Length - 1; |
| | 9895 | 413 | | int inputWords = coefficients.Length / signs.Length; |
| | 9895 | 414 | | int words = Math.Max(1, (GetFiniteRootCoefficientBits(coefficients, signs.Length) + 63) / 64); |
| | 9895 | 415 | | offsets[0] = 0; |
| | 9895 | 416 | | widths[0] = words; |
| | 9895 | 417 | | degrees[0] = degree; |
| | 9895 | 418 | | used = (degree + 1) * words; |
| | 114548 | 419 | | for (int index = 0; index <= degree; index++) |
| | 47379 | 420 | | CopyFiniteRootMagnitude(coefficients.Slice(index * inputWords, inputWords), |
| | 47379 | 421 | | arena.Slice(index * words, words)); |
| | 9895 | 422 | | signs.CopyTo(sturmSigns); |
| | 9895 | 423 | | NormalizeFiniteValueContent(arena[..used], degree, words, sturmSigns[..9], arena[used..]); |
| | 9895 | 424 | | offsets[1] = used; |
| | 9895 | 425 | | widths[1] = words + 1; |
| | 9895 | 426 | | degrees[1] = degree - 1; |
| | 9895 | 427 | | used += degree * (words + 1); |
| | 9895 | 428 | | DifferentiateFiniteAxisPolynomial(arena[..offsets[1]], sturmSigns[..(degree + 1)], |
| | 9895 | 429 | | arena.Slice(offsets[1], degree * (words + 1)), sturmSigns.Slice(9, degree)); |
| | 9895 | 430 | | NormalizeFiniteValueContent(arena.Slice(offsets[1], degree * (words + 1)), |
| | 9895 | 431 | | degree - 1, words + 1, sturmSigns.Slice(9, 9), arena[used..]); |
| | 9895 | 432 | | int count = 2; |
| | 9895 | 433 | | if (degree == 1) |
| | 689 | 434 | | return count; |
| | | 435 | | |
| | | 436 | | // Brown's subresultant PRS uses exact scalar divisors instead of a |
| | | 437 | | // content GCD after every pseudo-remainder. Keep its raw signs until |
| | | 438 | | // construction ends, then apply the separate Sturm row orientations. |
| | | 439 | | // For C=B+64 and h(j)=2*n-1-2*j, the retained minors need <=221C |
| | | 440 | | // bits. In a pair of actual degrees d>e, the second row is the |
| | | 441 | | // subresultant of natural index d-1 even after an abnormal drop; |
| | | 442 | | // its height is <=h(d-1)C, not the looser h(e)C. Graded prem slots |
| | | 443 | | // plus two products peak at215C; exact-division scratch peaks172C. |
| | | 444 | | // Scalar powers/divisors need <=28C each. Eight such scratch slots |
| | | 445 | | // give a total <=660C, within880C plus the existing rounding padding. |
| | | 446 | | // Store raw Brown rows as S_i=2^rowPower[i]*T_i, and its scalars |
| | | 447 | | // as c=2^cPower*cOdd, beta=2^betaPower*betaOdd. Only positive |
| | | 448 | | // powers of two are removed, so row signs and Sturm variations stay |
| | | 449 | | // unchanged. Normalized retained rows/scalars only shrink; transient |
| | | 450 | | // quotients still use the existing graded pseudo-remainder storage. |
| | 9206 | 451 | | int unitBits = GetFiniteRootCoefficientBits(coefficients, signs.Length) + 64; |
| | 9206 | 452 | | int scalarWords = (28 * unitBits + 63) / 64; |
| | 9206 | 453 | | int scalarOffset = arena.Length - (8 * scalarWords + 1); |
| | 9206 | 454 | | Span<ulong> scalarWork = arena[scalarOffset..]; |
| | 9206 | 455 | | Span<ulong> c = scalarWork[..scalarWords]; |
| | 9206 | 456 | | Span<ulong> beta = scalarWork.Slice(scalarWords, scalarWords); |
| | 9206 | 457 | | CopyFiniteRootMagnitude(arena.Slice(offsets[1] + (degree - 1) * widths[1], widths[1]), c); |
| | 9206 | 458 | | int cPower = CountRoundedCylinderTrailingZeroes(c); |
| | 9206 | 459 | | ShiftRoundedCylinderWideRight(c, cPower); |
| | 9206 | 460 | | int cSign = -sturmSigns[9 + degree - 1]; |
| | 9206 | 461 | | beta.Clear(); |
| | 9206 | 462 | | beta[0] = 1; // Initial degree gap is one: beta=(-1)^(gap+1). |
| | 9206 | 463 | | int betaSign = 1; |
| | 9206 | 464 | | int betaPower = 0; |
| | 9206 | 465 | | Span<int> rowPowers = stackalloc int[9]; |
| | 9206 | 466 | | rowPowers[0] = rowPowers[1] = 0; |
| | 9206 | 467 | | Span<sbyte> rowSigns = stackalloc sbyte[9]; |
| | 9206 | 468 | | rowSigns[0] = rowSigns[1] = 1; |
| | 31244 | 469 | | while (degrees[count - 1] > 0) |
| | | 470 | | { |
| | 24905 | 471 | | int left = count - 2; |
| | 24905 | 472 | | int right = count - 1; |
| | 24905 | 473 | | int power = degrees[left] - degrees[right] + 1; |
| | 24905 | 474 | | int pseudoSign = (power & 1) == 0 ? 1 : sturmSigns[right * 9 + degrees[right]]; |
| | 24905 | 475 | | int nextDegree = BuildFiniteValueSubresultant( |
| | 24905 | 476 | | arena.Slice(offsets[left], (degrees[left] + 1) * widths[left]), |
| | 24905 | 477 | | sturmSigns.Slice(left * 9, degrees[left] + 1), |
| | 24905 | 478 | | arena.Slice(offsets[right], (degrees[right] + 1) * widths[right]), |
| | 24905 | 479 | | sturmSigns.Slice(right * 9, degrees[right] + 1), |
| | 24905 | 480 | | beta, betaSign * pseudoSign, arena.Slice(used, scalarOffset - used), |
| | 24905 | 481 | | sturmSigns.Slice(count * 9, 9), out int nextWords, out int removedPower); |
| | 24905 | 482 | | if (nextDegree < 0) |
| | | 483 | | break; |
| | | 484 | | // prem is homogeneous: the omitted input scale is |
| | | 485 | | // rowPower[left]+power*rowPower[right]. After division by |
| | | 486 | | // betaOdd, stripping 2^removedPower leaves the raw row exponent |
| | | 487 | | // below. The difference before adding removedPower may be |
| | | 488 | | // negative; integrality guarantees only the final sum is >=0. |
| | 22038 | 489 | | rowPowers[count] = rowPowers[left] + power * rowPowers[right] |
| | 22038 | 490 | | - betaPower + removedPower; |
| | | 491 | | System.Diagnostics.Debug.Assert(rowPowers[count] >= 0); |
| | 22038 | 492 | | offsets[count] = used; |
| | 22038 | 493 | | widths[count] = nextWords; |
| | 22038 | 494 | | degrees[count] = nextDegree; |
| | | 495 | | // prem=lc(right)^(gap+1)*rem. Multiplying each completed raw |
| | | 496 | | // row by this sign makes it a positive multiple of -rem(Sturm). |
| | 22038 | 497 | | rowSigns[count] = (sbyte)(-rowSigns[left] * betaSign * pseudoSign); |
| | 22038 | 498 | | used += (nextDegree + 1) * nextWords; |
| | 22038 | 499 | | if (nextDegree > 0) |
| | 15699 | 500 | | UpdateFiniteValueSubresultantScalars( |
| | 15699 | 501 | | arena.Slice(offsets[right] + degrees[right] * widths[right], widths[right]), |
| | 15699 | 502 | | sturmSigns[right * 9 + degrees[right]], rowPowers[right], |
| | 15699 | 503 | | arena.Slice(offsets[count] + nextDegree * nextWords, nextWords), |
| | 15699 | 504 | | sturmSigns[count * 9 + nextDegree], rowPowers[count], degrees[right] - nextDegree, |
| | 15699 | 505 | | scalarWork, scalarWords, ref cSign, ref cPower, out betaSign, out betaPower); |
| | 22038 | 506 | | count++; |
| | | 507 | | } |
| | 62488 | 508 | | for (int row = 2; row < count; row++) |
| | 176514 | 509 | | for (int index = 0; index <= degrees[row]; index++) |
| | 66219 | 510 | | sturmSigns[row * 9 + index] *= rowSigns[row]; |
| | 9206 | 511 | | return count; |
| | | 512 | | } |
| | | 513 | | |
| | | 514 | | private static void UpdateFiniteValueSubresultantScalars(ReadOnlySpan<ulong> previousLeading, |
| | | 515 | | int previousLeadingSign, int previousRowPower, ReadOnlySpan<ulong> leading, |
| | | 516 | | int leadingSign, int rowPower, int gap, Span<ulong> work, int words, |
| | | 517 | | ref int cSign, ref int cPower, out int betaSign, out int betaPower) |
| | | 518 | | { |
| | 15699 | 519 | | Span<ulong> c = work[..words]; |
| | 15699 | 520 | | Span<ulong> beta = work.Slice(words, words); |
| | 15699 | 521 | | Span<ulong> numerator = work.Slice(2 * words, words); |
| | 15699 | 522 | | Span<ulong> denominator = work.Slice(3 * words, words); |
| | 15699 | 523 | | Span<ulong> product = work.Slice(4 * words, words); |
| | | 524 | | // beta=-previousLeading*c^gap; c_next=(-leading)^gap/c^(gap-1). |
| | | 525 | | // Each leading coefficient has both the row's omitted power and |
| | | 526 | | // its stored leading term's power. Track those separately and use |
| | | 527 | | // only odd magnitudes in the unchanged scalar recurrence. The odd |
| | | 528 | | // denominator divides the odd numerator exactly, since it is coprime |
| | | 529 | | // to every omitted power of two, including abnormal degree drops. |
| | 15699 | 530 | | int previousLeadingPower = CountRoundedCylinderTrailingZeroes(previousLeading); |
| | 15699 | 531 | | int leadingPower = CountRoundedCylinderTrailingZeroes(leading); |
| | 15699 | 532 | | betaPower = previousRowPower + previousLeadingPower + gap * cPower; |
| | 15699 | 533 | | int nextCPower = gap * (rowPower + leadingPower) - (gap - 1) * cPower; |
| | | 534 | | System.Diagnostics.Debug.Assert(betaPower >= 0 && nextCPower >= 0); |
| | | 535 | | // Reuse the denominator slot for leading odd parts before it is |
| | | 536 | | // needed for cOdd^(gap-1); no extra coefficient buffer is retained. |
| | 15699 | 537 | | CopyFiniteRootMagnitude(previousLeading, denominator); |
| | 15699 | 538 | | ShiftRoundedCylinderWideRight(denominator, previousLeadingPower); |
| | 15699 | 539 | | PowerFiniteValueMagnitude(c, gap, numerator, product); |
| | 15699 | 540 | | WideArithmetic.MultiplyMagnitudes(denominator, numerator, beta); |
| | 15699 | 541 | | betaSign = -previousLeadingSign * ((gap & 1) == 0 ? 1 : cSign); |
| | 15699 | 542 | | if (gap == 1) |
| | | 543 | | { |
| | 15671 | 544 | | CopyFiniteRootMagnitude(leading, c); |
| | 15671 | 545 | | ShiftRoundedCylinderWideRight(c, leadingPower); |
| | 15671 | 546 | | cSign = -leadingSign; |
| | 15671 | 547 | | cPower = nextCPower; |
| | 15671 | 548 | | return; |
| | | 549 | | } |
| | 28 | 550 | | CopyFiniteRootMagnitude(leading, denominator); |
| | 28 | 551 | | ShiftRoundedCylinderWideRight(denominator, leadingPower); |
| | 28 | 552 | | PowerFiniteValueMagnitude(denominator, gap, numerator, product); |
| | 28 | 553 | | PowerFiniteValueMagnitude(c, gap - 1, denominator, product); |
| | 28 | 554 | | int nextSign = ((gap & 1) == 0 ? 1 : -leadingSign) |
| | 28 | 555 | | * ((gap & 1) == 0 ? cSign : 1); |
| | 28 | 556 | | Span<ulong> remainder = work.Slice(5 * words, words); |
| | 28 | 557 | | WideArithmetic.DivideMagnitudes(numerator, denominator, c, remainder, work[(6 * words)..]); |
| | | 558 | | System.Diagnostics.Debug.Assert(GetRoundedCylinderWideLength(remainder) == 0); |
| | 28 | 559 | | cSign = nextSign; |
| | 28 | 560 | | cPower = nextCPower; |
| | 28 | 561 | | } |
| | | 562 | | |
| | | 563 | | private static void PowerFiniteValueMagnitude(ReadOnlySpan<ulong> value, int power, |
| | | 564 | | Span<ulong> result, Span<ulong> product) |
| | | 565 | | { |
| | 15755 | 566 | | result.Clear(); |
| | 15755 | 567 | | result[0] = 1; |
| | 63276 | 568 | | for (int index = 0; index < power; index++) |
| | | 569 | | { |
| | 15883 | 570 | | WideArithmetic.MultiplyMagnitudes(result, value, product); |
| | 15883 | 571 | | product.CopyTo(result); |
| | | 572 | | } |
| | 15755 | 573 | | } |
| | | 574 | | |
| | | 575 | | private static int BuildFiniteValueSubresultant(ReadOnlySpan<ulong> left, |
| | | 576 | | ReadOnlySpan<sbyte> leftSigns, ReadOnlySpan<ulong> right, |
| | | 577 | | ReadOnlySpan<sbyte> rightSigns, ReadOnlySpan<ulong> beta, int resultSign, |
| | | 578 | | Span<ulong> work, Span<sbyte> resultSigns, |
| | | 579 | | out int resultWords, out int removedPower) |
| | | 580 | | { |
| | 24905 | 581 | | int d = leftSigns.Length - 1; |
| | 24905 | 582 | | int e = rightSigns.Length - 1; |
| | 24905 | 583 | | int leftWords = left.Length / (d + 1); |
| | 24905 | 584 | | int rightWords = right.Length / (e + 1); |
| | 24905 | 585 | | int leftBits = GetFiniteRootCoefficientBits(left, d + 1); |
| | 24905 | 586 | | int rightBits = GetFiniteRootCoefficientBits(right, e + 1); |
| | 24905 | 587 | | Span<int> offsets = stackalloc int[9]; |
| | 24905 | 588 | | Span<int> widths = stackalloc int[9]; |
| | 24905 | 589 | | Span<sbyte> signs = stackalloc sbyte[9]; |
| | 24905 | 590 | | int used = 0; |
| | 24905 | 591 | | int maximumWords = 0; |
| | 293074 | 592 | | for (int index = 0; index <= d; index++) |
| | | 593 | | { |
| | 121632 | 594 | | int cancellations = Math.Min(d - e + 1, d - index); |
| | 121632 | 595 | | int width = (leftBits + cancellations * (rightBits + 1) + 127) / 64; |
| | 121632 | 596 | | offsets[index] = used; |
| | 121632 | 597 | | widths[index] = width; |
| | 121632 | 598 | | used += width; |
| | 121632 | 599 | | maximumWords = Math.Max(maximumWords, width); |
| | 121632 | 600 | | CopyFiniteRootMagnitude(left.Slice(index * leftWords, leftWords), |
| | 121632 | 601 | | work.Slice(offsets[index], width)); |
| | 121632 | 602 | | signs[index] = leftSigns[index]; |
| | | 603 | | } |
| | 24905 | 604 | | Span<ulong> product = work.Slice(used, maximumWords); |
| | 24905 | 605 | | Span<ulong> subtract = work.Slice(used + maximumWords, maximumWords); |
| | 24905 | 606 | | ReadOnlySpan<ulong> divisorLeading = right.Slice(e * rightWords, rightWords); |
| | 149534 | 607 | | for (int degree = d; degree >= e; degree--) |
| | | 608 | | { |
| | | 609 | | // Even a zero leading term must multiply every remaining slot: |
| | | 610 | | // Brown's prem uses exactly d-e+1 leading-coefficient factors. |
| | 49862 | 611 | | ReadOnlySpan<ulong> leading = work.Slice(offsets[degree], widths[degree]); |
| | 49862 | 612 | | int shift = degree - e; |
| | 437294 | 613 | | for (int index = 0; index < degree; index++) |
| | | 614 | | { |
| | 168785 | 615 | | Span<ulong> coefficient = work.Slice(offsets[index], widths[index]); |
| | 168785 | 616 | | WideArithmetic.MultiplyMagnitudes(coefficient, divisorLeading, product); |
| | 168785 | 617 | | int sign = signs[index]; |
| | 168785 | 618 | | if (signs[degree] != 0 && index >= shift && rightSigns[index - shift] != 0) |
| | | 619 | | { |
| | 139739 | 620 | | WideArithmetic.MultiplyMagnitudes(leading, |
| | 139739 | 621 | | right.Slice((index - shift) * rightWords, rightWords), subtract); |
| | 139739 | 622 | | WideArithmetic.AddShiftedSignedMagnitude(subtract, |
| | 139739 | 623 | | -signs[degree] * rightSigns[e] * rightSigns[index - shift], 0, product, ref sign); |
| | | 624 | | } |
| | 168785 | 625 | | product[..widths[index]].CopyTo(coefficient); |
| | 168785 | 626 | | signs[index] = (sbyte)sign; |
| | | 627 | | } |
| | 49862 | 628 | | signs[degree] = 0; |
| | | 629 | | } |
| | 24905 | 630 | | int remainderDegree = e - 1; |
| | 30456 | 631 | | while (remainderDegree >= 0 && signs[remainderDegree] == 0) |
| | 5551 | 632 | | remainderDegree--; |
| | 24905 | 633 | | resultWords = 0; |
| | 24905 | 634 | | removedPower = 0; |
| | 24905 | 635 | | if (remainderDegree < 0) |
| | 2867 | 636 | | return -1; |
| | | 637 | | |
| | | 638 | | // All surviving slots have the same width. Higher discarded slots |
| | | 639 | | // become exact-division scratch; no pseudo-division products remain live. |
| | 22038 | 640 | | int remainderWords = widths[0]; |
| | 22038 | 641 | | int remainderLength = (remainderDegree + 1) * remainderWords; |
| | 22038 | 642 | | Span<ulong> remainder = work[..remainderLength]; |
| | | 643 | | // beta stores only its odd part. It divides this smaller prem |
| | | 644 | | // exactly: divisibility of the raw Brown numerator cannot depend on |
| | | 645 | | // its omitted power of two when the divisor is odd. |
| | 22038 | 646 | | int betaWords = GetRoundedCylinderWideLength(beta); |
| | 22038 | 647 | | if (betaWords != 1 || beta[0] != 1) |
| | | 648 | | { |
| | 10730 | 649 | | Span<ulong> integerRemainder = work.Slice(remainderLength, remainderWords); |
| | 10730 | 650 | | Span<ulong> division = work.Slice(remainderLength + remainderWords, 2 * remainderWords + 1); |
| | 77482 | 651 | | for (int index = 0; index <= remainderDegree; index++) |
| | | 652 | | { |
| | 28011 | 653 | | Span<ulong> coefficient = remainder.Slice(index * remainderWords, remainderWords); |
| | 28011 | 654 | | WideArithmetic.DivideMagnitudes(coefficient, beta[..betaWords], coefficient, integerRemainder, division) |
| | | 655 | | System.Diagnostics.Debug.Assert(GetRoundedCylinderWideLength(integerRemainder) == 0); |
| | | 656 | | } |
| | | 657 | | } |
| | | 658 | | // U=prem(T_left,T_right)/betaOdd can exceed the raw next row when |
| | | 659 | | // its omitted scale exponent is negative. It fits the prem scratch, |
| | | 660 | | // but must lose its common power of two BEFORE retained width/used |
| | | 661 | | // accounting; only the normalized T_next is bounded by the raw row. |
| | 22038 | 662 | | removedPower = NormalizeFiniteAxisPolynomialPowerOfTwo( |
| | 22038 | 663 | | remainder, signs[..(remainderDegree + 1)]); |
| | 176514 | 664 | | for (int index = 0; index <= remainderDegree; index++) |
| | 66219 | 665 | | resultWords = Math.Max(resultWords, |
| | 66219 | 666 | | GetRoundedCylinderWideLength(remainder.Slice(index * remainderWords, remainderWords))); |
| | 176514 | 667 | | for (int index = 0; index <= remainderDegree; index++) |
| | | 668 | | { |
| | 66219 | 669 | | remainder.Slice(index * remainderWords, resultWords) |
| | 66219 | 670 | | .CopyTo(work.Slice(index * resultWords, resultWords)); |
| | 66219 | 671 | | resultSigns[index] = (sbyte)(resultSign * signs[index]); |
| | | 672 | | } |
| | 22038 | 673 | | return remainderDegree; |
| | | 674 | | } |
| | | 675 | | |
| | | 676 | | private static void NormalizeFiniteValueContent(Span<ulong> coefficients, int degree, |
| | | 677 | | int words, ReadOnlySpan<sbyte> signs, Span<ulong> scratch) |
| | | 678 | | { |
| | 20812 | 679 | | if (degree == 0) |
| | | 680 | | { |
| | 1080 | 681 | | coefficients.Clear(); |
| | 1080 | 682 | | coefficients[0] = 1; |
| | 1080 | 683 | | return; |
| | | 684 | | } |
| | | 685 | | // Variable scaling introduces large shared powers of two in later |
| | | 686 | | // remainders too. Strip those exactly before Euclid works on the |
| | | 687 | | // remaining integer content; roots and signed variations are unchanged. |
| | 19732 | 688 | | NormalizeFiniteAxisPolynomialPowerOfTwo(coefficients, signs[..(degree + 1)]); |
| | 19732 | 689 | | Span<ulong> content = scratch[..words]; |
| | 19732 | 690 | | Span<ulong> other = scratch.Slice(words, words); |
| | 19732 | 691 | | Span<ulong> division = scratch.Slice(2 * words, 2 * words + 1); |
| | 19732 | 692 | | content.Clear(); |
| | 146804 | 693 | | for (int index = 0; index <= degree; index++) |
| | | 694 | | { |
| | 68541 | 695 | | if (signs[index] == 0) |
| | | 696 | | continue; |
| | 66524 | 697 | | coefficients.Slice(index * words, words).CopyTo(other); |
| | 66524 | 698 | | WideArithmetic.GetMagnitudeGreatestCommonDivisor(content, other, division); |
| | 66524 | 699 | | if (content[0] == 1 && GetRoundedCylinderWideLength(content) == 1) |
| | 14871 | 700 | | return; |
| | | 701 | | } |
| | 64780 | 702 | | for (int index = 0; index <= degree; index++) |
| | | 703 | | { |
| | 27529 | 704 | | Span<ulong> coefficient = coefficients.Slice(index * words, words); |
| | 27529 | 705 | | WideArithmetic.DivideMagnitudes(coefficient, content, coefficient, other, division); |
| | | 706 | | System.Diagnostics.Debug.Assert(GetRoundedCylinderWideLength(other) == 0); |
| | | 707 | | } |
| | 4861 | 708 | | } |
| | | 709 | | |
| | | 710 | | private static int GetFiniteValueVariations(ReadOnlySpan<ulong> chain, |
| | | 711 | | ReadOnlySpan<int> offsets, ReadOnlySpan<int> widths, ReadOnlySpan<int> degrees, |
| | | 712 | | ReadOnlySpan<sbyte> signs, int count, ReadOnlySpan<ulong> numerator, int shift, |
| | | 713 | | Span<ulong> evaluation) |
| | | 714 | | { |
| | 113637 | 715 | | int previous = 0; |
| | 113637 | 716 | | int variations = 0; |
| | 113637 | 717 | | int numeratorBits = GetFiniteRootBits(numerator); |
| | 1347948 | 718 | | for (int index = 0; index < count; index++) |
| | | 719 | | { |
| | 560337 | 720 | | ReadOnlySpan<ulong> row = chain.Slice(offsets[index], (degrees[index] + 1) * widths[index]); |
| | 560337 | 721 | | ReadOnlySpan<sbyte> rowSigns = signs.Slice(index * 9, degrees[index] + 1); |
| | 560337 | 722 | | int sign = 0; |
| | 560337 | 723 | | if (degrees[index] > 0 && numeratorBits > 0 && numeratorBits <= shift) |
| | | 724 | | { |
| | | 725 | | // A certified point sign can replace homogeneous evaluation. |
| | | 726 | | // Reuse the caller's returning evaluation scratch; cap only |
| | | 727 | | // work precision to fit, never the sign acceptance threshold. |
| | 432295 | 728 | | int countBits = GetFiniteValueCeilingLog2(degrees[index] + 1); |
| | 432295 | 729 | | int productWords = (shift + 64) / 64; |
| | 432295 | 730 | | int maximumPrecision = 64 * ((evaluation.Length - productWords) / 2) - countBits - 64; |
| | 432295 | 731 | | int precision = Math.Min(maximumPrecision, |
| | 432295 | 732 | | shift + FiniteValuePointGuardBits + degrees[index] * (shift - numeratorBits + 1)); |
| | | 733 | | System.Diagnostics.Debug.Assert(precision > 0); |
| | 432295 | 734 | | int words = (precision + countBits + 127) / 64; |
| | 432295 | 735 | | sign = GetFiniteValueApproximateSignCore(numerator, shift, row, rowSigns, |
| | 432295 | 736 | | precision, GetFiniteRootCoefficientBits(row, rowSigns.Length), |
| | 432295 | 737 | | evaluation[..words], evaluation.Slice(words, words + productWords), out _); |
| | | 738 | | } |
| | | 739 | | // Uncertainty is not zero evidence. At a vanishing Sturm row, |
| | | 740 | | // its first nonzero derivative supplies the exact right-limit sign. |
| | 1462876 | 741 | | for (int derivative = 0; derivative <= degrees[index] && sign == 0; derivative++) |
| | 171101 | 742 | | sign = EvaluateFiniteRootPolynomial(row, rowSigns, numerator, shift, derivative, evaluation); |
| | 560337 | 743 | | if (previous != 0 && previous != sign) |
| | 179135 | 744 | | variations++; |
| | 560337 | 745 | | previous = sign; |
| | | 746 | | } |
| | 113637 | 747 | | return variations; |
| | | 748 | | } |
| | | 749 | | } |