| | | 1 | | //======================================================================= |
| | | 2 | | // ConePlaneRayPointExits.cs |
| | | 3 | | //======================================================================= |
| | | 4 | | // MIT License, Copyright (c) 2024-present David Oravsky (mrdav30) |
| | | 5 | | // See LICENSE file in the project root for full license information. |
| | | 6 | | //======================================================================= |
| | | 7 | | |
| | | 8 | | using System; |
| | | 9 | | |
| | | 10 | | namespace FixedMathSharp.Geometry; |
| | | 11 | | |
| | | 12 | | /// <summary>Exact closed admission and first finite exits of retained plane-section points.</summary> |
| | | 13 | | internal static class ConePlaneRayPointExits |
| | | 14 | | { |
| | | 15 | | // A point certificate stores four homogeneous quadratic-field integers; |
| | | 16 | | // their rational/radical coefficients are <3072 bits, root <2498 bits. |
| | | 17 | | // One cone product needs <2*3072+2498+126+3=8771 bits (144 words). |
| | | 18 | | // Dot/axial/wall predicates need <3072+526+3=3601 bits (64 words). |
| | | 19 | | // Depth fields are <4096 bits. Stationary exit verification multiplies |
| | | 20 | | // one depth, one point, A<650 and the shared root; <10320 bits, 176 words. |
| | | 21 | | // These ordinary-path products are transient, never retained as point/depth |
| | | 22 | | // fields. Exact affine/wide-clip cohorts select the frame's wider banks. |
| | | 23 | | |
| | | 24 | | /// <summary>Tests closed shared-plane and finite-source admission in the cone.</summary> |
| | | 25 | | /// <remarks>A side certificate proves F(point)=0 exactly in the supplied root; only that solid predicate may then b |
| | | 26 | | internal static bool ContainsPoint(in ConePlaneRayEventSource source, in ConePlaneRayFrame frame, |
| | | 27 | | scoped ConePlaneRayPoint point, scoped ReadOnlySpan<ulong> root, bool sideCertificate = false) |
| | | 28 | | { |
| | 12450 | 29 | | if (frame.NormalSquared.IsZero || point.Denominator.Sign(root) <= 0) |
| | 2 | 30 | | return false; |
| | 12448 | 31 | | Span<ulong> work = stackalloc ulong[4 * frame.FieldWords]; |
| | 12448 | 32 | | Span<int> signs = stackalloc int[4]; |
| | 12448 | 33 | | ContactQuadratic value = ContactQuadratic.At(work, signs, 0, frame.FieldWords); |
| | 12448 | 34 | | ContactQuadratic term = ContactQuadratic.At(work, signs, 1, frame.FieldWords); |
| | 12448 | 35 | | Dot(point, frame.Normal.X, frame.Normal.Y, |
| | 12448 | 36 | | frame.Normal.Z, value, term); |
| | 24896 | 37 | | ContactQuadratic.Scale(point.Denominator, frame.PlaneConstant, term); value.Add(term, -1); |
| | 12448 | 38 | | if (value.Sign(root) != 0 || point.Y.Sign(root) < 0) |
| | 1131 | 39 | | return false; |
| | 11317 | 40 | | ContactQuadratic.Scale(point.Denominator, Signed576.ExtendValue(frame.Finite.FullHeight), value); |
| | 11317 | 41 | | value.Add(point.Y, -1); |
| | 11317 | 42 | | if (value.Sign(root) < 0) |
| | 1591 | 43 | | return false; |
| | 9726 | 44 | | if (!sideCertificate) |
| | | 45 | | { |
| | 12604 | 46 | | ContactQuadratic x = point.X, z = point.Z; |
| | | 47 | | // The admitted axial interval already proves an axis point is in |
| | | 48 | | // the cone: F(0,Y,0)=-R^2Y^2<=0, also when R=0. Use identically |
| | | 49 | | // zero fields only; root-dependent cancellation keeps the full test. |
| | 6302 | 50 | | if ((x.Signs[0] | x.Signs[1] | z.Signs[0] | z.Signs[1]) != 0) |
| | | 51 | | { |
| | 5563 | 52 | | Span<ulong> coneWork = stackalloc ulong[2 * frame.AdmissionWords]; |
| | 5563 | 53 | | Span<int> coneSigns = stackalloc int[2]; |
| | 5563 | 54 | | ContactQuadratic cone = ContactQuadratic.At(coneWork, coneSigns, 0, frame.AdmissionWords); |
| | 5563 | 55 | | ConeProduct(point, point, frame, root, cone); |
| | 9642 | 56 | | if (cone.Sign(root) > 0) return false; |
| | | 57 | | } |
| | | 58 | | } |
| | | 59 | | // Native segment bounds are checked by the sink using its same exact |
| | | 60 | | // borrowed endpoints, before a candidate can enter the retained maximum. |
| | 10875 | 61 | | if (source.Scope != ConePlaneRayEventScope.Segment) return true; |
| | 419 | 62 | | frame.GetEdgeLine(new FixedSegment(source.A, source.B), out Signed576 wx, out Signed576 wy, out Signed576 wz, ou |
| | 419 | 63 | | Dot(point, wx, wy, wz, value, term); |
| | 838 | 64 | | ContactQuadratic.Scale(point.Denominator, offset, term); value.Add(term, -1); |
| | 420 | 65 | | if (value.Sign(root) != 0) return false; |
| | 836 | 66 | | WideAxis3 start = frame.Transform(source.A), end = frame.Transform(source.B); |
| | 418 | 67 | | WideAxis3 delta = new(WideArithmetic.SubtractSigned320(end.X, start.X), |
| | 418 | 68 | | WideArithmetic.SubtractSigned320(end.Y, start.Y), WideArithmetic.SubtractSigned320(end.Z, start.Z)); |
| | 418 | 69 | | if (delta.IsZero) |
| | | 70 | | { |
| | | 71 | | // A point segment has no line or projection direction; its |
| | | 72 | | // finite domain is exact homogeneous coordinate equality. |
| | 352 | 73 | | for (int axis = 0; axis < 3; axis++) |
| | | 74 | | { |
| | 133 | 75 | | (axis == 0 ? point.X : axis == 1 ? point.Y : point.Z).CopyTo(value); |
| | 133 | 76 | | ContactQuadratic.Scale(point.Denominator, Signed576.ExtendValue(axis == 0 ? start.X : axis == 1 ? start. |
| | 133 | 77 | | value.Add(term, -1); |
| | 135 | 78 | | if (value.Sign(root) != 0) return false; |
| | | 79 | | } |
| | 43 | 80 | | return true; |
| | | 81 | | } |
| | | 82 | | // Point coefficients <3072 and transformed edge <198 keep the |
| | | 83 | | // complete projection and squared-edge bounds <3472 bits (64 words). |
| | 373 | 84 | | Dot(point, Signed576.ExtendValue(delta.X), Signed576.ExtendValue(delta.Y), Signed576.ExtendValue(delta.Z), value |
| | 746 | 85 | | ContactQuadratic.Scale(point.Denominator, WideDirection3d.Dot(new WideDirection3d(start), delta), term); value.A |
| | 376 | 86 | | if (value.Sign(root) < 0) return false; |
| | 740 | 87 | | ContactQuadratic.Scale(point.Denominator, new WideDirection3d(delta).SquaredLength, term); term.Add(value, -1); |
| | 370 | 88 | | return term.Sign(root) >= 0; |
| | | 89 | | } |
| | | 90 | | |
| | | 91 | | /// <summary>Tests only the closed triangle walls of an already certified finite-cone point in the same plane.</summ |
| | | 92 | | internal static bool ContainsTrianglePoint(FixedTriangle triangle, in ConePlaneRayFrame frame, |
| | | 93 | | scoped ConePlaneRayPoint point, scoped ReadOnlySpan<ulong> root) |
| | | 94 | | { |
| | 4076 | 95 | | triangle.GetExactNormal(out Signed192 nx, out Signed192 ny, out Signed192 nz, out _); |
| | 4094 | 96 | | if (nx.IsZero && ny.IsZero && nz.IsZero) return false; |
| | 8116 | 97 | | Span<ulong> work = stackalloc ulong[4 * frame.FieldWords]; Span<int> signs = stackalloc int[4]; |
| | 8116 | 98 | | ContactQuadratic value = ContactQuadratic.At(work, signs, 0, frame.FieldWords), term = ContactQuadratic.At(work, |
| | 31930 | 99 | | for (int edge = 0; edge < 3; edge++) |
| | | 100 | | { |
| | 12032 | 101 | | frame.GetTriangleEdgeWall(triangle, edge, out Signed576 x, out Signed576 y, out Signed576 z, out Signed832 o |
| | 12032 | 102 | | Dot(point, x, y, z, value, term); |
| | 24064 | 103 | | ContactQuadratic.Scale(point.Denominator, offset, term); value.Add(term, -1); |
| | 12157 | 104 | | if (value.Sign(root) < 0) return false; |
| | | 105 | | } |
| | 3933 | 106 | | return true; |
| | | 107 | | } |
| | | 108 | | |
| | | 109 | | /// <summary>Tests a certified source-frame point against triangle walls in another frame of the same cone.</summary |
| | | 110 | | /// <remarks>The target triangle and the certified point are coplanar; both frames describe the same physical cone.< |
| | | 111 | | internal static bool ContainsTrianglePoint(FixedTriangle triangle, in ConePlaneRayFrame frame, |
| | | 112 | | scoped ConePlaneRayPoint point, scoped ReadOnlySpan<ulong> root, in ConePlaneRayFrame sourceFrame) |
| | | 113 | | { |
| | 50 | 114 | | Signed320 sourceScale = sourceFrame.Finite.ShapeFrame.CoordinateDenominator; |
| | 50 | 115 | | Signed320 targetScale = frame.Finite.ShapeFrame.CoordinateDenominator; |
| | 75 | 116 | | if (sourceScale.Equals(targetScale)) return ContainsTrianglePoint(triangle, frame, point, root); |
| | | 117 | | // Apex-to-base coordinates share the same physical cone frame but |
| | | 118 | | // carry 2D native scale. Cross multiplication preserves every root |
| | | 119 | | // and closed wall sign without materializing a rounded point. Point |
| | | 120 | | // fields <3072 plus D<132 and wall offsets <526 stay <3733 bits. |
| | 25 | 121 | | Span<ulong> values = stackalloc ulong[frame.PointStorageWords]; |
| | 25 | 122 | | Span<int> signs = stackalloc int[ConePlaneRayPoint.SignCount]; |
| | 25 | 123 | | ConePlaneRayPoint transformed = new(values, signs, frame.FieldWords); |
| | 25 | 124 | | point.CopyTo(transformed, sourceScale, targetScale); |
| | 25 | 125 | | return ContainsTrianglePoint(triangle, frame, transformed, root); |
| | | 126 | | } |
| | | 127 | | |
| | | 128 | | internal static bool AccumulateRationalPoint(in ConePlaneRayEventSource source, in ConePlaneRayFrame frame, |
| | | 129 | | scoped ConePlaneRayPoint point, scoped ref ConePlaneRaySelection positive, scoped ref ConePlaneRaySelection nega |
| | | 130 | | bool pointOnly = false, int requestedOrientation = 0) |
| | | 131 | | { |
| | 5888 | 132 | | Span<ulong> root = stackalloc ulong[frame.RootWords]; root.Clear(); |
| | 2944 | 133 | | if (!ContainsPoint(source, frame, point, root)) |
| | 2295 | 134 | | return false; |
| | 698 | 135 | | if (pointOnly) return true; |
| | | 136 | | // Only vertices, rational edge endpoints and the axis-plane point use |
| | | 137 | | // this path: coordinate numerators <526, denominator <326. Then the |
| | | 138 | | // ray equation A<1302, B<1241, C<1181 has discriminant <2485, |
| | | 139 | | // and its root fits the 40-word chart. |
| | | 140 | | // Rational cuts with positive <319-bit denominators give coordinates |
| | | 141 | | // <198+319+1=518 bits and remain inside the same endpoint bounds. |
| | 1172 | 142 | | if (requestedOrientation >= 0) AccumulateRationalOrientation(point, frame, 1, ref positive); |
| | 1172 | 143 | | if (requestedOrientation <= 0) AccumulateRationalOrientation(point, frame, -1, ref negative); |
| | 600 | 144 | | return true; |
| | | 145 | | } |
| | | 146 | | |
| | | 147 | | internal static bool AccumulateSidePoint(in ConePlaneRayEventSource source, in ConePlaneRayFrame frame, |
| | | 148 | | scoped ConePlaneRayPoint point, scoped ReadOnlySpan<ulong> root, |
| | | 149 | | scoped ref ConePlaneRaySelection positive, scoped ref ConePlaneRaySelection negative, |
| | | 150 | | bool pointOnly = false, int requestedOrientation = 0) |
| | | 151 | | { |
| | 3985 | 152 | | if (!ContainsPoint(source, frame, point, root, sideCertificate: true)) |
| | 1007 | 153 | | return false; |
| | 3027 | 154 | | if (pointOnly) return true; |
| | | 155 | | // The chart certifies F(point)=0. The nonzero crossing -2B/A is |
| | | 156 | | // quadratic-field linear; solving a second quadratic would introduce |
| | | 157 | | // a needless nested radical. Inward zero crossings are rejected. |
| | 5770 | 158 | | if (requestedOrientation >= 0) AccumulateSideOrientation(point, frame, root, 1, ref positive); |
| | 5770 | 159 | | if (requestedOrientation <= 0) AccumulateSideOrientation(point, frame, root, -1, ref negative); |
| | 2929 | 160 | | return true; |
| | | 161 | | } |
| | | 162 | | |
| | | 163 | | internal static bool AccumulateCertifiedExit(in ConePlaneRayEventSource source, in ConePlaneRayFrame frame, |
| | | 164 | | scoped ConePlaneRayPoint point, scoped ReadOnlySpan<ulong> root, scoped ContactQuadratic numerator, |
| | | 165 | | scoped ContactQuadratic denominator, int orientation, |
| | | 166 | | scoped ref ConePlaneRaySelection positive, scoped ref ConePlaneRaySelection negative, |
| | | 167 | | bool pointOnly = false, bool sideCertificate = false) |
| | | 168 | | { |
| | 3160 | 169 | | if (!ContainsPoint(source, frame, point, root, sideCertificate) || numerator.Sign(root) < 0 || denominator.Sign( |
| | 1791 | 170 | | return false; |
| | | 171 | | // The constructing chart certifies an upper support/base/rim point |
| | | 172 | | // and its first-exit branch. A side certificate skips only F(point); |
| | | 173 | | // closed plane/height/source admission and exit signs remain checked. |
| | 1393 | 174 | | if (pointOnly) return true; |
| | 2000 | 175 | | if (orientation > 0) positive.Keep(numerator, denominator, root, frame, point); |
| | 690 | 176 | | else negative.Keep(numerator, denominator, root, frame, point); |
| | 1345 | 177 | | return true; |
| | | 178 | | } |
| | | 179 | | |
| | | 180 | | internal static bool AccumulateStationarySideExit(in ConePlaneRayEventSource source, in ConePlaneRayFrame frame, |
| | | 181 | | scoped ConePlaneRayPoint point, scoped ReadOnlySpan<ulong> root, scoped ContactQuadratic numerator, |
| | | 182 | | scoped ContactQuadratic denominator, int orientation, |
| | | 183 | | scoped ref ConePlaneRaySelection positive, scoped ref ConePlaneRaySelection negative, bool pointOnly = false) |
| | | 184 | | { |
| | 2171 | 185 | | if (!ContainsPoint(source, frame, point, root) || numerator.Sign(root) < 0 || denominator.Sign(root) <= 0) |
| | 1687 | 186 | | return false; |
| | 484 | 187 | | Span<ulong> work = stackalloc ulong[10 * frame.VerificationWords]; |
| | 484 | 188 | | Span<int> signs = stackalloc int[10]; |
| | 484 | 189 | | ContactQuadratic y = ContactQuadratic.At(work, signs, 0, frame.VerificationWords); |
| | 484 | 190 | | ContactQuadratic d = ContactQuadratic.At(work, signs, 1, frame.VerificationWords); |
| | 484 | 191 | | ContactQuadratic term = ContactQuadratic.At(work, signs, 2, frame.VerificationWords); |
| | 484 | 192 | | ContactQuadratic derivative = ContactQuadratic.At(work, signs, 3, frame.VerificationWords); |
| | 484 | 193 | | ContactQuadratic a = ContactQuadratic.At(work, signs, 4, frame.VerificationWords); |
| | 484 | 194 | | ContactQuadratic.Multiply(point.Y, denominator, root, y); |
| | 484 | 195 | | ContactQuadratic.Multiply(point.Denominator, numerator, root, term); |
| | 1452 | 196 | | ContactQuadratic.Scale(term, frame.Normal.Y, d); d.MultiplySign(orientation); y.Add(d); |
| | 484 | 197 | | ContactQuadratic.Multiply(point.Denominator, denominator, root, d); |
| | 484 | 198 | | if (y.Sign(root) < 0) |
| | 11 | 199 | | return false; |
| | 946 | 200 | | ContactQuadratic.Scale(d, Signed576.ExtendValue(frame.Finite.FullHeight), term); term.Add(y, -1); |
| | 473 | 201 | | if (term.Sign(root) < 0) |
| | 166 | 202 | | return false; |
| | 307 | 203 | | GetRayCoefficients(point, frame, root, a, derivative); |
| | 307 | 204 | | derivative.MultiplySign(orientation); |
| | 307 | 205 | | ContactQuadratic.Multiply(derivative, denominator, root, y); |
| | 307 | 206 | | ContactQuadratic.Multiply(a, numerator, root, d); |
| | 614 | 207 | | ContactQuadratic.Multiply(d, point.Denominator, root, term); y.Add(term); |
| | | 208 | | // F(q)=0 is certified by the stationary elimination polynomial. A |
| | | 209 | | // nonnegative derivative selects the exit rather than opposite nappe. |
| | 307 | 210 | | if (y.Sign(root) < 0 || y.Sign(root) == 0 && a.Sign(root) < 0) |
| | 20 | 211 | | return false; |
| | 295 | 212 | | if (pointOnly) return true; |
| | 441 | 213 | | if (orientation > 0) positive.Keep(numerator, denominator, root, frame, point); |
| | 117 | 214 | | else negative.Keep(numerator, denominator, root, frame, point); |
| | 279 | 215 | | return true; |
| | | 216 | | } |
| | | 217 | | |
| | | 218 | | private static void Dot(scoped ConePlaneRayPoint point, Signed576 x, Signed576 y, Signed576 z, |
| | | 219 | | ContactQuadratic result, ContactQuadratic term) |
| | | 220 | | { |
| | 25272 | 221 | | ContactQuadratic.Scale(point.X, x, result); |
| | 50544 | 222 | | ContactQuadratic.Scale(point.Y, y, term); result.Add(term); |
| | 50544 | 223 | | ContactQuadratic.Scale(point.Z, z, term); result.Add(term); |
| | 25272 | 224 | | } |
| | | 225 | | |
| | | 226 | | private static void ConeProduct(ConePlaneRayPoint left, ConePlaneRayPoint right, |
| | | 227 | | in ConePlaneRayFrame frame, scoped ReadOnlySpan<ulong> root, ContactQuadratic result) |
| | | 228 | | { |
| | 19846 | 229 | | int words = result.FieldWords; |
| | 39692 | 230 | | Span<ulong> work = stackalloc ulong[4 * words]; Span<int> signs = stackalloc int[4]; |
| | 19846 | 231 | | ContactQuadratic term = ContactQuadratic.At(work, signs, 0, words); |
| | 19846 | 232 | | ContactQuadratic scaled = ContactQuadratic.At(work, signs, 1, words); |
| | 19846 | 233 | | Signed576 hh = Signed576.ExtendValue(WideArithmetic.MultiplySigned192(Signed192.Raw(frame.Height), Signed192.Raw |
| | 19846 | 234 | | Signed576 rr = Signed576.ExtendValue(WideArithmetic.MultiplySigned192(Signed192.Raw(frame.Radius), Signed192.Raw |
| | 19846 | 235 | | ContactQuadratic.Multiply(left.X, right.X, root, result); |
| | 39692 | 236 | | ContactQuadratic.Multiply(left.Z, right.Z, root, term); result.Add(term); |
| | 39692 | 237 | | ContactQuadratic.Scale(result, hh, scaled); scaled.CopyTo(result); |
| | 19846 | 238 | | ContactQuadratic.Multiply(left.Y, right.Y, root, term); |
| | 39692 | 239 | | ContactQuadratic.Scale(term, rr, scaled); result.Add(scaled, -1); |
| | 19846 | 240 | | } |
| | | 241 | | |
| | | 242 | | private static void GetRayCoefficients(scoped ConePlaneRayPoint point, in ConePlaneRayFrame frame, |
| | | 243 | | scoped ReadOnlySpan<ulong> root, ContactQuadratic a, ContactQuadratic b) |
| | | 244 | | { |
| | 6704 | 245 | | Span<ulong> normalWords = stackalloc ulong[frame.PointStorageWords]; |
| | 6704 | 246 | | Span<int> normalSigns = stackalloc int[ConePlaneRayPoint.SignCount]; |
| | 13408 | 247 | | ConePlaneRayPoint normal = new(normalWords, normalSigns, frame.FieldWords); normal.Set(frame.Normal); |
| | 6704 | 248 | | ConeProduct(normal, normal, frame, root, a); |
| | 6704 | 249 | | ConeProduct(point, normal, frame, root, b); |
| | 6704 | 250 | | } |
| | | 251 | | |
| | | 252 | | private static void AccumulateRationalOrientation(scoped ConePlaneRayPoint point, in ConePlaneRayFrame frame, |
| | | 253 | | int orientation, scoped ref ConePlaneRaySelection selection) |
| | | 254 | | { |
| | 1144 | 255 | | if (TryAccumulateAxialOrientation(point, frame, ReadOnlySpan<ulong>.Empty, |
| | 1413 | 256 | | orientation, sidePoint: false, ref selection)) return; |
| | 1750 | 257 | | Span<ulong> work = stackalloc ulong[20 * frame.FieldWords]; Span<int> signs = stackalloc int[20]; |
| | 875 | 258 | | ContactQuadratic a = ContactQuadratic.At(work, signs, 0, frame.FieldWords); |
| | 875 | 259 | | ContactQuadratic b = ContactQuadratic.At(work, signs, 1, frame.FieldWords); |
| | 875 | 260 | | ContactQuadratic c = ContactQuadratic.At(work, signs, 2, frame.FieldWords); |
| | 875 | 261 | | ContactQuadratic term = ContactQuadratic.At(work, signs, 3, frame.FieldWords); |
| | 875 | 262 | | ContactQuadratic polynomialA = ContactQuadratic.At(work, signs, 4, frame.FieldWords); |
| | 875 | 263 | | ContactQuadratic polynomialB = ContactQuadratic.At(work, signs, 5, frame.FieldWords); |
| | 875 | 264 | | ContactQuadratic numerator = ContactQuadratic.At(work, signs, 6, frame.FieldWords); |
| | 875 | 265 | | ContactQuadratic denominator = ContactQuadratic.At(work, signs, 7, frame.FieldWords); |
| | 875 | 266 | | ContactQuadratic bestN = ContactQuadratic.At(work, signs, 8, frame.FieldWords); |
| | 875 | 267 | | ContactQuadratic bestD = ContactQuadratic.At(work, signs, 9, frame.FieldWords); |
| | 1750 | 268 | | Span<ulong> root = stackalloc ulong[frame.RootWords]; root.Clear(); |
| | 1750 | 269 | | Span<ulong> bestRoot = stackalloc ulong[frame.RootWords]; bestRoot.Clear(); |
| | 875 | 270 | | bool hasValue = false; |
| | 1750 | 271 | | GetRayCoefficients(point, frame, root, a, b); b.MultiplySign(orientation); |
| | 875 | 272 | | ConeProduct(point, point, frame, root, c); |
| | 875 | 273 | | ContactQuadratic.Multiply(point.Denominator, point.Denominator, root, term); |
| | 875 | 274 | | ContactQuadratic.Multiply(a, term, root, polynomialA); |
| | 875 | 275 | | ContactQuadratic.Multiply(b, point.Denominator, root, polynomialB); |
| | 5250 | 276 | | for (int branch = -1; branch <= 1; branch += 2) |
| | | 277 | | { |
| | 1750 | 278 | | if (!ContactQuadratic.TryGetParameter(polynomialA, polynomialB, c, branch, root, numerator, denominator) |
| | 1750 | 279 | | || numerator.Sign(root) < 0) |
| | | 280 | | continue; |
| | | 281 | | // Derivative b*d+A*n at a root picks the first outward crossing. |
| | | 282 | | // A concave repeated root touches F=0 from inside and is not exit. |
| | 907 | 283 | | ContactQuadratic.Multiply(polynomialB, denominator, root, term); |
| | 1814 | 284 | | ContactQuadratic.Multiply(polynomialA, numerator, root, b); term.Add(b); |
| | 907 | 285 | | int derivative = term.Sign(root); |
| | 907 | 286 | | if (derivative < 0 || derivative == 0 && polynomialA.Signs[0] < 0) |
| | | 287 | | continue; |
| | 804 | 288 | | KeepFirstExit(numerator, denominator, root, bestN, bestD, bestRoot, ref hasValue); |
| | | 289 | | } |
| | 875 | 290 | | root.Clear(); |
| | 875 | 291 | | if (GetAxialExit(point, frame, orientation, numerator, denominator)) |
| | 681 | 292 | | KeepFirstExit(numerator, denominator, root, bestN, bestD, bestRoot, ref hasValue); |
| | | 293 | | // A nonzero ray from an admitted point in the compact finite cone |
| | | 294 | | // must meet a side or axial boundary, including zero-length exits. |
| | | 295 | | System.Diagnostics.Debug.Assert(hasValue); |
| | 875 | 296 | | selection.Keep(bestN, bestD, bestRoot, frame, point); |
| | 875 | 297 | | } |
| | | 298 | | |
| | | 299 | | private static void AccumulateSideOrientation(scoped ConePlaneRayPoint point, in ConePlaneRayFrame frame, |
| | | 300 | | scoped ReadOnlySpan<ulong> root, int orientation, scoped ref ConePlaneRaySelection selection) |
| | | 301 | | { |
| | 5682 | 302 | | if (TryAccumulateAxialOrientation(point, frame, root, |
| | 5842 | 303 | | orientation, sidePoint: true, ref selection)) return; |
| | 11044 | 304 | | Span<ulong> work = stackalloc ulong[12 * frame.FieldWords]; Span<int> signs = stackalloc int[12]; |
| | 5522 | 305 | | ContactQuadratic a = ContactQuadratic.At(work, signs, 0, frame.FieldWords); |
| | 5522 | 306 | | ContactQuadratic b = ContactQuadratic.At(work, signs, 1, frame.FieldWords); |
| | 5522 | 307 | | ContactQuadratic numerator = ContactQuadratic.At(work, signs, 2, frame.FieldWords); |
| | 5522 | 308 | | ContactQuadratic denominator = ContactQuadratic.At(work, signs, 3, frame.FieldWords); |
| | 5522 | 309 | | ContactQuadratic bestN = ContactQuadratic.At(work, signs, 4, frame.FieldWords); |
| | 5522 | 310 | | ContactQuadratic bestD = ContactQuadratic.At(work, signs, 5, frame.FieldWords); |
| | 11044 | 311 | | Span<ulong> bestRoot = stackalloc ulong[frame.RootWords]; bestRoot.Clear(); |
| | 5522 | 312 | | bool hasValue = false; |
| | 11044 | 313 | | GetRayCoefficients(point, frame, root, a, b); b.MultiplySign(orientation); |
| | 11044 | 314 | | int bSign = b.Sign(root), aSign = a.Sign(root); |
| | 5522 | 315 | | if (bSign > 0 || bSign == 0 && aSign > 0) |
| | | 316 | | { |
| | 6698 | 317 | | numerator.Clear(); denominator.Set(Signed576.ExtendValue(Signed320.ExtendValue(Signed192.Signed(1)))); |
| | 3349 | 318 | | KeepFirstExit(numerator, denominator, root, bestN, bestD, bestRoot, ref hasValue); |
| | | 319 | | } |
| | 2173 | 320 | | else if (bSign < 0 && aSign > 0) |
| | | 321 | | { |
| | 3688 | 322 | | b.CopyTo(numerator, -1); numerator.Add(numerator); |
| | 1844 | 323 | | ContactQuadratic.Multiply(a, point.Denominator, root, denominator); |
| | 1844 | 324 | | KeepFirstExit(numerator, denominator, root, bestN, bestD, bestRoot, ref hasValue); |
| | | 325 | | } |
| | | 326 | | // A<=0 with B<=0 remains in this nappe until the finite axial cap. |
| | | 327 | | // A=B=0 is a generator interval and likewise keeps its finite endpoint. |
| | 5522 | 328 | | if (GetAxialExit(point, frame, orientation, numerator, denominator)) |
| | 4286 | 329 | | KeepFirstExit(numerator, denominator, root, bestN, bestD, bestRoot, ref hasValue); |
| | | 330 | | // With Ny=0, A=HÂ^2(NxÂ^2+NzÂ^2)>0 and a side exit exists. Otherwise |
| | | 331 | | // the admitted finite axial interval supplies an endpoint. |
| | | 332 | | System.Diagnostics.Debug.Assert(hasValue); |
| | 5522 | 333 | | selection.Keep(bestN, bestD, bestRoot, frame, point); |
| | 5522 | 334 | | } |
| | | 335 | | |
| | | 336 | | /// <summary>Applies the axial first-exit construction to an already admitted point.</summary> |
| | | 337 | | /// <remarks>The orientation is +1 or -1. The caller proves exact finite admission; sidePoint also certifies the zer |
| | | 338 | | internal static bool TryAccumulateAxialOrientation(scoped ConePlaneRayPoint point, in ConePlaneRayFrame frame, |
| | | 339 | | scoped ReadOnlySpan<ulong> root, int orientation, bool sidePoint, scoped ref ConePlaneRaySelection selection) |
| | | 340 | | { |
| | 13684 | 341 | | if (!frame.Normal.X.IsZero || !frame.Normal.Z.IsZero) return false; |
| | 940 | 342 | | bool axisPoint = point.X.Sign(root) == 0 && point.Z.Sign(root) == 0; |
| | 940 | 343 | | int direction = frame.Normal.Y.Sign * orientation; |
| | 965 | 344 | | if (!sidePoint && direction < 0 && !axisPoint) return false; |
| | 1830 | 345 | | Span<ulong> work = stackalloc ulong[4 * frame.FieldWords]; Span<int> signs = stackalloc int[4]; |
| | 915 | 346 | | ContactQuadratic n = ContactQuadratic.At(work, signs, 0, frame.FieldWords); |
| | 915 | 347 | | ContactQuadratic d = ContactQuadratic.At(work, signs, 1, frame.FieldWords); |
| | | 348 | | // Along the axis the widening direction stays inside until the base. |
| | | 349 | | // A point on the axis reaches either axial endpoint (also for R=0). |
| | | 350 | | // A nonaxial side point exits immediately in the narrowing direction. |
| | 915 | 351 | | if (direction > 0 || axisPoint) |
| | | 352 | | { |
| | 648 | 353 | | bool found = GetAxialExit(point, frame, orientation, n, d); |
| | | 354 | | System.Diagnostics.Debug.Assert(found); |
| | | 355 | | } |
| | 534 | 356 | | else { n.Clear(); d.Set(Signed576.ExtendValue(Signed320.ExtendValue(Signed192.Signed(1)))); } |
| | 915 | 357 | | selection.Keep(n, d, root, frame, point); |
| | 915 | 358 | | return true; |
| | | 359 | | } |
| | | 360 | | |
| | | 361 | | private static bool GetAxialExit(scoped ConePlaneRayPoint point, in ConePlaneRayFrame frame, |
| | | 362 | | int orientation, scoped ContactQuadratic numerator, scoped ContactQuadratic denominator) |
| | | 363 | | { |
| | 7045 | 364 | | int direction = frame.Normal.Y.Sign * orientation; |
| | 7045 | 365 | | if (direction == 0) |
| | 1430 | 366 | | return false; |
| | 5615 | 367 | | Signed576 ny = frame.Normal.Y; |
| | 5615 | 368 | | if (direction > 0) |
| | | 369 | | { |
| | 2942 | 370 | | ContactQuadratic.Scale(point.Denominator, Signed576.ExtendValue(frame.Finite.FullHeight), numerator); |
| | 2942 | 371 | | numerator.Add(point.Y, -1); |
| | | 372 | | } |
| | 2673 | 373 | | else point.Y.CopyTo(numerator); |
| | 8875 | 374 | | if (ny.Sign < 0) ny = WideArithmetic.SubtractSigned576(default, ny); |
| | 5615 | 375 | | ContactQuadratic.Scale(point.Denominator, ny, denominator); |
| | 5615 | 376 | | return true; |
| | | 377 | | } |
| | | 378 | | |
| | | 379 | | private static void KeepFirstExit(scoped ContactQuadratic numerator, scoped ContactQuadratic denominator, |
| | | 380 | | scoped ReadOnlySpan<ulong> root, ContactQuadratic bestN, ContactQuadratic bestD, |
| | | 381 | | Span<ulong> bestRoot, ref bool hasValue) |
| | | 382 | | { |
| | | 383 | | // The three private constructions establish these signs: admitted |
| | | 384 | | // quadratic roots, inward side crossings, and finite axial endpoints. |
| | | 385 | | #if DEBUG |
| | | 386 | | System.Diagnostics.Debug.Assert(numerator.Sign(root) >= 0 && denominator.Sign(root) > 0); |
| | | 387 | | #endif |
| | 10964 | 388 | | if (hasValue && ContactQuadratic.CompareRatios(numerator, denominator, root, bestN, bestD, bestRoot) >= 0) |
| | 3626 | 389 | | return; |
| | 14676 | 390 | | numerator.CopyTo(bestN); denominator.CopyTo(bestD); |
| | 22014 | 391 | | root.CopyTo(bestRoot); bestRoot[root.Length..].Clear(); hasValue = true; |
| | 7338 | 392 | | } |
| | | 393 | | } |