| | | 1 | | //======================================================================= |
| | | 2 | | // FiniteSlabRayExtrema.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 | | |
| | | 13 | | |
| | | 14 | | |
| | | 15 | | |
| | | 16 | | |
| | | 17 | | |
| | | 18 | | /// <summary>Finite-domain clearance of one flat-capped planar convex primitive.</summary> |
| | | 19 | | /// <content>Solves shared finite-domain linear and radial active sets.</content> |
| | | 20 | | internal static partial class FiniteSlabRayExtrema |
| | | 21 | | { |
| | | 22 | | // Width ledger: exact rigid points/edges <198 bits; conditioned normal |
| | | 23 | | // components <268; prism linear rows <334; chart radial coordinates <532; |
| | | 24 | | // conic coefficients <1065; PlaneDisk elimination retains quadratic |
| | | 25 | | // coefficients <2133. The single radial constraint produces only |
| | | 26 | | // quadratic roots; reconstructed point coefficients and final depth/sign |
| | | 27 | | // products remain <3400 bits. |
| | | 28 | | // Rigid calls retain the 64-word bank. For exact affine ingress, the |
| | | 29 | | // full clipping rows have <832 bits and coordinates/rays <352. A plane |
| | | 30 | | // chart has <1186-bit coordinates, <2374-bit conics and <4751-bit |
| | | 31 | | // eliminated quadratics. Reconstructed witnesses stay <3563 bits; |
| | | 32 | | // their squared physical metric stays <7834, within 128 words. |
| | | 33 | | private const int Words = FiniteSlabPolynomial.Words; |
| | | 34 | | private readonly struct Linear |
| | | 35 | | { |
| | | 36 | | internal readonly Signed832 A, B, C, D; |
| | | 37 | | internal Linear(Signed832 a, Signed832 b, Signed832 c, Signed832 d) |
| | 32210 | 38 | | { A = a; B = b; C = c; D = d; } |
| | | 39 | | internal Linear(Signed576 a, Signed576 b, Signed576 c, Signed576 d) |
| | 12676 | 40 | | : this(Signed832.ExtendValue(a), Signed832.ExtendValue(b), Signed832.ExtendValue(c), Signed832.ExtendValue(d |
| | 769590 | 41 | | internal Signed832 this[int index] => index == 0 ? A : index == 1 ? B : index == 2 ? C : D; |
| | | 42 | | } |
| | | 43 | | private readonly struct Radial |
| | | 44 | | { |
| | | 45 | | internal readonly WideAxis3 A, E, F; |
| | | 46 | | internal readonly WideDirection3d N; |
| | | 47 | | internal readonly Signed320 Radius; |
| | | 48 | | internal Radial(WideAxis3 a, WideAxis3 e, WideAxis3 f, WideDirection3d n, Signed320 radius) |
| | 5808 | 49 | | { A = a; E = e; F = f; N = n; Radius = radius; } |
| | | 50 | | } |
| | | 51 | | private readonly ref struct Problem |
| | | 52 | | { |
| | | 53 | | internal readonly ReadOnlySpan<Linear> Lines; |
| | | 54 | | internal readonly ReadOnlySpan<Radial> Disks; |
| | | 55 | | internal readonly FixedTriangle Triangle; |
| | | 56 | | internal readonly int Orientation; |
| | | 57 | | internal readonly FiniteSlabRayCandidateVisitor? Visitor; |
| | | 58 | | private readonly ReadOnlySpan<Radial> _shape; |
| | 7714 | 59 | | internal ref readonly Radial Shape => ref _shape[0]; |
| | | 60 | | internal readonly Signed320 HalfThickness, CoreCap; |
| | | 61 | | internal readonly bool HasCore, IsPrism; |
| | | 62 | | internal readonly int ExitLinearStart; |
| | | 63 | | internal readonly int DomainEdgeMask; |
| | | 64 | | internal readonly int ActiveFirst, ActiveSecond, ActiveThird, SolvedDisks; |
| | | 65 | | internal readonly int DomainLineCount; |
| | | 66 | | internal readonly int ExitRegion; |
| | | 67 | | internal readonly int FieldWords; |
| | | 68 | | internal readonly Span<ulong> VertexValues; |
| | | 69 | | internal readonly Span<sbyte> VertexSigns; |
| | | 70 | | internal Problem(ReadOnlySpan<Linear> lines, ReadOnlySpan<Radial> disks, |
| | | 71 | | FixedTriangle triangle, int orientation, FiniteSlabRayCandidateVisitor? visitor, |
| | | 72 | | ReadOnlySpan<Radial> shape, Signed320 halfThickness, Signed320 coreCap, bool hasCore, bool isPrism, int exit |
| | | 73 | | Span<ulong> vertexValues, Span<sbyte> vertexSigns, |
| | | 74 | | int domainEdgeMask = 0, int activeFirst = -1, int activeSecond = -1, int activeThird = -1, int solvedDisks = |
| | | 75 | | int domainLineCount = 3, int exitRegion = 0, int fieldWords = Words) |
| | | 76 | | { |
| | 120770 | 77 | | Lines = lines; Disks = disks; Triangle = triangle; Orientation = orientation; Visitor = visitor; |
| | 96616 | 78 | | _shape = shape; HalfThickness = halfThickness; CoreCap = coreCap; HasCore = hasCore; |
| | 48308 | 79 | | IsPrism = isPrism; ExitLinearStart = exitLinearStart; |
| | 24154 | 80 | | DomainEdgeMask = domainEdgeMask; |
| | 96616 | 81 | | ActiveFirst = activeFirst; ActiveSecond = activeSecond; ActiveThird = activeThird; SolvedDisks = solvedDisks |
| | 24154 | 82 | | DomainLineCount = domainLineCount; |
| | 48308 | 83 | | ExitRegion = exitRegion; FieldWords = fieldWords; |
| | 48308 | 84 | | VertexValues = vertexValues; VertexSigns = vertexSigns; |
| | 24154 | 85 | | } |
| | 19998 | 86 | | internal Problem WithActive(int first, int second = -1, int third = -1) => new(Lines, Disks, Triangle, Orientati |
| | 19998 | 87 | | _shape, HalfThickness, CoreCap, HasCore, IsPrism, ExitLinearStart, VertexValues, VertexSigns, |
| | 19998 | 88 | | DomainEdge(first, DomainLineCount) | DomainEdge(second, DomainLineCount) | DomainEdge(third, DomainLineCount |
| | 19998 | 89 | | first, second, third, SolvedDisks, DomainLineCount, ExitRegion, FieldWords); |
| | 3146 | 90 | | internal Problem WithDisk(int disk) => new(Lines, Disks, Triangle, Orientation, Visitor, |
| | 3146 | 91 | | _shape, HalfThickness, CoreCap, HasCore, IsPrism, ExitLinearStart, VertexValues, VertexSigns, DomainEdgeMask |
| | 3146 | 92 | | ActiveFirst, ActiveSecond, ActiveThird, SolvedDisks | (1 << disk), DomainLineCount, ExitRegion, FieldWords); |
| | 138 | 93 | | internal Problem WithoutSolvedDisks() => new(Lines, Disks, Triangle, Orientation, Visitor, |
| | 138 | 94 | | _shape, HalfThickness, CoreCap, HasCore, IsPrism, ExitLinearStart, VertexValues, VertexSigns, DomainEdgeMask |
| | 138 | 95 | | ActiveFirst, ActiveSecond, ActiveThird, domainLineCount: DomainLineCount, |
| | 138 | 96 | | exitRegion: ExitRegion, fieldWords: FieldWords); |
| | 47789 | 97 | | internal bool IsActive(int index) => index == ActiveFirst || index == ActiveSecond || index == ActiveThird; |
| | | 98 | | } |
| | | 99 | | |
| | | 100 | | /// <summary> |
| | | 101 | | /// Accumulates the maximum nonnegative intersection parameter q+tN, |
| | | 102 | | /// where q belongs to the finite triangle and q+tN belongs to the slab. |
| | | 103 | | /// The source q need not intersect the slab at parameter zero. Empty offsets |
| | | 104 | | /// select a circle/stadium; nonempty offsets select a convex prism. The |
| | | 105 | | /// plane metric belongs to frame, never to the domain triangle's winding |
| | | 106 | | /// or cross-product magnitude. Candidate callbacks are synchronous. |
| | | 107 | | /// </summary> |
| | | 108 | | internal static void AccumulateTriangle(in FiniteSlabRayFrame frame, FixedTriangle triangle, |
| | | 109 | | ReadOnlySpan<Vector2d> prismOffsets, int orientation, ref FiniteSlabRaySelection maximum, |
| | | 110 | | FiniteSlabRayCandidateVisitor? visitor = null) => |
| | 130 | 111 | | AccumulateDomain(frame, triangle, prismOffsets, orientation, true, false, ref maximum, visitor); |
| | | 112 | | |
| | | 113 | | /// <summary>Clips the original segment parameter without rounding either rational cut point.</summary> |
| | | 114 | | internal static void AccumulateSegment(in FiniteSlabRayFrame frame, FixedSegment segment, in RationalSegmentInterval |
| | | 115 | | ReadOnlySpan<Vector2d> prismOffsets, int orientation, ref FiniteSlabRaySelection maximum, |
| | | 116 | | FiniteSlabRayCandidateVisitor? visitor = null) |
| | | 117 | | { |
| | 17 | 118 | | FixedTriangle triangle = new(segment.Start, segment.End, segment.End); |
| | 19 | 119 | | if (!clip.IsClipped) { AccumulateTriangle(frame, triangle, prismOffsets, orientation, ref maximum, visitor); ret |
| | 16 | 120 | | Span<RationalSegmentInterval> clips = stackalloc RationalSegmentInterval[1] { clip }; |
| | 16 | 121 | | AccumulateDomain(frame, triangle, prismOffsets, orientation, true, false, ref maximum, visitor, clips); |
| | 16 | 122 | | } |
| | | 123 | | |
| | | 124 | | /// <summary> |
| | | 125 | | /// Clips physical support charts to a nondegenerate coplanar triangle. |
| | | 126 | | /// Nonaxial faces use their point/line charts; free-height supports also |
| | | 127 | | /// clip against original domain edges. Compare with the unbounded plane |
| | | 128 | | /// peak to certify admission; a lower stationary root is not support. |
| | | 129 | | /// </summary> |
| | | 130 | | internal static void AccumulatePhysicalSupportTriangle(in FiniteSlabRayFrame frame, FixedTriangle triangle, |
| | | 131 | | ReadOnlySpan<Vector2d> prismOffsets, int orientation, ref FiniteSlabRaySelection maximum, |
| | | 132 | | FiniteSlabRayCandidateVisitor? visitor = null) => |
| | 44 | 133 | | AccumulateDomain(frame, triangle, prismOffsets, orientation, true, true, ref maximum, visitor); |
| | | 134 | | |
| | | 135 | | /// <summary> |
| | | 136 | | /// Accumulates the intrinsic slab/plane family without any finite domain |
| | | 137 | | /// inequalities. Hosts must independently certify its retained source |
| | | 138 | | /// projected clearance footprint is covered before using this in a finite |
| | | 139 | | /// mesh domain. Original-section coverage suffices only for axial rays. |
| | | 140 | | /// Axial consumers with representable boundary poles may omit the redundant |
| | | 141 | | /// interior sample; otherwise it remains a potentially sole representable witness. |
| | | 142 | | /// </summary> |
| | | 143 | | internal static void AccumulatePlane(in FiniteSlabRayFrame frame, |
| | | 144 | | ReadOnlySpan<Vector2d> prismOffsets, int orientation, ref FiniteSlabRaySelection maximum, |
| | | 145 | | FiniteSlabRayCandidateVisitor? visitor = null, bool includeInteriorSample = true) => |
| | 355 | 146 | | AccumulateDomain(frame, frame.SourcePlane, prismOffsets, orientation, false, false, ref maximum, visitor, |
| | 355 | 147 | | ReadOnlySpan<RationalSegmentInterval>.Empty, includeInteriorSample); |
| | | 148 | | |
| | | 149 | | private static void AccumulateDomain(in FiniteSlabRayFrame frame, FixedTriangle triangle, |
| | | 150 | | ReadOnlySpan<Vector2d> prismOffsets, int orientation, bool finiteDomain, bool physicalSupportOnly, |
| | | 151 | | ref FiniteSlabRaySelection maximum, FiniteSlabRayCandidateVisitor? visitor) => |
| | 174 | 152 | | AccumulateDomain(frame, triangle, prismOffsets, orientation, finiteDomain, physicalSupportOnly, ref maximum, vis |
| | 174 | 153 | | ReadOnlySpan<RationalSegmentInterval>.Empty); |
| | | 154 | | |
| | | 155 | | private static void AccumulateDomain(in FiniteSlabRayFrame frame, FixedTriangle triangle, |
| | | 156 | | ReadOnlySpan<Vector2d> prismOffsets, int orientation, bool finiteDomain, bool physicalSupportOnly, |
| | | 157 | | ref FiniteSlabRaySelection maximum, FiniteSlabRayCandidateVisitor? visitor, scoped ReadOnlySpan<RationalSegmentI |
| | | 158 | | bool includeInteriorSample = true, scoped ReadOnlySpan<WideAxis3> nativeSource = default, |
| | | 159 | | scoped ReadOnlySpan<Signed832> nativeRows = default) |
| | | 160 | | { |
| | 565 | 161 | | if (orientation != -1 && orientation != 1) throw new ArgumentOutOfRangeException(nameof(orientation)); |
| | 560 | 162 | | if (frame.Plane.NormalSquared.IsZero) return; |
| | 558 | 163 | | WideAxis3 a = nativeSource.IsEmpty ? frame.Section.Transform(triangle.A) : nativeSource[0]; |
| | 558 | 164 | | WideAxis3 e = nativeSource.IsEmpty ? Difference(frame.Section.Transform(triangle.B), a) : nativeSource[1]; |
| | 558 | 165 | | WideAxis3 f = nativeSource.IsEmpty ? Difference(frame.Section.Transform(triangle.C), a) : nativeSource[2]; |
| | 558 | 166 | | WideDirection3d n = Signed(frame.RayNormal, orientation); |
| | 558 | 167 | | bool axial = n.X.IsZero && n.Z.IsZero; |
| | 558 | 168 | | if (!finiteDomain && axial) |
| | | 169 | | { |
| | 31 | 170 | | AccumulateHorizontalPlane(frame, triangle, prismOffsets, orientation, a, e, f, n, ref maximum, visitor, incl |
| | 31 | 171 | | return; |
| | | 172 | | } |
| | 527 | 173 | | int pieces = prismOffsets.IsEmpty && frame.HasCore ? 3 : 1; |
| | 527 | 174 | | int capacity = 10 + prismOffsets.Length + (clips.IsEmpty ? 0 : 2) + nativeRows.Length / 4; |
| | 1054 | 175 | | Span<Linear> lines = stackalloc Linear[capacity]; Span<Radial> disks = stackalloc Radial[1]; |
| | | 176 | | // Active vertices are processed sequentially. Each determinant writes its |
| | | 177 | | // complete bank; retained selections copy their coefficients before reuse. |
| | | 178 | | // Native callers provide a line (dummy mu=0). Its active-set |
| | | 179 | | // coefficients fit the ordinary 4096-bit bank even when the exact |
| | | 180 | | // source denominator requires the section's wider coordinate path. |
| | 527 | 181 | | int fieldWords = nativeSource.IsEmpty ? frame.FieldWords : FiniteSlabPolynomial.Words; |
| | 527 | 182 | | Span<ulong> vertexValues = stackalloc ulong[7 * fieldWords]; |
| | 527 | 183 | | Span<sbyte> vertexSigns = stackalloc sbyte[7]; |
| | | 184 | | // Every active-set view borrows the same immutable shape for this call. |
| | 527 | 185 | | Span<Radial> shape = stackalloc Radial[1] { new(a, e, f, n, frame.Radius) }; |
| | 2736 | 186 | | for (int exitPiece = 0; exitPiece < pieces; exitPiece++) |
| | | 187 | | { |
| | 1682 | 188 | | int count = 0, diskCount = 0; |
| | 841 | 189 | | if (!nativeRows.IsEmpty) |
| | | 190 | | { |
| | 170 | 191 | | for (int row = 0; row < nativeRows.Length; row += 4) |
| | 68 | 192 | | lines[count++] = new(nativeRows[row], nativeRows[row + 1], nativeRows[row + 2], nativeRows[row + 3]) |
| | | 193 | | } |
| | 824 | 194 | | else if (finiteDomain) |
| | | 195 | | { |
| | 300 | 196 | | lines[count++] = new(S(1), default, default, default); |
| | 300 | 197 | | lines[count++] = new(default, S(1), default, default); |
| | 300 | 198 | | lines[count++] = new(S(-1), S(-1), default, S(1)); |
| | | 199 | | } |
| | 841 | 200 | | if (!clips.IsEmpty) |
| | | 201 | | { |
| | 18 | 202 | | ref readonly RationalSegmentInterval clip = ref clips[0]; |
| | | 203 | | // In the degenerate triangle A,B,B the original parameter is |
| | | 204 | | // lambda+mu. Cut numerators/positive denominators are <319 bits, |
| | | 205 | | // below the existing <334-bit prism rows. Active sets still |
| | | 206 | | // use at most three linear rows: LP determinants <1005 bits, |
| | | 207 | | // radial chart coefficients <1065, eliminated quadratics |
| | | 208 | | // <2133 and final witness/sign products <3400 fit the existing |
| | | 209 | | // 4096-bit retained fields. Constraint count adds no width. |
| | 18 | 210 | | lines[count++] = new(clip.WideLowerDenominator, clip.WideLowerDenominator, default, Negate(clip.WideLowe |
| | 18 | 211 | | lines[count++] = new(Negate(clip.WideUpperDenominator), Negate(clip.WideUpperDenominator), default, clip |
| | | 212 | | } |
| | 841 | 213 | | int domainLineCount = count; |
| | 841 | 214 | | lines[count++] = new(default, default, S(1), default); |
| | 841 | 215 | | int exitLinearStart = count; |
| | | 216 | | // Clearance concerns every finite-domain source point that can |
| | | 217 | | // meet the translated slab. Restricting q to the original slab |
| | | 218 | | // misses points entering while it moves, and under-corrects even |
| | | 219 | | // an ordinary tall cylinder against a complete tilted plane. |
| | 841 | 220 | | AddInterval(lines, ref count, a.Y, e.Y, f.Y, n.Y, frame.Section.Cap); |
| | 841 | 221 | | if (!prismOffsets.IsEmpty) |
| | | 222 | | { |
| | 192 | 223 | | int winding = FiniteSlabSectionFrame.GetPrismWinding(prismOffsets); |
| | 192 | 224 | | Signed320 scale = WideArithmetic.AddSigned320(frame.Section.RadialFrame.ShapeFrame.CoordinateDenominator |
| | 192 | 225 | | frame.Section.RadialFrame.ShapeFrame.CoordinateDenominator); |
| | 1824 | 226 | | for (int edge = 0; edge < prismOffsets.Length; edge++) |
| | | 227 | | { |
| | 1440 | 228 | | Vector2d first = prismOffsets[edge], last = prismOffsets[(edge + 1) % prismOffsets.Length]; |
| | 720 | 229 | | Signed320 dx = Signed320.ExtendValue(WideArithmetic.SubtractSigned192(Signed192.Raw(last.X), Signed1 |
| | 720 | 230 | | Signed320 dz = Signed320.ExtendValue(WideArithmetic.SubtractSigned192(Signed192.Raw(last.Y), Signed1 |
| | 720 | 231 | | WideAxis3 offset = new(WideArithmetic.SubtractSigned320(a.X, ScaledOffset(first.X, scale)), |
| | 720 | 232 | | a.Y, WideArithmetic.SubtractSigned320(a.Z, ScaledOffset(first.Y, scale))); |
| | 720 | 233 | | lines[count++] = new(Signed(Wall(e, dx, dz), winding), Signed(Wall(f, dx, dz), winding), |
| | 720 | 234 | | Signed(Wall(n, dx, dz), winding), Signed(Wall(offset, dx, dz), winding)); |
| | | 235 | | } |
| | | 236 | | } |
| | 649 | 237 | | else if (frame.HasCore && exitPiece == 0) |
| | | 238 | | { |
| | 157 | 239 | | AddInterval(lines, ref count, a.X, e.X, f.X, n.X, frame.Radius); |
| | 157 | 240 | | AddInterval(lines, ref count, a.Z, e.Z, f.Z, n.Z, frame.CoreCap); |
| | | 241 | | } |
| | | 242 | | else |
| | | 243 | | { |
| | 492 | 244 | | Signed320 center = frame.HasCore ? Signed(frame.CoreCap, exitPiece == 1 ? -1 : 1) : default; |
| | 492 | 245 | | WideAxis3 offset = new(a.X, a.Y, WideArithmetic.SubtractSigned320(a.Z, center)); |
| | 492 | 246 | | if (frame.Radius.IsZero) |
| | | 247 | | { |
| | 82 | 248 | | AddInterval(lines, ref count, offset.X, e.X, f.X, n.X, default); |
| | 82 | 249 | | AddInterval(lines, ref count, offset.Z, e.Z, f.Z, n.Z, default); |
| | | 250 | | } |
| | 410 | 251 | | else disks[diskCount++] = new(offset, e, f, n, frame.Radius); |
| | | 252 | | } |
| | 841 | 253 | | Problem problem = new(lines[..count], disks[..diskCount], triangle, orientation, visitor, |
| | 841 | 254 | | shape, frame.Section.Cap, frame.CoreCap, frame.HasCore, |
| | 841 | 255 | | !prismOffsets.IsEmpty, exitLinearStart, vertexValues, vertexSigns, domainLineCount: domainLineCount, |
| | 841 | 256 | | exitRegion: frame.HasCore && exitPiece != 0 ? exitPiece == 1 ? -1 : 1 : 0, fieldWords: fieldWords); |
| | 841 | 257 | | if (!finiteDomain || physicalSupportOnly && !axial) |
| | 588 | 258 | | AccumulatePhysicalSupport(problem, frame, prismOffsets, exitPiece, ref maximum); |
| | 301 | 259 | | else if (axial) AccumulateAxial(problem, n.Y.Sign, ref maximum); |
| | 205 | 260 | | else Accumulate(problem, ref maximum); |
| | | 261 | | } |
| | 527 | 262 | | } |
| | | 263 | | |
| | | 264 | | private static void AccumulatePhysicalSupport(scoped Problem problem, in FiniteSlabRayFrame frame, |
| | | 265 | | ReadOnlySpan<Vector2d> offsets, int exitPiece, ref FiniteSlabRaySelection maximum) |
| | | 266 | | { |
| | 588 | 267 | | WideDirection3d n = problem.Shape.N; |
| | 588 | 268 | | if (offsets.IsEmpty && frame.HasCore) |
| | | 269 | | { |
| | 372 | 270 | | if (exitPiece == 0 && !n.Z.IsZero) return; |
| | 336 | 271 | | if (exitPiece != 0 && !n.Z.IsZero && n.Z.Sign != problem.ExitRegion) return; |
| | 284 | 272 | | if (problem.DomainLineCount != 0 && n.Z.IsZero && exitPiece != 0) return; |
| | | 273 | | } |
| | | 274 | | // Without a finite domain, support of the extrusion is the product |
| | | 275 | | // of an outward Y cap and footprint support; free Y retains both caps. |
| | 496 | 276 | | int firstCap = problem.ExitLinearStart + (n.Y.Sign > 0 ? 1 : 0); |
| | 496 | 277 | | int lastCap = n.Y.IsZero ? firstCap + 1 : firstCap; |
| | 496 | 278 | | Signed576 footprintMaximum = default; |
| | 992 | 279 | | int firstVertex = 0, lastVertex = 0; |
| | 496 | 280 | | if (!offsets.IsEmpty) |
| | 1464 | 281 | | for (int vertex = 0; vertex < offsets.Length; vertex++) |
| | | 282 | | { |
| | 576 | 283 | | Signed576 projection = FootprintProjection(n, offsets[vertex]); |
| | 576 | 284 | | int order = WideArithmetic.SubtractSigned576(projection, footprintMaximum).Sign; |
| | 576 | 285 | | if (vertex == 0 || order > 0) |
| | 524 | 286 | | { footprintMaximum = projection; firstVertex = lastVertex = vertex; } |
| | 314 | 287 | | else if (order == 0) |
| | | 288 | | { |
| | 342 | 289 | | Vector2d offset = offsets[vertex], first = offsets[firstVertex], last = offsets[lastVertex]; |
| | 138 | 290 | | if (offset.X < first.X || offset.X == first.X && offset.Y < first.Y) firstVertex = vertex; |
| | 174 | 291 | | if (offset.X > last.X || offset.X == last.X && offset.Y > last.Y) lastVertex = vertex; |
| | | 292 | | } |
| | | 293 | | } |
| | 496 | 294 | | int supportWall = -1; |
| | | 295 | | // A free-height support can end on a triangle edge before either Y cap. |
| | | 296 | | // Keep those original domain constraints as well as the two real caps; |
| | | 297 | | // the same point/line charts then clip the entire physical support face. |
| | 496 | 298 | | int domainChoices = n.Y.IsZero ? problem.DomainLineCount : 0; |
| | 496 | 299 | | if (!offsets.IsEmpty && problem.DomainLineCount != 0 && firstVertex != lastVertex) |
| | 120 | 300 | | for (int edge = 0; edge < offsets.Length; edge++) |
| | | 301 | | { |
| | 120 | 302 | | Vector2d first = offsets[edge], last = offsets[(edge + 1) % offsets.Length]; |
| | 60 | 303 | | if (first == last) continue; |
| | 60 | 304 | | if (WideArithmetic.SubtractSigned576(FootprintProjection(n, first), footprintMaximum).Sign != 0 |
| | 60 | 305 | | || WideArithmetic.SubtractSigned576(FootprintProjection(n, last), footprintMaximum).Sign != 0) conti |
| | 20 | 306 | | supportWall = problem.ExitLinearStart + 2 + edge; |
| | 20 | 307 | | break; |
| | | 308 | | } |
| | 2476 | 309 | | for (int choice = 0; choice < domainChoices + lastCap - firstCap + 1; choice++) |
| | | 310 | | { |
| | 742 | 311 | | int capIndex = choice < domainChoices ? choice : firstCap + choice - domainChoices; |
| | 742 | 312 | | ref readonly Linear cap = ref problem.Lines[capIndex]; |
| | 742 | 313 | | if (!offsets.IsEmpty) |
| | | 314 | | { |
| | 228 | 315 | | if (supportWall >= 0) |
| | | 316 | | { |
| | 20 | 317 | | AccumulateSupportLine(problem, capIndex, supportWall, ref maximum); |
| | 20 | 318 | | continue; |
| | | 319 | | } |
| | 208 | 320 | | Signed320 scale = WideArithmetic.AddSigned320(frame.Section.RadialFrame.ShapeFrame.CoordinateDenominator |
| | 208 | 321 | | frame.Section.RadialFrame.ShapeFrame.CoordinateDenominator); |
| | | 322 | | // A convex footprint's tied support is one straight face; |
| | | 323 | | // its lexical endpoints omit duplicate/collinear authored samples. |
| | 972 | 324 | | for (int vertex = 0; vertex < (firstVertex == lastVertex ? 1 : 2); vertex++) |
| | | 325 | | { |
| | 278 | 326 | | Vector2d offset = offsets[vertex == 0 ? firstVertex : lastVertex]; |
| | 278 | 327 | | Linear x = new(S(problem.Shape.E.X), S(problem.Shape.F.X), S(n.X), |
| | 278 | 328 | | S(WideArithmetic.SubtractSigned320(problem.Shape.A.X, |
| | 278 | 329 | | ScaledOffset(offset.X, scale)))); |
| | 278 | 330 | | Linear z = new(S(problem.Shape.E.Z), S(problem.Shape.F.Z), S(n.Z), |
| | 278 | 331 | | S(WideArithmetic.SubtractSigned320(problem.Shape.A.Z, |
| | 278 | 332 | | ScaledOffset(offset.Y, scale)))); |
| | | 333 | | // Direct coordinates retain authored collinear/duplicate |
| | | 334 | | // corners whose adjacent wall equations are dependent. |
| | 278 | 335 | | ThreePlanes(problem.WithActive(capIndex), cap, x, z, ref maximum); |
| | | 336 | | } |
| | | 337 | | } |
| | 514 | 338 | | else if (frame.HasCore && exitPiece == 0) |
| | | 339 | | { |
| | 108 | 340 | | int xIndex = problem.ExitLinearStart + 2 + (n.X.Sign > 0 ? 1 : 0); |
| | 108 | 341 | | if (problem.DomainLineCount != 0) |
| | | 342 | | { |
| | 10 | 343 | | AccumulateSupportLine(problem, capIndex, xIndex, ref maximum); |
| | 10 | 344 | | continue; |
| | | 345 | | } |
| | 588 | 346 | | for (int zIndex = problem.ExitLinearStart + 4; zIndex < problem.ExitLinearStart + 6; zIndex++) |
| | 196 | 347 | | ThreePlanes(problem.WithActive(capIndex, xIndex, zIndex), |
| | 196 | 348 | | cap, problem.Lines[xIndex], problem.Lines[zIndex], ref maximum); |
| | | 349 | | } |
| | 406 | 350 | | else if (frame.Radius.IsZero) |
| | | 351 | | { |
| | 104 | 352 | | int xIndex = problem.ExitLinearStart + 2, zIndex = problem.ExitLinearStart + 4; |
| | 52 | 353 | | ThreePlanes(problem.WithActive(capIndex, xIndex, zIndex), |
| | 52 | 354 | | cap, problem.Lines[xIndex], problem.Lines[zIndex], ref maximum); |
| | | 355 | | } |
| | | 356 | | else |
| | | 357 | | { |
| | 354 | 358 | | Radial disk = problem.Disks[0]; |
| | | 359 | | // The existing cap-constrained stationary chart supplies the |
| | | 360 | | // radial support pole without enumerating other active strata. |
| | 354 | 361 | | PlaneDisk(problem.WithActive(capIndex).WithDisk(0), cap, disk, ref maximum); |
| | | 362 | | } |
| | | 363 | | } |
| | 496 | 364 | | } |
| | | 365 | | |
| | | 366 | | private static void AccumulateSupportLine(scoped Problem problem, int cap, int wall, |
| | | 367 | | ref FiniteSlabRaySelection maximum) |
| | | 368 | | { |
| | | 369 | | // Cap and support wall fix a constant-depth line. Its clipped ends |
| | | 370 | | // occur on an existing footprint/domain constraint, including tangency. |
| | 660 | 371 | | for (int third = 0; third < problem.Lines.Length; third++) |
| | | 372 | | { |
| | 300 | 373 | | if (third == cap || third == wall) continue; |
| | 240 | 374 | | ThreePlanes(problem.WithActive(cap, wall, third), problem.Lines[cap], |
| | 240 | 375 | | problem.Lines[wall], problem.Lines[third], ref maximum); |
| | | 376 | | } |
| | 30 | 377 | | } |
| | | 378 | | |
| | 666 | 379 | | private static Signed576 FootprintProjection(WideDirection3d n, Vector2d offset) => WideArithmetic.AddSigned576( |
| | 666 | 380 | | WideArithmetic.MultiplySigned576(n.X, offset.X.m_rawValue), |
| | 666 | 381 | | WideArithmetic.MultiplySigned576(n.Z, offset.Y.m_rawValue)); |
| | | 382 | | |
| | | 383 | | private static void AccumulateHorizontalPlane(in FiniteSlabRayFrame frame, FixedTriangle triangle, |
| | | 384 | | ReadOnlySpan<Vector2d> offsets, int orientation, WideAxis3 a, WideAxis3 e, WideAxis3 f, WideDirection3d n, |
| | | 385 | | ref FiniteSlabRaySelection maximum, FiniteSlabRayCandidateVisitor? visitor, bool includeInteriorSample) |
| | | 386 | | { |
| | | 387 | | // The intrinsic horizontal section is the complete footprint at one |
| | | 388 | | // Y coordinate. Every point has the same rational clearance depth; |
| | | 389 | | // canonical footprint poles/corners cover the constant-depth family |
| | | 390 | | // directly, without solving redundant circle/line quadratics. |
| | 62 | 391 | | Span<Linear> lines = stackalloc Linear[3]; int count = 0; |
| | 34 | 392 | | if (!offsets.IsEmpty) _ = FiniteSlabSectionFrame.GetPrismWinding(offsets); |
| | 31 | 393 | | lines[count++] = new(default, default, S(1), default); |
| | 31 | 394 | | AddInterval(lines, ref count, a.Y, e.Y, f.Y, n.Y, frame.Section.Cap); |
| | 31 | 395 | | int capIndex = 1 + (n.Y.Sign > 0 ? 1 : 0); |
| | 31 | 396 | | Span<Radial> shape = stackalloc Radial[1] { new(a, e, f, n, frame.Radius) }; |
| | 31 | 397 | | Span<ulong> vertexValues = stackalloc ulong[7 * frame.FieldWords]; |
| | 31 | 398 | | Span<sbyte> vertexSigns = stackalloc sbyte[7]; |
| | 31 | 399 | | Problem problem = new(lines, ReadOnlySpan<Radial>.Empty, triangle, orientation, visitor, |
| | 31 | 400 | | shape, frame.Section.Cap, frame.CoreCap, frame.HasCore, |
| | 31 | 401 | | !offsets.IsEmpty, 1, vertexValues, vertexSigns, activeFirst: capIndex, domainLineCount: 0, fieldWords: frame |
| | 31 | 402 | | int samples = !offsets.IsEmpty ? offsets.Length : (frame.HasCore ? 6 : 4) + (includeInteriorSample ? 1 : 0); |
| | 31 | 403 | | Signed320 scale = WideArithmetic.AddSigned320(frame.Section.RadialFrame.ShapeFrame.CoordinateDenominator, |
| | 31 | 404 | | frame.Section.RadialFrame.ShapeFrame.CoordinateDenominator); |
| | 376 | 405 | | for (int index = 0; index < samples; index++) |
| | | 406 | | { |
| | 314 | 407 | | Signed320 x = default, z = default; |
| | 157 | 408 | | if (!offsets.IsEmpty) |
| | | 409 | | { |
| | 12 | 410 | | x = ScaledOffset(offsets[index].X, scale); |
| | 12 | 411 | | z = ScaledOffset(offsets[index].Y, scale); |
| | | 412 | | } |
| | 145 | 413 | | else if (frame.HasCore) |
| | | 414 | | { |
| | 33 | 415 | | if (index < 4) |
| | 40 | 416 | | { x = Signed(frame.Radius, index % 2 == 0 ? -1 : 1); z = Signed(frame.CoreCap, index < 2 ? -1 : 1); } |
| | 13 | 417 | | else if (index < 6) |
| | 10 | 418 | | z = Signed(WideArithmetic.AddSigned320(frame.CoreCap, frame.Radius), index == 4 ? -1 : 1); |
| | | 419 | | } |
| | 112 | 420 | | else if (index < 4) |
| | | 421 | | { |
| | 138 | 422 | | if (index < 2) x = Signed(frame.Radius, index == 0 ? -1 : 1); |
| | 46 | 423 | | else z = Signed(frame.Radius, index == 2 ? -1 : 1); |
| | | 424 | | } |
| | 157 | 425 | | Linear qx = new(S(e.X), S(f.X), default, S(WideArithmetic.SubtractSigned320(a.X, x))); |
| | 157 | 426 | | Linear qz = new(S(e.Z), S(f.Z), default, S(WideArithmetic.SubtractSigned320(a.Z, z))); |
| | 157 | 427 | | ThreePlanes(problem, lines[capIndex], qx, qz, ref maximum); |
| | | 428 | | } |
| | 31 | 429 | | } |
| | | 430 | | |
| | | 431 | | private static void AddInterval(Span<Linear> lines, ref int count, |
| | | 432 | | Signed320 a, Signed320 e, Signed320 f, Signed576 n, Signed320 cap) |
| | | 433 | | { |
| | 1350 | 434 | | lines[count++] = new(S(e), S(f), S(n), S(WideArithmetic.AddSigned320(a, cap))); |
| | 1350 | 435 | | lines[count++] = new(S(Negate(e)), S(Negate(f)), S(Signed(n, -1)), S(WideArithmetic.SubtractSigned320(cap, a))); |
| | 1350 | 436 | | } |
| | | 437 | | |
| | | 438 | | private static void Accumulate(scoped Problem problem, ref FiniteSlabRaySelection maximum) |
| | | 439 | | { |
| | | 440 | | // There is at most one radial quadratic. An isolated maximum has |
| | | 441 | | // three active planes, two planes plus that quadratic, or one plane |
| | | 442 | | // plus its constrained stationary point. A free radial family has a |
| | | 443 | | // null direction preserving depth; the finite Y caps/domain force it |
| | | 444 | | // to one of these real planes, so no separate zero-plane chart is needed. |
| | 3788 | 445 | | for (int first = 0; first < problem.Lines.Length; first++) |
| | | 446 | | { |
| | 1689 | 447 | | ref readonly Linear a = ref problem.Lines[first]; |
| | 16394 | 448 | | for (int second = first + 1; second < problem.Lines.Length; second++) |
| | | 449 | | { |
| | 6508 | 450 | | ref readonly Linear b = ref problem.Lines[second]; |
| | 44126 | 451 | | for (int third = second + 1; third < problem.Lines.Length; third++) |
| | 15555 | 452 | | ThreePlanes(problem.WithActive(first, second, third), |
| | 15555 | 453 | | a, b, problem.Lines[third], ref maximum); |
| | 16294 | 454 | | for (int disk = 0; disk < problem.Disks.Length; disk++) |
| | 1639 | 455 | | TwoPlanesDisk(problem.WithActive(first, second).WithDisk(disk), |
| | 1639 | 456 | | a, b, problem.Disks[disk], ref maximum); |
| | | 457 | | } |
| | 4608 | 458 | | for (int disk = 0; disk < problem.Disks.Length; disk++) |
| | 615 | 459 | | PlaneDisk(problem.WithActive(first).WithDisk(disk), a, problem.Disks[disk], ref maximum); |
| | | 460 | | } |
| | 205 | 461 | | } |
| | | 462 | | |
| | | 463 | | private static void AccumulateAxial(scoped Problem problem, int direction, |
| | | 464 | | ref FiniteSlabRaySelection maximum) |
| | | 465 | | { |
| | | 466 | | // Every feasible ray leaves through the outward cap. Enumerating only |
| | | 467 | | // active sets containing it retains all clipped-region extrema and |
| | | 468 | | // cardinal witnesses, without solving irrelevant radial exit charts. |
| | 48 | 469 | | int capIndex = problem.ExitLinearStart + (direction > 0 ? 1 : 0); |
| | 48 | 470 | | ref readonly Linear cap = ref problem.Lines[capIndex]; |
| | 780 | 471 | | for (int first = 0; first < problem.Lines.Length; first++) |
| | | 472 | | { |
| | 342 | 473 | | if (first == capIndex) continue; |
| | 2688 | 474 | | for (int second = first + 1; second < problem.Lines.Length; second++) |
| | | 475 | | { |
| | 1050 | 476 | | if (second == capIndex) continue; |
| | 807 | 477 | | ThreePlanes(problem.WithActive(capIndex, first, second), |
| | 807 | 478 | | cap, problem.Lines[first], problem.Lines[second], ref maximum); |
| | | 479 | | } |
| | 1032 | 480 | | for (int disk = 0; disk < problem.Disks.Length; disk++) |
| | 222 | 481 | | TwoPlanesDisk(problem.WithActive(capIndex, first).WithDisk(disk), cap, problem.Lines[first], problem.Dis |
| | | 482 | | } |
| | 176 | 483 | | for (int disk = 0; disk < problem.Disks.Length; disk++) |
| | 40 | 484 | | PlaneDisk(problem.WithActive(capIndex).WithDisk(disk), cap, problem.Disks[disk], ref maximum); |
| | 48 | 485 | | } |
| | | 486 | | |
| | | 487 | | private static void ThreePlanes(scoped Problem problem, in Linear a, in Linear b, in Linear c, |
| | | 488 | | ref FiniteSlabRaySelection maximum) |
| | | 489 | | { |
| | 17423 | 490 | | int words = problem.FieldWords; |
| | 17423 | 491 | | Span<ulong> values = problem.VertexValues; |
| | 17423 | 492 | | Span<sbyte> signs = problem.VertexSigns; |
| | 34846 | 493 | | FiniteSlabPolynomial l = P(values, signs, 0), m = P(values, signs, 1), |
| | 69692 | 494 | | t = P(values, signs, 2), d = P(values, signs, 3), n = P(values, signs, 4), u = P(values, signs, 5, 2); |
| | 17423 | 495 | | Determinant(d, a, b, c, 0, 1, 2); |
| | 26402 | 496 | | if (d.IsZero) return; |
| | 16888 | 497 | | Determinant(l, a, b, c, 3, 1, 2); l.Signs[0] *= -1; |
| | 16888 | 498 | | Determinant(m, a, b, c, 0, 3, 2); m.Signs[0] *= -1; |
| | 16888 | 499 | | Determinant(t, a, b, c, 0, 1, 3); t.Signs[0] *= -1; |
| | 25496 | 500 | | if (d.Signs[0] < 0) { l.Signs[0] *= -1; m.Signs[0] *= -1; t.Signs[0] *= -1; d.Signs[0] = 1; } |
| | 10384 | 501 | | if (t.Signs[0] < 0) return; |
| | | 502 | | // The LP witness is rational. Its source has no algebraic parameter; |
| | | 503 | | // retain a known half-unit root and the independent rational depth. |
| | 6504 | 504 | | WideFiniteAxisIntersection.NormalizeFiniteAxisPolynomialPowerOfTwo(values[..(4 * words)], signs[..4]); |
| | 32520 | 505 | | n.Copy(t); n.Add(d); u.Clear(); u.Add(t); u.Add(d, 2); |
| | 6504 | 506 | | ProcessRational(problem, l, m, t, d, n, u.Slice(0, 1), ref maximum); |
| | 6504 | 507 | | } |
| | | 508 | | |
| | | 509 | | private static void Determinant(FiniteSlabPolynomial result, in Linear a, in Linear b, in Linear c, |
| | | 510 | | int x, int y, int z) |
| | | 511 | | { |
| | | 512 | | // Rigid rows retain their original 576-bit factors and scratch. Only |
| | | 513 | | // affine coordinates or wide clips require complete 832-bit factors. |
| | 42755 | 514 | | bool compact = result.FieldWords == Words; |
| | 42755 | 515 | | int limbs = compact ? 9 : 13; |
| | 42755 | 516 | | Span<ulong> values = compact ? stackalloc ulong[72] : stackalloc ulong[104]; |
| | 128265 | 517 | | Span<ulong> first = values[..limbs], second = values.Slice(limbs, limbs), third = values.Slice(2 * limbs, limbs) |
| | 85510 | 518 | | pair = values.Slice(3 * limbs, 2 * limbs), term = values.Slice(5 * limbs, 3 * limbs); |
| | 42755 | 519 | | result.Clear(); |
| | 42755 | 520 | | int sign = 0; |
| | 598570 | 521 | | for (int permutation = 0; permutation < 6; permutation++) |
| | | 522 | | { |
| | 256530 | 523 | | int i = permutation % 3; |
| | 256530 | 524 | | int j = permutation < 3 ? (i + 1) % 3 : (i + 2) % 3; |
| | 256530 | 525 | | int k = 3 - i - j; |
| | 513060 | 526 | | Signed832 av = a[i == 0 ? x : i == 1 ? y : z], bv = b[j == 0 ? x : j == 1 ? y : z], |
| | 256530 | 527 | | cv = c[k == 0 ? x : k == 1 ? y : z]; |
| | | 528 | | // Sparse cap/footprint rows have many exact zero products; they |
| | | 529 | | // need no polynomial scratch initialization or multiplication. |
| | 256530 | 530 | | if (av.IsZero || bv.IsZero || cv.IsZero) continue; |
| | 47107 | 531 | | if (compact) |
| | | 532 | | { |
| | 46325 | 533 | | WideArithmetic.GetMagnitude(Signed576.NarrowValue(av), first); |
| | 46325 | 534 | | WideArithmetic.GetMagnitude(Signed576.NarrowValue(bv), second); |
| | 46325 | 535 | | WideArithmetic.GetMagnitude(Signed576.NarrowValue(cv), third); |
| | | 536 | | } |
| | | 537 | | else |
| | | 538 | | { |
| | 2346 | 539 | | WideArithmetic.GetMagnitude(av, first); WideArithmetic.GetMagnitude(bv, second); WideArithmetic.GetMagni |
| | | 540 | | } |
| | 47107 | 541 | | WideArithmetic.MultiplyMagnitudes(first, second, pair); |
| | 47107 | 542 | | WideArithmetic.MultiplyMagnitudes(pair, third, term); |
| | 47107 | 543 | | WideArithmetic.AddShiftedSignedMagnitude(term, av.Sign * bv.Sign * cv.Sign * (permutation < 3 ? 1 : -1), |
| | 47107 | 544 | | 0, result.Values, ref sign); |
| | | 545 | | } |
| | 42755 | 546 | | result.Signs[0] = (sbyte)sign; |
| | 42755 | 547 | | } |
| | | 548 | | |
| | | 549 | | private static void Process(scoped Problem problem, scoped FiniteSlabPolynomial equation, |
| | | 550 | | scoped FiniteSlabPolynomial l, scoped FiniteSlabPolynomial m, scoped FiniteSlabPolynomial t, |
| | | 551 | | scoped FiniteSlabPolynomial d, bool parameterIsDepth, |
| | | 552 | | scoped FiniteSlabPolynomial depthN, scoped FiniteSlabPolynomial depthD, |
| | | 553 | | ref FiniteSlabRaySelection maximum) |
| | | 554 | | { |
| | 3154 | 555 | | int words = problem.FieldWords; |
| | 3237 | 556 | | if (equation.IsZero) return; |
| | 3071 | 557 | | WideFiniteAxisIntersection.NormalizeFiniteAxisPolynomialPowerOfTwo(equation.Values, equation.Signs); |
| | 6142 | 558 | | Span<ulong> cells = stackalloc ulong[64]; Span<int> shifts = stackalloc int[8]; |
| | 3071 | 559 | | FiniteAxisValueRoots roots = WideFiniteAxisIntersection.GetFiniteValueRoots(equation.Values, equation.Signs, cel |
| | 3071 | 560 | | Span<ulong> cell = stackalloc ulong[WideFiniteAxisIntersection.GetFiniteValueRootCellWords(equation.Values, equa |
| | 3071 | 561 | | Span<ulong> candidateValues = stackalloc ulong[problem.Visitor is null ? 0 : FiniteSlabRaySelection.GetStorageWo |
| | 3071 | 562 | | Span<sbyte> candidateSigns = stackalloc sbyte[problem.Visitor is null ? 0 : FiniteSlabRaySelection.SignCount]; |
| | 3071 | 563 | | Span<ulong> half = stackalloc ulong[1] { 1 }; |
| | 11096 | 564 | | for (int ordinal = 0; ordinal < roots.Count; ordinal++) |
| | | 565 | | { |
| | 2477 | 566 | | FiniteAxisValueRoot root = FiniteAxisValueRoots.GetRoot(roots, ordinal, equation.Values, equation.Signs, cel |
| | | 567 | | // The depth chart maps t>=0 to u>=1/2. Reject its negative |
| | | 568 | | // branch using the already isolated cell before full admission. |
| | 2477 | 569 | | if (parameterIsDepth && WideFiniteAxisIntersection.CompareFiniteValueRootToDyadic(root, half, 1) < 0) contin |
| | 1752 | 570 | | if (!Admitted(problem, ref root, l, m, t, d)) continue; |
| | 895 | 571 | | maximum.Keep(root, roots.Count, l, m, t, d, parameterIsDepth, depthN, depthD, |
| | 895 | 572 | | problem.Triangle, problem.Orientation, problem.DomainEdgeMask); |
| | 895 | 573 | | if (problem.Visitor is not null && IsClearanceExit(problem, ref root, l, m, t, d)) |
| | | 574 | | { |
| | 187 | 575 | | FiniteSlabRaySelection candidate = new(candidateValues, candidateSigns, words); |
| | 187 | 576 | | candidate.Keep(root, roots.Count, l, m, t, d, parameterIsDepth, depthN, depthD, |
| | 187 | 577 | | problem.Triangle, problem.Orientation, problem.DomainEdgeMask); |
| | 187 | 578 | | problem.Visitor(candidate); |
| | | 579 | | } |
| | | 580 | | } |
| | 3071 | 581 | | } |
| | | 582 | | |
| | | 583 | | private static void ProcessRational(scoped Problem problem, |
| | | 584 | | scoped FiniteSlabPolynomial l, scoped FiniteSlabPolynomial m, scoped FiniteSlabPolynomial t, |
| | | 585 | | scoped FiniteSlabPolynomial d, scoped FiniteSlabPolynomial depthN, scoped FiniteSlabPolynomial depthD, |
| | | 586 | | ref FiniteSlabRaySelection maximum) |
| | | 587 | | { |
| | 6504 | 588 | | int words = problem.FieldWords; |
| | 13008 | 589 | | Span<ulong> values = stackalloc ulong[2 * words]; Span<sbyte> signs = stackalloc sbyte[2]; |
| | 26016 | 590 | | FiniteSlabPolynomial equation = new(values, signs); equation.Clear(); equation.Set(-1); equation.Set(2, 1); |
| | 13008 | 591 | | Span<ulong> cell = stackalloc ulong[1]; cell[0] = 1; |
| | 6504 | 592 | | FiniteAxisValueRoot root = new() { Coefficients = values, Signs = signs, |
| | 6504 | 593 | | LowerNumerator = cell, DenominatorShift = 1, IsRational = true, Ordinal = 0 }; |
| | 11033 | 594 | | if (!Admitted(problem, ref root, l, m, t, d)) return; |
| | 1975 | 595 | | maximum.Keep(root, 1, l, m, t, d, false, depthN, depthD, |
| | 1975 | 596 | | problem.Triangle, problem.Orientation, problem.DomainEdgeMask); |
| | 3440 | 597 | | if (problem.Visitor is null || !IsClearanceExit(problem, ref root, l, m, t, d)) return; |
| | 510 | 598 | | Span<ulong> candidateValues = stackalloc ulong[FiniteSlabRaySelection.GetStorageWords(words)]; |
| | 510 | 599 | | Span<sbyte> candidateSigns = stackalloc sbyte[FiniteSlabRaySelection.SignCount]; |
| | 510 | 600 | | FiniteSlabRaySelection candidate = new(candidateValues, candidateSigns, words); |
| | 510 | 601 | | candidate.Keep(root, 1, l, m, t, d, false, depthN, depthD, |
| | 510 | 602 | | problem.Triangle, problem.Orientation, problem.DomainEdgeMask); |
| | 510 | 603 | | problem.Visitor(candidate); |
| | 510 | 604 | | } |
| | | 605 | | |
| | | 606 | | private static bool Admitted(scoped Problem problem, scoped ref FiniteAxisValueRoot root, |
| | | 607 | | scoped FiniteSlabPolynomial l, scoped FiniteSlabPolynomial m, |
| | | 608 | | scoped FiniteSlabPolynomial t, scoped FiniteSlabPolynomial d) |
| | | 609 | | { |
| | 8256 | 610 | | int words = problem.FieldWords; |
| | 8256 | 611 | | int denominatorSign = FiniteSlabPolynomial.Sign(ref root, d); |
| | | 612 | | System.Diagnostics.Debug.Assert(denominatorSign != 0); |
| | 16512 | 613 | | Span<ulong> values = stackalloc ulong[10 * words]; Span<sbyte> signs = stackalloc sbyte[10]; |
| | 16512 | 614 | | FiniteSlabPolynomial q = P(values, signs, 0, 5), term = P(values, signs, 5, 5); |
| | 8256 | 615 | | if (problem.ExitRegion != 0) |
| | | 616 | | { |
| | | 617 | | // End disks are admitted only on their genuine hemispheres. |
| | | 618 | | // Their cut boundaries z=+/-core already belong to the central |
| | | 619 | | // rectangle's linear strata, so no extra chart is required. |
| | | 620 | | // This exact partition removes covered disk families and their |
| | | 621 | | // duplicate root/depth comparisons without changing the union. |
| | 2775 | 622 | | Coordinate(problem.Shape, 2, l, m, t, d, q, term); |
| | 5550 | 623 | | term.Scale(d, Signed(problem.CoreCap, problem.ExitRegion)); q.Add(term, -1); |
| | 3847 | 624 | | if (FiniteSlabPolynomial.Sign(ref root, q) * denominatorSign * problem.ExitRegion < 0) return false; |
| | | 625 | | } |
| | 97312 | 626 | | for (int index = 0; index < problem.Lines.Length; index++) |
| | | 627 | | { |
| | | 628 | | // Defining planes are exact algebraic identities of this chart. |
| | | 629 | | // Rechecking them would invoke the full resultant equality bound. |
| | 45471 | 630 | | if (problem.IsActive(index)) continue; |
| | 31281 | 631 | | ref readonly Linear line = ref problem.Lines[index]; |
| | 156405 | 632 | | q.Scale(l, line.A); term.Scale(m, line.B); q.Add(term); term.Scale(t, line.C); q.Add(term); |
| | 62562 | 633 | | term.Scale(d, line.D); q.Add(term); |
| | 35280 | 634 | | if (FiniteSlabPolynomial.Sign(ref root, q) * denominatorSign < 0) return false; |
| | | 635 | | } |
| | | 636 | | // Prisms and already separated candidates never need radial scratch. |
| | 4958 | 637 | | if (problem.Disks.IsEmpty) return true; |
| | 2824 | 638 | | Span<ulong> radialValues = stackalloc ulong[24 * words]; Span<sbyte> radialSigns = stackalloc sbyte[24]; |
| | 2824 | 639 | | FiniteSlabPolynomial x = P(radialValues, radialSigns, 0, 5), z = P(radialValues, radialSigns, 5, 5), |
| | 2824 | 640 | | square = P(radialValues, radialSigns, 10, 9), second = P(radialValues, radialSigns, 19, 5); |
| | 2824 | 641 | | Span<ulong> scratchValues = stackalloc ulong[18 * words]; Span<sbyte> scratchSigns = stackalloc sbyte[18]; |
| | 1412 | 642 | | FiniteSlabPolynomial other = new(scratchValues[..(9 * words)], scratchSigns[..9]); |
| | 1412 | 643 | | FiniteSlabPolynomial radiusSquare = new(scratchValues[(9 * words)..], scratchSigns[9..]); |
| | 5018 | 644 | | for (int index = 0; index < problem.Disks.Length; index++) |
| | | 645 | | { |
| | 1412 | 646 | | if ((problem.SolvedDisks & (1 << index)) != 0) continue; |
| | 517 | 647 | | Radial disk = problem.Disks[index]; |
| | 1034 | 648 | | Coordinate(disk, 0, l, m, t, d, x, term); Coordinate(disk, 2, l, m, t, d, z, term); |
| | 517 | 649 | | square.Multiply(x, x); |
| | 1034 | 650 | | other.Multiply(z, z); square.Add(other); |
| | 1551 | 651 | | second.Scale(d, disk.Radius); radiusSquare.Multiply(second, second); square.Add(radiusSquare, -1); |
| | 832 | 652 | | if (FiniteSlabPolynomial.Sign(ref root, square) > 0) return false; |
| | | 653 | | } |
| | 1097 | 654 | | return true; |
| | | 655 | | } |
| | | 656 | | |
| | | 657 | | private static void Coordinate(Radial disk, int axis, FiniteSlabPolynomial l, FiniteSlabPolynomial m, |
| | | 658 | | FiniteSlabPolynomial t, FiniteSlabPolynomial d, FiniteSlabPolynomial result, FiniteSlabPolynomial term) |
| | | 659 | | { |
| | 23973 | 660 | | result.Scale(d, Component(disk.A, axis)); term.Scale(l, Component(disk.E, axis)); result.Add(term); |
| | 31964 | 661 | | term.Scale(m, Component(disk.F, axis)); result.Add(term); term.Scale(t, Component(disk.N, axis)); result.Add(ter |
| | 7991 | 662 | | } |
| | | 663 | | |
| | | 664 | | private static bool IsClearanceExit(scoped Problem problem, scoped ref FiniteAxisValueRoot root, |
| | | 665 | | scoped FiniteSlabPolynomial l, scoped FiniteSlabPolynomial m, |
| | | 666 | | scoped FiniteSlabPolynomial t, scoped FiniteSlabPolynomial d) |
| | | 667 | | { |
| | 717 | 668 | | int words = problem.FieldWords; |
| | 5376 | 669 | | for (int index = problem.ExitLinearStart; index < problem.Lines.Length; index++) |
| | 2318 | 670 | | if (problem.IsActive(index) && problem.Lines[index].C.Sign < 0 |
| | 2665 | 671 | | && (problem.IsPrism || index < problem.ExitLinearStart + 2)) return true; |
| | 740 | 672 | | Span<ulong> values = stackalloc ulong[48 * words]; Span<sbyte> signs = stackalloc sbyte[48]; |
| | 740 | 673 | | FiniteSlabPolynomial x = P(values, signs, 0, 5), y = P(values, signs, 5, 5), |
| | 1110 | 674 | | z = P(values, signs, 10, 5), offset = P(values, signs, 15, 5), term = P(values, signs, 20, 5), |
| | 1110 | 675 | | square = P(values, signs, 25, 9), other = P(values, signs, 34, 9), dot = P(values, signs, 43, 5); |
| | 370 | 676 | | int denominatorSign = FiniteSlabPolynomial.Sign(ref root, d); |
| | 370 | 677 | | Coordinate(problem.Shape, 0, l, m, t, d, x, term); |
| | 370 | 678 | | Coordinate(problem.Shape, 1, l, m, t, d, y, term); |
| | 370 | 679 | | Coordinate(problem.Shape, 2, l, m, t, d, z, term); |
| | 2190 | 680 | | for (int cap = -1; cap <= 1; cap += 2) |
| | | 681 | | { |
| | 735 | 682 | | if (problem.Shape.N.Y.Sign * cap <= 0) continue; |
| | 30 | 683 | | offset.Copy(y); term.Scale(d, Signed(problem.HalfThickness, cap)); offset.Add(term, -1); |
| | 20 | 684 | | if (FiniteSlabPolynomial.Sign(ref root, offset) == 0) return true; |
| | | 685 | | } |
| | 360 | 686 | | if (problem.IsPrism) |
| | | 687 | | { |
| | 856 | 688 | | foreach (ref readonly Linear wall in problem.Lines[problem.ExitLinearStart..]) |
| | | 689 | | { |
| | 380 | 690 | | if (wall.C.Sign >= 0) continue; |
| | 264 | 691 | | offset.Scale(l, wall.A); term.Scale(m, wall.B); offset.Add(term); |
| | 352 | 692 | | term.Scale(t, wall.C); offset.Add(term); term.Scale(d, wall.D); offset.Add(term); |
| | 168 | 693 | | if (FiniteSlabPolynomial.Sign(ref root, offset) == 0) return true; |
| | | 694 | | } |
| | 8 | 695 | | return false; |
| | | 696 | | } |
| | | 697 | | // Axial charts always mark their outward cap active and returned above. |
| | | 698 | | System.Diagnostics.Debug.Assert(!problem.Shape.N.X.IsZero || !problem.Shape.N.Z.IsZero); |
| | 272 | 699 | | if (problem.HasCore) |
| | | 700 | | { |
| | | 701 | | // A straight side only blocks a strictly outward X ray. Sliding |
| | | 702 | | // along it reaches a farther rim and is not an outward-exit sample. |
| | 346 | 703 | | if (problem.Shape.Radius.IsZero && !problem.Shape.N.X.IsZero) return true; |
| | 328 | 704 | | square.Multiply(x, x); term.Scale(d, problem.Shape.Radius); other.Multiply(term, term); square.Add(other, -1 |
| | 82 | 705 | | dot.Scale(x, problem.Shape.N.X); |
| | 82 | 706 | | if (FiniteSlabPolynomial.Sign(ref root, square) == 0 && FiniteSlabPolynomial.Sign(ref root, dot) * denominat |
| | | 707 | | { |
| | 96 | 708 | | offset.Copy(z); term.Scale(d, problem.CoreCap); offset.Add(term, -1); |
| | 32 | 709 | | int upper = FiniteSlabPolynomial.Sign(ref root, offset) * denominatorSign; |
| | 64 | 710 | | offset.Copy(z); offset.Add(term); |
| | 64 | 711 | | if (upper <= 0 && FiniteSlabPolynomial.Sign(ref root, offset) * denominatorSign >= 0) return true; |
| | | 712 | | } |
| | | 713 | | } |
| | 324 | 714 | | for (int side = -1; side <= 1; side += 2) |
| | | 715 | | { |
| | 150 | 716 | | offset.Copy(z); |
| | 326 | 717 | | if (problem.HasCore) { term.Scale(d, Signed(problem.CoreCap, side)); offset.Add(term, -1); } |
| | 150 | 718 | | int hemisphere = FiniteSlabPolynomial.Sign(ref root, offset) * denominatorSign * side; |
| | 150 | 719 | | if (problem.HasCore && hemisphere < 0) continue; |
| | 336 | 720 | | square.Multiply(x, x); other.Multiply(offset, offset); square.Add(other); |
| | 336 | 721 | | term.Scale(d, problem.Shape.Radius); other.Multiply(term, term); square.Add(other, -1); |
| | 112 | 722 | | if (FiniteSlabPolynomial.Sign(ref root, square) != 0) continue; |
| | 312 | 723 | | dot.Scale(x, problem.Shape.N.X); term.Scale(offset, problem.Shape.N.Z); dot.Add(term); |
| | 104 | 724 | | int derivative = FiniteSlabPolynomial.Sign(ref root, dot) * denominatorSign; |
| | 104 | 725 | | if (derivative > 0 || derivative == 0 && (!problem.HasCore || hemisphere > 0 || problem.Shape.N.Z.Sign * sid |
| | 96 | 726 | | return true; |
| | | 727 | | } |
| | 12 | 728 | | return false; |
| | | 729 | | } |
| | | 730 | | |
| | | 731 | | private static void TwoPlanesDisk(scoped Problem problem, Linear a, in Linear b, Radial disk, |
| | | 732 | | ref FiniteSlabRaySelection maximum) |
| | | 733 | | { |
| | 2137 | 734 | | int words = problem.FieldWords; |
| | 4274 | 735 | | Span<ulong> values = stackalloc ulong[22 * words]; Span<sbyte> signs = stackalloc sbyte[22]; |
| | 4274 | 736 | | FiniteSlabPolynomial l = P(values, signs, 0, 2), m = P(values, signs, 2, 2), |
| | 6411 | 737 | | t = P(values, signs, 4, 2), d = P(values, signs, 6), n = P(values, signs, 7), |
| | 6411 | 738 | | equation = P(values, signs, 8, 3), x = P(values, signs, 11, 2), z = P(values, signs, 13, 2), |
| | 8548 | 739 | | term = P(values, signs, 15, 2), square = P(values, signs, 17, 3), depthD = P(values, signs, 20), depthN = P( |
| | 2137 | 740 | | Cross(d, a.A, b.B, a.B, b.A); |
| | 2137 | 741 | | bool depthParameter = !d.IsZero; |
| | 6411 | 742 | | l.Clear(); m.Clear(); t.Clear(); |
| | 2137 | 743 | | if (depthParameter) |
| | | 744 | | { |
| | 1378 | 745 | | Cross(l.Slice(0, 1), a.B, b.D, a.D, b.B); Cross(l.Slice(1, 1), a.B, b.C, a.C, b.B); |
| | 1378 | 746 | | Cross(m.Slice(0, 1), a.D, b.A, a.A, b.D); Cross(m.Slice(1, 1), a.C, b.A, a.A, b.C); |
| | 2067 | 747 | | t.Add(d, 1, 1); depthN.Clear(); depthD.Clear(); |
| | | 748 | | } |
| | | 749 | | else |
| | | 750 | | { |
| | | 751 | | // Dependent lambda/mu rows fix t, leaving a line of source points. |
| | | 752 | | // Its endpoints can be quadratic even though the depth is rational. |
| | 1448 | 753 | | if (a.A.IsZero && a.B.IsZero) |
| | | 754 | | { |
| | 959 | 755 | | if (a.C.IsZero || b.A.IsZero && b.B.IsZero) return; |
| | 1425 | 756 | | n.Set(Negate(a.D)); d.Set(a.C); a = b; |
| | | 757 | | } |
| | | 758 | | else |
| | | 759 | | { |
| | 731 | 760 | | bool first = !a.A.IsZero; |
| | 731 | 761 | | Cross(d, first ? a.A : a.B, b.C, first ? b.A : b.B, a.C); |
| | 1035 | 762 | | if (d.IsZero) return; |
| | 427 | 763 | | Cross(n, a.D, first ? b.A : b.B, b.D, first ? a.A : a.B); |
| | | 764 | | } |
| | | 765 | | // These active planes already fix t=n/d. A negative value is |
| | | 766 | | // infeasible before any radial polynomial or root is constructed. |
| | 957 | 767 | | if (n.Signs[0] * d.Signs[0] < 0) return; |
| | 847 | 768 | | PlaneChart(a, out Signed832 l0, out Signed832 ls, out Signed832 lt, |
| | 847 | 769 | | out Signed832 m0, out Signed832 ms, out Signed832 mt, out Signed832 denominator); |
| | 4235 | 770 | | l.Scale(d, l0); term.Scale(n, lt); l.Add(term); term.Scale(d, ls); l.Add(term, 1, 1); |
| | 4235 | 771 | | m.Scale(d, m0); term.Scale(n, mt); m.Add(term); term.Scale(d, ms); m.Add(term, 1, 1); |
| | 2541 | 772 | | t.Scale(n, denominator); term.Scale(d, denominator); d.Copy(term.Slice(0, 1)); |
| | 3388 | 773 | | depthN.Copy(t.Slice(0, 1)); depthN.Add(d); depthD.Copy(t.Slice(0, 1)); depthD.Add(d, 2); |
| | 1517 | 774 | | if (depthD.Signs[0] < 0) { depthN.Signs[0] *= -1; depthD.Signs[0] *= -1; } |
| | | 775 | | } |
| | 3072 | 776 | | Coordinate(disk, 0, l, m, t, d, x, term); Coordinate(disk, 2, l, m, t, d, z, term); |
| | 4608 | 777 | | equation.Multiply(x, x); square.Multiply(z, z); equation.Add(square); |
| | 4608 | 778 | | term.Scale(d, disk.Radius); square.Multiply(term, term); equation.Add(square, -1); |
| | 1536 | 779 | | ProcessChart(problem, equation, l, m, t, d, depthParameter, depthN, depthD, ref maximum); |
| | 1536 | 780 | | } |
| | | 781 | | |
| | | 782 | | private static void PlaneDisk(scoped Problem problem, in Linear plane, Radial disk, |
| | | 783 | | ref FiniteSlabRaySelection maximum) |
| | | 784 | | { |
| | 1009 | 785 | | int words = problem.FieldWords; |
| | 1009 | 786 | | if (plane.A.IsZero && plane.B.IsZero) |
| | | 787 | | { |
| | 226 | 788 | | if (plane.C.IsZero) return; |
| | 146 | 789 | | if (-plane.D.Sign * plane.C.Sign < 0) return; |
| | | 790 | | // A constant-depth face is a clipped disk/polygon region. Include |
| | | 791 | | // its intrinsic center/cardinal witnesses, even when the complete |
| | | 792 | | // circle lies inside the triangle and no authored edge meets it. |
| | 138 | 793 | | Linear x = new(S(disk.E.X), S(disk.F.X), S(disk.N.X), S(disk.A.X)); |
| | 138 | 794 | | Linear z = new(S(disk.E.Z), S(disk.F.Z), S(disk.N.Z), S(disk.A.Z)); |
| | | 795 | | // The center is not on the disk that selected this family. |
| | 138 | 796 | | Problem centerProblem = problem.WithoutSolvedDisks(); |
| | 138 | 797 | | ThreePlanes(centerProblem, plane, x, z, ref maximum); |
| | 552 | 798 | | for (int index = 0; index < problem.Disks.Length; index++) |
| | | 799 | | { |
| | 138 | 800 | | TwoPlanesDisk(centerProblem.WithDisk(index), plane, x, problem.Disks[index], ref maximum); |
| | 138 | 801 | | TwoPlanesDisk(centerProblem.WithDisk(index), plane, z, problem.Disks[index], ref maximum); |
| | | 802 | | } |
| | 138 | 803 | | return; |
| | | 804 | | } |
| | 1650 | 805 | | Span<ulong> values = stackalloc ulong[25 * words]; Span<sbyte> signs = stackalloc sbyte[25]; |
| | 1650 | 806 | | FiniteSlabPolynomial conic = P(values, signs, 0, 6), equation = P(values, signs, 6, 3), |
| | 1650 | 807 | | sn = P(values, signs, 9, 2), sd = P(values, signs, 11), |
| | 1650 | 808 | | l = P(values, signs, 12, 2), m = P(values, signs, 14, 2), |
| | 1650 | 809 | | t = P(values, signs, 16, 2), d = P(values, signs, 18), |
| | 2475 | 810 | | first = P(values, signs, 19), second = P(values, signs, 20), term = P(values, signs, 21, 2), |
| | 1650 | 811 | | depthN = P(values, signs, 23), depthD = P(values, signs, 24); |
| | 825 | 812 | | PlaneChart(plane, out Signed832 l0, out Signed832 ls, out Signed832 lt, |
| | 825 | 813 | | out Signed832 m0, out Signed832 ms, out Signed832 mt, out Signed832 denominator); |
| | 825 | 814 | | BuildConic(disk, l0, ls, lt, m0, ms, mt, denominator, conic); |
| | 2475 | 815 | | FiniteSlabPolynomial a = conic.Slice(0, 1), b = conic.Slice(1, 1), c = conic.Slice(2, 1), |
| | 2475 | 816 | | e = conic.Slice(4, 1), f = conic.Slice(5, 1), h = conic.Slice(3, 1); |
| | 949 | 817 | | if (a.IsZero) return; // PSD then has no s dependence; a real wall owns the free family. |
| | 3505 | 818 | | sd.Clear(); sd.Add(a, 2); sn.Clear(); sn.Add(h, -1); sn.Add(b, -1, 1); |
| | 3505 | 819 | | l.Scale(sd, l0); term.Scale(sn, ls); l.Add(term); term.Scale(sd, lt); l.Add(term, 1, 1); |
| | 3505 | 820 | | m.Scale(sd, m0); term.Scale(sn, ms); m.Add(term); term.Scale(sd, mt); m.Add(term, 1, 1); |
| | 2103 | 821 | | d.Scale(sd, denominator); t.Clear(); t.Add(d, 1, 1); |
| | 701 | 822 | | equation.Clear(); |
| | 2804 | 823 | | first.Multiply(a, f); second.Multiply(h, h); equation.Add(first, 4); equation.Add(second, -1); |
| | 2804 | 824 | | first.Multiply(a, e); second.Multiply(b, h); equation.Add(first, 4, 1); equation.Add(second, -2, 1); |
| | 2804 | 825 | | first.Multiply(a, c); second.Multiply(b, b); equation.Add(first, 4, 2); equation.Add(second, -1, 2); |
| | 1402 | 826 | | depthN.Clear(); depthD.Clear(); |
| | 701 | 827 | | ProcessChart(problem, equation, l, m, t, d, true, depthN, depthD, ref maximum); |
| | 701 | 828 | | } |
| | | 829 | | |
| | | 830 | | private static void PlaneChart(in Linear plane, out Signed832 l0, out Signed832 ls, out Signed832 lt, |
| | | 831 | | out Signed832 m0, out Signed832 ms, out Signed832 mt, out Signed832 denominator) |
| | | 832 | | { |
| | 1672 | 833 | | if (!plane.B.IsZero) |
| | | 834 | | { |
| | 3393 | 835 | | denominator = plane.B; l0 = lt = default; ls = denominator; |
| | 3393 | 836 | | m0 = Negate(plane.D); ms = Negate(plane.A); mt = Negate(plane.C); |
| | | 837 | | } |
| | | 838 | | else |
| | | 839 | | { |
| | 1623 | 840 | | denominator = plane.A; m0 = mt = default; ms = denominator; |
| | 1623 | 841 | | l0 = Negate(plane.D); ls = Negate(plane.B); lt = Negate(plane.C); |
| | | 842 | | } |
| | 1672 | 843 | | if (denominator.Sign < 0) |
| | | 844 | | { |
| | 2732 | 845 | | denominator = Negate(denominator); l0 = Negate(l0); ls = Negate(ls); lt = Negate(lt); |
| | 2049 | 846 | | m0 = Negate(m0); ms = Negate(ms); mt = Negate(mt); |
| | | 847 | | } |
| | 1672 | 848 | | } |
| | | 849 | | |
| | | 850 | | private static void BuildConic(Radial disk, Signed832 l0, Signed832 ls, Signed832 lt, |
| | | 851 | | Signed832 m0, Signed832 ms, Signed832 mt, Signed832 denominator, FiniteSlabPolynomial conic) |
| | | 852 | | { |
| | 825 | 853 | | int words = conic.FieldWords; |
| | 1650 | 854 | | Span<ulong> values = stackalloc ulong[10 * words]; Span<sbyte> signs = stackalloc sbyte[10]; |
| | 1650 | 855 | | FiniteSlabPolynomial x = P(values, signs, 0, 3), z = P(values, signs, 3, 3), |
| | 3300 | 856 | | first = P(values, signs, 6), second = P(values, signs, 7), factor = P(values, signs, 8), term = P(values, si |
| | 4950 | 857 | | for (int axis = 0; axis < 2; axis++) |
| | | 858 | | { |
| | 3300 | 859 | | FiniteSlabPolynomial coordinate = axis == 0 ? x : z; int component = axis == 0 ? 0 : 2; |
| | 1650 | 860 | | coordinate.Clear(); |
| | 13200 | 861 | | for (int index = 0; index < 3; index++) |
| | | 862 | | { |
| | 9900 | 863 | | first.Set(index == 0 ? l0 : index == 1 ? ls : lt); factor.Set(Component(disk.E, component)); |
| | 9900 | 864 | | term.Multiply(first, factor); coordinate.Slice(index, 1).Copy(term); |
| | 9900 | 865 | | first.Set(index == 0 ? m0 : index == 1 ? ms : mt); factor.Set(Component(disk.F, component)); |
| | 9900 | 866 | | term.Multiply(first, factor); coordinate.Slice(index, 1).Add(term); |
| | 4950 | 867 | | if (index != 1) |
| | | 868 | | { |
| | 6600 | 869 | | first.Set(denominator); factor.Set(index == 0 ? Signed576.ExtendValue(Component(disk.A, component)) |
| | 6600 | 870 | | term.Multiply(first, factor); coordinate.Slice(index, 1).Add(term); |
| | | 871 | | } |
| | | 872 | | } |
| | | 873 | | } |
| | 825 | 874 | | conic.Clear(); |
| | 11550 | 875 | | for (int index = 0; index < 6; index++) |
| | | 876 | | { |
| | 4950 | 877 | | int i = index == 0 || index == 1 || index == 3 ? 1 : index == 2 || index == 4 ? 2 : 0; |
| | 4950 | 878 | | int j = index == 0 ? 1 : index == 1 || index == 2 ? 2 : 0; |
| | 9900 | 879 | | first.Multiply(x.Slice(i, 1), x.Slice(j, 1)); second.Multiply(z.Slice(i, 1), z.Slice(j, 1)); |
| | 9900 | 880 | | conic.Slice(index, 1).Add(first, i == j ? 1 : 2); conic.Slice(index, 1).Add(second, i == j ? 1 : 2); |
| | | 881 | | } |
| | 3300 | 882 | | first.Set(denominator); factor.Set(disk.Radius); term.Multiply(first, factor); first.Multiply(term, term); |
| | 825 | 883 | | conic.Slice(5, 1).Add(first, -1); |
| | 825 | 884 | | } |
| | | 885 | | |
| | | 886 | | private static void Cross(FiniteSlabPolynomial result, Signed832 a, Signed832 b, Signed832 c, Signed832 d) |
| | | 887 | | { |
| | 6051 | 888 | | int words = result.FieldWords; |
| | 12102 | 889 | | Span<ulong> values = stackalloc ulong[5 * words]; Span<sbyte> signs = stackalloc sbyte[5]; |
| | 12102 | 890 | | FiniteSlabPolynomial first = P(values, signs, 0), second = P(values, signs, 1), |
| | 18153 | 891 | | third = P(values, signs, 2), fourth = P(values, signs, 3), term = P(values, signs, 4); |
| | 24204 | 892 | | first.Set(a); second.Set(b); third.Set(c); fourth.Set(d); |
| | 18153 | 893 | | result.Multiply(first, second); term.Multiply(third, fourth); result.Add(term, -1); |
| | 6051 | 894 | | } |
| | | 895 | | |
| | | 896 | | private static void ProcessChart(scoped Problem problem, scoped FiniteSlabPolynomial equation, |
| | | 897 | | scoped FiniteSlabPolynomial l, scoped FiniteSlabPolynomial m, scoped FiniteSlabPolynomial t, |
| | | 898 | | scoped FiniteSlabPolynomial d, bool depthParameter, |
| | | 899 | | scoped FiniteSlabPolynomial depthN, scoped FiniteSlabPolynomial depthD, |
| | | 900 | | ref FiniteSlabRaySelection maximum) |
| | | 901 | | { |
| | 2237 | 902 | | int words = problem.FieldWords; |
| | 4474 | 903 | | Span<ulong> values = stackalloc ulong[11 * words]; Span<sbyte> signs = stackalloc sbyte[11]; |
| | 2237 | 904 | | FiniteSlabPolynomial transformedEquation = P(values, signs, 0, 3), |
| | 4474 | 905 | | tl = P(values, signs, 3, 2), tm = P(values, signs, 5, 2), |
| | 4474 | 906 | | tt = P(values, signs, 7, 2), td = P(values, signs, 9, 2); |
| | 2237 | 907 | | int degree = Math.Max(Math.Max(Degree(l), Degree(m)), Math.Max(Degree(t), Degree(d))); |
| | 10642 | 908 | | for (int chart = 0; chart < (depthParameter ? 1 : 2); chart++) |
| | | 909 | | { |
| | 3084 | 910 | | int sign = chart == 0 ? 1 : -1; |
| | 3084 | 911 | | Transform(equation, transformedEquation, Degree(equation), depthParameter, sign); |
| | 6168 | 912 | | Transform(l, tl, degree, depthParameter, sign); Transform(m, tm, degree, depthParameter, sign); |
| | 6168 | 913 | | Transform(t, tt, degree, depthParameter, sign); Transform(d, td, degree, depthParameter, sign); |
| | | 914 | | // These four homogeneous numerators share one denominator. Strip |
| | | 915 | | // only their common scale, preserving every exact ratio while |
| | | 916 | | // avoiding inherited transform powers in all admission products. |
| | 3084 | 917 | | WideFiniteAxisIntersection.NormalizeFiniteAxisPolynomialPowerOfTwo(values[(3 * words)..], signs[3..]); |
| | 3084 | 918 | | Process(problem, transformedEquation, tl, tm, tt, td, depthParameter, depthN, depthD, ref maximum); |
| | | 919 | | } |
| | 2237 | 920 | | if (!depthParameter && equation.Signs[0] == 0) |
| | | 921 | | { |
| | 210 | 922 | | transformedEquation.Clear(); transformedEquation.Set(-1); transformedEquation.Set(2, 1); |
| | 70 | 923 | | Process(problem, transformedEquation, l.Slice(0, 1), m.Slice(0, 1), t.Slice(0, 1), d.Slice(0, 1), |
| | 70 | 924 | | false, depthN, depthD, ref maximum); |
| | | 925 | | } |
| | 2237 | 926 | | } |
| | | 927 | | |
| | | 928 | | private static void Transform(FiniteSlabPolynomial source, FiniteSlabPolynomial result, |
| | | 929 | | int degree, bool depthParameter, int parameterSign) |
| | | 930 | | { |
| | | 931 | | // t=(2u-1)/(1-u) keeps nonnegative depth in [1/2,1). |
| | | 932 | | // A free line uses s=+/-u/(1-u); its zero endpoint is separate. |
| | 15420 | 933 | | result.Clear(); |
| | 80316 | 934 | | for (int index = 0; index <= Degree(source); index++) |
| | 111048 | 935 | | for (int first = 0; first <= (depthParameter ? index : 0); first++) |
| | 164976 | 936 | | for (int second = 0; second <= degree - index; second++) |
| | | 937 | | { |
| | 51702 | 938 | | int exponent = depthParameter ? first + second : index + second; |
| | 51702 | 939 | | int factor = Choose(degree - index, second) * (second % 2 == 0 ? 1 : -1); |
| | 80322 | 940 | | if (depthParameter) factor *= Choose(index, first) * (1 << first) * ((index - first) % 2 == 0 ? 1 : |
| | 25386 | 941 | | else if (parameterSign < 0 && index % 2 != 0) factor = -factor; |
| | 51702 | 942 | | result.Add(source.Slice(index, 1), factor, exponent); |
| | | 943 | | } |
| | 15420 | 944 | | } |
| | | 945 | | private static int Choose(int n, int k) |
| | 417450 | 946 | | { int value = 1; for (int i = 1; i <= k; i++) value = value * (n - i + 1) / i; return value; } |
| | | 947 | | private static int Degree(FiniteSlabPolynomial polynomial) |
| | 190556 | 948 | | { int degree = polynomial.Count - 1; while (degree > 0 && polynomial.Signs[degree] == 0) degree--; return degree; } |
| | | 949 | | |
| | | 950 | | private static FiniteSlabPolynomial P(Span<ulong> values, Span<sbyte> signs, int index, int count = 1) => |
| | 212417 | 951 | | FiniteSlabPolynomial.At(values, signs, index, count); |
| | 2372 | 952 | | private static Signed576 S(int value) => S(Signed320.ExtendValue(Signed192.Signed(value))); |
| | 13910 | 953 | | private static Signed576 S(Signed320 value) => Signed576.ExtendValue(value); |
| | 3532 | 954 | | private static Signed576 S(Signed576 value) => value; |
| | | 955 | | private static Signed320 ScaledOffset(Fixed64 value, Signed320 scale) => |
| | 2020 | 956 | | Signed320.NarrowValue(WideArithmetic.MultiplySigned320(scale, Signed192.Raw(value))); |
| | 4358 | 957 | | private static Signed320 Negate(Signed320 value) => WideArithmetic.Negate(value); |
| | 3329 | 958 | | private static Signed320 Signed(Signed320 value, int sign) => sign < 0 ? Negate(value) : value; |
| | 10326 | 959 | | private static Signed832 Negate(Signed832 value) => WideArithmetic.SubtractSigned832(default, value); |
| | 5904 | 960 | | private static Signed576 Signed(Signed576 value, int sign) => sign < 0 ? WideArithmetic.SubtractSigned576(default, v |
| | 558 | 961 | | private static WideDirection3d Signed(WideDirection3d p, int sign) => new(Signed(p.X, sign), Signed(p.Y, sign), Sign |
| | 9641 | 962 | | private static Signed576 Component(WideDirection3d p, int axis) => axis == 0 ? p.X : axis == 1 ? p.Y : p.Z; |
| | 35523 | 963 | | private static Signed320 Component(WideAxis3 p, int axis) => axis == 0 ? p.X : axis == 1 ? p.Y : p.Z; |
| | 1082 | 964 | | private static WideAxis3 Difference(WideAxis3 p, WideAxis3 q) => new( |
| | 1082 | 965 | | WideArithmetic.SubtractSigned320(p.X, q.X), WideArithmetic.SubtractSigned320(p.Y, q.Y), |
| | 1082 | 966 | | WideArithmetic.SubtractSigned320(p.Z, q.Z)); |
| | 2160 | 967 | | private static Signed576 Wall(WideAxis3 p, Signed320 x, Signed320 z) => WideArithmetic.SubtractSigned576( |
| | 2160 | 968 | | WideArithmetic.MultiplySigned320(x, p.Z), WideArithmetic.MultiplySigned320(z, p.X)); |
| | | 969 | | private static Signed576 Wall(WideDirection3d p, Signed320 x, Signed320 z) => |
| | 720 | 970 | | Signed576.NarrowValue(WideArithmetic.SubtractSigned704( |
| | 720 | 971 | | WideArithmetic.MultiplySigned576ToSigned704(p.Z, x), |
| | 720 | 972 | | WideArithmetic.MultiplySigned576ToSigned704(p.X, z))); |
| | 59994 | 973 | | private static int DomainEdge(int line, int count) => count == 0 ? 0 : line == 0 ? 4 : line == 1 ? 1 : line == 2 ? 2 |
| | | 974 | | } |