< Summary

Information
Class: FixedMathSharp.Geometry.CylinderPairValuePolynomial
Assembly: FixedMathSharp
File(s): /home/runner/work/FixedMathSharp/FixedMathSharp/src/FixedMathSharp/Geometry/Wide/Convex/CylinderPairValuePolynomial.cs
Line coverage
100%
Covered lines: 186
Uncovered lines: 0
Coverable lines: 186
Total lines: 312
Line coverage: 100%
Branch coverage
100%
Covered branches: 30
Total branches: 30
Branch coverage: 100%
Method coverage

Feature is only available for sponsors

Upgrade to PRO version

Metrics

MethodBranch coverage Crap Score Cyclomatic complexity Line coverage
Build(...)100%2222100%
.ctor(...)100%11100%
Clear()100%11100%
.ctor(...)100%11100%
Slice(...)100%44100%
Multiply(...)100%22100%
Add(...)100%22100%

File(s)

/home/runner/work/FixedMathSharp/FixedMathSharp/src/FixedMathSharp/Geometry/Wide/Convex/CylinderPairValuePolynomial.cs

#LineLine coverage
 1//=======================================================================
 2// CylinderPairValuePolynomial.cs
 3//=======================================================================
 4// MIT License, Copyright (c) 2024-present David Oravsky (mrdav30)
 5// See LICENSE file in the project root for full license information.
 6//=======================================================================
 7
 8using System;
 9
 10namespace FixedMathSharp.Geometry;
 11
 12/// <summary>
 13/// Constructs the exact, unnormalized squared-distance discriminant for two rims.
 14/// </summary>
 15internal static class CylinderPairValuePolynomial
 16{
 17    internal const int InvariantCount = 11;
 18    internal const int CoefficientCount = 9;
 19    internal const int ScratchCoefficientCount = 73;
 20
 21    /// <summary>
 22    /// Invariant slots are u², w², b², ra², rb², c², c.u, c.w, b.u,
 23    /// b.w, b.c. Every scalar has <paramref name="words"/> little-endian
 24    /// magnitude words and a separate canonical sign. Returns the degree,
 25    /// or -1 for zero, of three times the cubic-pencil discriminant.
 26    /// </summary>
 27    /// <remarks>
 28    /// The caller supplies disjoint input, output, and scratch storage: eleven,
 29    /// nine, and seventy-three coefficient slots, respectively. The basis u,w
 30    /// is orthogonal in the first circle's plane; independently reduced basis
 31    /// vectors are allowed. Geometry admission and physical S scaling belong
 32    /// to the caller. No coefficient content or repeated factor is removed.
 33    ///
 34    /// Every intermediate must fit words*64 bits. A conservative sufficient
 35    /// bound for arbitrary input magnitudes below 2^b, b at least one, is
 36    /// words*64 at least 32*b+128. A caller may instead use a tighter proven
 37    /// geometry-specific bound, including convolution sum carries. Magnitude
 38    /// operations do not detect truncation. There is no nested stack scratch;
 39    /// construction storage need not remain live during root isolation.
 40    /// An optional direction adds seventy-two scratch slots and includes its
 41    /// input magnitudes in b. The raw scale is W=3*discriminant for both value
 42    /// and derivative: allow up to two additional coefficient bits relative
 43    /// to a bound for the unscaled discriminant. Default direction storage
 44    /// performs no derivative convolutions.
 45    /// Scratch slots use min(words, ceil((32*b+128)/64)) words; input and
 46    /// retained output slots keep the caller's original words stride.
 47    /// </remarks>
 48    internal static int Build(ReadOnlySpan<ulong> invariants,
 49        ReadOnlySpan<sbyte> invariantSigns, int words,
 50        Span<ulong> coefficients, Span<sbyte> signs,
 51        Span<ulong> scratch, Span<sbyte> scratchSigns,
 52        CylinderPairValueDerivative derivative = default)
 53    {
 640454        bool differentiate = !derivative.Invariants.IsEmpty;
 640455        int inputBits = 1;
 15369656        for (int index = 0; index < InvariantCount; index++)
 7044457            inputBits = Math.Max(inputBits,
 7044458                WideArithmetic.GetMagnitudeBitLength(invariants.Slice(index * words, words)));
 640459        if (differentiate)
 7179260            for (int index = 0; index < CylinderPairValueDerivative.InvariantCount; index++)
 3076861                inputBits = Math.Max(inputBits,
 3076862                    WideArithmetic.GetMagnitudeBitLength(derivative.Invariants.Slice(index * words, words)));
 63
 64        // Both the generic bound and the caller's geometry-specific capacity
 65        // bound every intermediate. Their minimum avoids scanning padded words
 66        // without changing retained coefficients or truncating dual products.
 640467        int workWords = Math.Min(words, (32 * inputBits + 128 + 63) / 64);
 640468        scratch = scratch.Slice(0, ScratchCoefficientCount * workWords);
 640469        scratchSigns = scratchSigns.Slice(0, ScratchCoefficientCount);
 640470        scratch.Clear();
 640471        scratchSigns.Clear();
 15369672        for (int index = 0; index < InvariantCount; index++)
 7044473            invariants.Slice(index * words, workWords)
 7044474                .CopyTo(scratch.Slice(index * workWords, workWords));
 640475        invariantSigns.Slice(0, InvariantCount).CopyTo(scratchSigns);
 640476        if (differentiate)
 77        {
 512878            derivative.Scratch.Slice(0, CylinderPairValueDerivative.ScratchCoefficientCount * workWords).Clear();
 512879            derivative.ScratchSigns.Slice(0, CylinderPairValueDerivative.ScratchCoefficientCount).Clear();
 7179280            for (int index = 0; index < CylinderPairValueDerivative.InvariantCount; index++)
 81            {
 3076882                int target = index < 5 ? index + 3 : 10;
 3076883                derivative.Invariants.Slice(index * words, workWords)
 3076884                    .CopyTo(derivative.Scratch.Slice(target * workWords, workWords));
 3076885                derivative.ScratchSigns[target] = derivative.InvariantSigns[index];
 86            }
 87        }
 88
 640489        ConstructionScratch workspace = new(scratch, scratchSigns, workWords, derivative);
 640490        Polynomial g = workspace.Slice(0, 1);
 640491        Polynomial h = workspace.Slice(1, 1);
 640492        Polynomial delta = workspace.Slice(2, 1);
 640493        Polynomial rho = workspace.Slice(3, 1);
 640494        Polynomial tau = workspace.Slice(4, 1);
 640495        Polynomial cSquared = workspace.Slice(5, 1);
 640496        Polynomial c0 = workspace.Slice(6, 1);
 640497        Polynomial c1 = workspace.Slice(7, 1);
 640498        Polynomial b0 = workspace.Slice(8, 1);
 640499        Polynomial b1 = workspace.Slice(9, 1);
 6404100        Polynomial j = workspace.Slice(10, 1);
 101
 6404102        Polynomial e = workspace.Slice(11, 1);
 6404103        Polynomial bigK = workspace.Slice(12, 1);
 6404104        Polynomial s0 = workspace.Slice(13, 1);
 6404105        Polynomial s1 = workspace.Slice(14, 1);
 6404106        Polynomial s2 = workspace.Slice(15, 1);
 6404107        Polynomial s3 = workspace.Slice(16, 1);
 6404108        Polynomial s4 = workspace.Slice(17, 1);
 6404109        Polynomial k = workspace.Slice(18, 2);
 6404110        Polynomial t0 = workspace.Slice(20, 2);
 6404111        Polynomial t1 = workspace.Slice(22, 2);
 6404112        Polynomial bigL = workspace.Slice(24, 3);
 6404113        Polynomial d = workspace.Slice(27, 3);
 6404114        Polynomial a = workspace.Slice(30, 1);
 6404115        Polynomial b = workspace.Slice(31, 3);
 6404116        Polynomial c = workspace.Slice(34, 3);
 6404117        Polynomial constant = workspace.Slice(37, 2);
 6404118        Polynomial p = workspace.Slice(39, 5);
 6404119        Polynomial r = workspace.Slice(44, 5);
 6404120        Polynomial n = workspace.Slice(49, 5);
 6404121        Polynomial temporary0 = workspace.Slice(54, 9);
 6404122        Polynomial temporary1 = workspace.Slice(63, 9);
 6404123        Span<ulong> product = scratch.Slice(72 * workWords, workWords);
 124
 125        // e=(c0*b1-c1*b0)^2.
 6404126        Multiply(s0, c0, b1, product);
 6404127        Multiply(s1, c1, b0, product);
 6404128        Add(s0, s1, -1);
 6404129        Multiply(e, s0, s0, product);
 130
 131        // K=delta*(g*c1^2+h*c0^2)+tau*(g*b1^2+h*b0^2).
 6404132        Multiply(s0, c1, c1, product);
 6404133        Multiply(s1, g, s0, product);
 6404134        Multiply(s0, c0, c0, product);
 6404135        Multiply(s2, h, s0, product);
 6404136        Add(s1, s2, 1);
 6404137        Multiply(s3, delta, s1, product);
 6404138        Multiply(s0, b1, b1, product);
 6404139        Multiply(s1, g, s0, product);
 6404140        Multiply(s0, b0, b0, product);
 6404141        Multiply(s2, h, s0, product);
 6404142        Add(s1, s2, 1);
 6404143        Multiply(s4, tau, s1, product);
 6404144        Add(bigK, s3, 1);
 6404145        Add(bigK, s4, 1);
 146
 147        // k=c^2+rho-tau-S.
 6404148        Add(k, cSquared, 1);
 6404149        Add(k, rho, 1);
 6404150        Add(k, tau, -1);
 6404151        k.Magnitudes[workWords] = 1UL;
 6404152        k.Signs[1] = -1;
 153
 154        // L=g*(2*j*c1-k*b1)^2+h*(2*j*c0-k*b0)^2.
 6404155        Multiply(s0, j, c1, product);
 6404156        Add(t0, s0, 2);
 6404157        Multiply(t1, k, b1, product);
 6404158        Add(t0, t1, -1);
 6404159        Multiply(temporary0, t0, t0, product);
 6404160        Multiply(bigL, g, temporary0, product);
 6404161        t0.Clear();
 6404162        Multiply(s0, j, c0, product);
 6404163        Add(t0, s0, 2);
 6404164        Multiply(t1, k, b0, product);
 6404165        Add(t0, t1, -1);
 6404166        Multiply(temporary0, t0, t0, product);
 6404167        Multiply(temporary1, h, temporary0, product);
 6404168        Add(bigL, temporary1, 1);
 169
 170        // Tzz=d=delta*(k^2-4*tau*S)+4*tau*j^2.
 6404171        Multiply(temporary0, k, k, product);
 6404172        Multiply(d, delta, temporary0, product);
 6404173        Multiply(s0, delta, tau, product);
 6404174        Add(d, s0, -4, 1);
 6404175        Multiply(s0, j, j, product);
 6404176        Multiply(s1, tau, s0, product);
 6404177        Add(d, s1, 4);
 178
 179        // det(T-lambda*C)=A*lambda^3+B*lambda^2+C1*lambda+D.
 6404180        Multiply(s0, g, h, product);
 6404181        Multiply(a, s0, rho, product);
 6404182        Multiply(b, s0, d, product);
 6404183        Multiply(s1, rho, bigK, product);
 6404184        Add(b, s1, -4);
 185
 6404186        temporary0.Clear();
 6404187        Multiply(s0, rho, e, product);
 6404188        Add(temporary0, s0, 4);
 6404189        Add(temporary0, bigK, 4, 1);
 6404190        Add(temporary0, bigL, -1);
 6404191        Multiply(s0, delta, tau, product);
 6404192        Multiply(temporary1, s0, temporary0, product);
 6404193        Add(c, temporary1, 4);
 194
 195        // The completed-square form makes D linear, not quadratic.
 6404196        Multiply(s0, delta, delta, product);
 6404197        Multiply(s1, tau, tau, product);
 6404198        Multiply(s2, s0, s1, product);
 6404199        Multiply(s3, s2, e, product);
 6404200        Add(constant, s3, -64, 1);
 201
 202        // P=B^2-3*A*C1; R=C1^2-3*B*D; N=9*A*D-B*C1.
 6404203        Multiply(p, b, b, product);
 6404204        Multiply(temporary0, a, c, product);
 6404205        Add(p, temporary0, -3);
 6404206        Multiply(r, c, c, product);
 6404207        Multiply(temporary0, b, constant, product);
 6404208        Add(r, temporary0, -3);
 6404209        Multiply(temporary0, a, constant, product);
 6404210        Add(n, temporary0, 9);
 6404211        Multiply(temporary0, b, c, product);
 6404212        Add(n, temporary0, -1);
 213
 214        // Three times the discriminant; preserve this exact envelope scale.
 6404215        Multiply(temporary0, p, r, product);
 6404216        Multiply(temporary1, n, n, product);
 6404217        Polynomial output = new(coefficients.Slice(0, CoefficientCount * words),
 6404218            signs.Slice(0, CoefficientCount),
 6404219            differentiate ? derivative.Coefficients.Slice(0, CoefficientCount * words) : default,
 6404220            differentiate ? derivative.Signs.Slice(0, CoefficientCount) : default);
 6404221        output.Clear();
 6404222        Add(output, temporary0, 4);
 6404223        Add(output, temporary1, -1);
 6404224        int degree = CoefficientCount - 1;
 6424225        while (degree >= 0 && output.Signs[degree] == 0)
 20226            degree--;
 6404227        return degree;
 228    }
 229
 230    // These are construction views, not retained roots or a second arithmetic
 231    // implementation. Scalar products and signed sums use WideArithmetic;
 232    // every destination in the schedule is disjoint from its sources.
 233    private readonly ref struct Polynomial
 234    {
 235        internal readonly Span<ulong> Magnitudes;
 236        internal readonly Span<sbyte> Signs;
 237        internal readonly Span<ulong> DerivativeMagnitudes;
 238        internal readonly Span<sbyte> DerivativeSigns;
 239
 240        internal Polynomial(Span<ulong> magnitudes, Span<sbyte> signs,
 241            Span<ulong> derivativeMagnitudes = default, Span<sbyte> derivativeSigns = default)
 242        {
 211332243            Magnitudes = magnitudes;
 211332244            Signs = signs;
 211332245            DerivativeMagnitudes = derivativeMagnitudes;
 211332246            DerivativeSigns = derivativeSigns;
 211332247        }
 248
 249        internal void Clear()
 250        {
 19212251            Magnitudes.Clear();
 19212252            Signs.Clear();
 19212253            DerivativeMagnitudes.Clear();
 19212254            DerivativeSigns.Clear();
 19212255        }
 256    }
 257
 258    private readonly ref struct ConstructionScratch
 259    {
 260        private readonly Span<ulong> _magnitudes;
 261        private readonly Span<sbyte> _signs;
 262        private readonly int _words;
 263        private readonly CylinderPairValueDerivative _derivative;
 264
 265        internal ConstructionScratch(Span<ulong> magnitudes, Span<sbyte> signs,
 266            int words, CylinderPairValueDerivative derivative)
 267        {
 6404268            _magnitudes = magnitudes;
 6404269            _signs = signs;
 6404270            _words = words;
 6404271            _derivative = derivative;
 6404272        }
 273
 274        internal Polynomial Slice(int first, int count) =>
 204928275            new(_magnitudes.Slice(first * _words, count * _words),
 204928276                _signs.Slice(first, count),
 204928277                _derivative.Invariants.IsEmpty ? default :
 204928278                    _derivative.Scratch.Slice(first * _words, count * _words),
 204928279                _derivative.Invariants.IsEmpty ? default :
 204928280                    _derivative.ScratchSigns.Slice(first, count));
 281    }
 282
 283    private static void Multiply(Polynomial destination, Polynomial left,
 284        Polynomial right, Span<ulong> product)
 285    {
 288180286        WideFiniteAxisIntersection.MultiplyFiniteAxisPolynomials(
 288180287            left.Magnitudes, left.Signs, right.Magnitudes, right.Signs,
 288180288            destination.Magnitudes, destination.Signs, product);
 288180289        if (!destination.DerivativeSigns.IsEmpty)
 290        {
 230760291            WideFiniteAxisIntersection.MultiplyFiniteAxisPolynomials(
 230760292                left.DerivativeMagnitudes, left.DerivativeSigns, right.Magnitudes, right.Signs,
 230760293                destination.DerivativeMagnitudes, destination.DerivativeSigns, product);
 230760294            WideFiniteAxisIntersection.MultiplyFiniteAxisPolynomials(
 230760295                left.Magnitudes, left.Signs, right.DerivativeMagnitudes, right.DerivativeSigns,
 230760296                destination.DerivativeMagnitudes, destination.DerivativeSigns, product, accumulate: true);
 297        }
 288180298    }
 299
 300    private static void Add(Polynomial destination, Polynomial source,
 301        int multiplier, int degreeOffset = 0)
 302    {
 172908303        WideFiniteAxisIntersection.AddFiniteAxisPolynomial(
 172908304            source.Magnitudes, source.Signs, destination.Magnitudes, destination.Signs,
 172908305            multiplier, degreeOffset);
 172908306        if (!destination.DerivativeSigns.IsEmpty)
 138456307            WideFiniteAxisIntersection.AddFiniteAxisPolynomial(
 138456308                source.DerivativeMagnitudes, source.DerivativeSigns,
 138456309                destination.DerivativeMagnitudes, destination.DerivativeSigns,
 138456310                multiplier, degreeOffset);
 172908311    }
 312}