nx_fwddiff_equiv_gate.nx source
↩ module page · 216 lines · 11345 B
1// nx_fwddiff_equiv_gate.nx -- THE ARITHMETIC CORE OF THE FORWARD-DIFFERENCED CONIC EVALUATION,
2// PROVEN INDEPENDENT OF THE RENDERER.
3//
4// WHY THIS GATE EXISTS. The splat rasterizer's inner loop evaluated, per pixel per splat:
5// pnum = cc*dx*dx - 2*cb*dx*dy + ca*dy*dy ; k = pnum * MUL / det
6// which is six multiplies and ONE 64-BIT INTEGER DIVISION on the hottest path in the estate --
7// measured at 25.9 ns per conic evaluation, an idiv's price rather than a multiply's. The published
8// remedy (SHARP-GS, SIGGRAPH 2026: "Forward Differencing, which replaces expensive probability
9// density function evaluations with efficient incremental updates") is to walk the quadratic by its
10// constant second difference and to carry the quotient in a bracket instead of dividing.
11//
12// THAT REWRITE IS ONLY ADMISSIBLE IF IT IS EXACT. A renderer that is faster and changes one pixel is
13// a broken renderer. Two identities carry the whole change, and this gate proves BOTH over a swept
14// domain rather than asserting them from the algebra:
15//
16// IDENTITY 1 (forward difference): walking P by P += D ; D += DD reproduces
17// P(dx) = MUL*(cc*dx^2 - 2*cb*dx*dy + ca*dy^2) exactly,
18// because P is quadratic in dx and DD = 2*cc*MUL is its
19// constant second difference. Integer arithmetic, so "exactly"
20// means bit-for-bit, not "to within rounding".
21//
22// IDENTITY 2 (bracket = division): maintaining k*det <= P < (k+1)*det by stepping k yields
23// exactly floor(P/det) = the integer the division produced,
24// for every P >= 0 and det >= 1.
25//
26// THE SUBJECT IS THE IDENTITY, NOT THE RENDERER, so this gate deliberately imports NO renderer
27// constants: it sweeps the multiplier, the LUT cut and the covariance terms, and our shipped values
28// are one point inside that sweep. A gate that mirrored GLUTU/GEXPN would be a second copy of them --
29// this one is strictly stronger and owns nothing it could drift from.
30//
31// THE SWEEP IS DERIVED FROM THE GEOMETRY IT MODELS, never picked: ca and cc are built as
32// e^2 + e^2 + 1 exactly as the projector builds them, cb is swept across the Cauchy-Schwarz bound
33// that makes the form positive definite, and the pixel rectangle is the 3-sigma box the rasterizer
34// actually walks.
35
36import "nx_syscalls.nx"
37import "nx_gate_verdict.nx"
38
39// Sweep extents. Each is a COUNT of cases to visit, not a limit on what is representable: the gate
40// walks every combination it declares and prints the total, so a shrinking sweep cannot hide.
41const FD_NE: i64 = 6 // distinct projected-axis magnitudes
42const FD_NCB: i64 = 5 // cross-term steps across the positive-definite range
43const FD_NDY: i64 = 7 // scanline offsets from the splat centre
44const FD_NMUL: i64 = 3 // LUT sub-step multipliers (ours is one of these)
45const FD_SIGMA: i64 = 3 // AABB half-extent in sigma -- the box the rasterizer walks
46
47func fd_isqrt(v: i64) -> i64 {
48 if v <= 0 { return 0 }
49 var x: i64 = v
50 var y: i64 = (x + 1) / 2
51 while y < x { x = y; y = (x + v / x) / 2 }
52 return x
53}
54
55func fd_e_at(i: i64) -> i64 {
56 if i == 0 { return 1 }
57 if i == 1 { return 3 }
58 if i == 2 { return 11 }
59 if i == 3 { return 40 }
60 if i == 4 { return 137 }
61 return 512
62}
63func fd_mul_at(i: i64) -> i64 {
64 if i == 0 { return 1 }
65 if i == 1 { return 16 } // the shipped LUT sub-step; present as a POINT in the sweep, not as a mirror
66 return 64
67}
68
69func main() -> i64 {
70 let ctr: *i64 = gv_ctr()
71 gv_head("nx_fwddiff_equiv_gate -- the incremental conic evaluation is EXACTLY the division it replaces" as *u8)
72
73 var pixels: i64 = 0 // pixels visited across the whole sweep
74 var kmoves: i64 = 0 // times the bracket actually stepped -- the anti-vacuity signal
75 var kdistinct: i64 = 0 // scanlines on which k took more than one value
76 var negp: i64 = 0 // pixels where P < 0 (the positive-definiteness check)
77 var bad_p: i64 = 0 // forward-difference disagreements
78 var bad_k: i64 = 0 // bracket-vs-division disagreements
79 var cases: i64 = 0
80
81 var ie: i64 = 0
82 while ie < FD_NE {
83 let e1: i64 = fd_e_at(ie)
84 var je: i64 = 0
85 while je < FD_NE {
86 let e2: i64 = fd_e_at(je)
87 // ca and cc are built the way the projector builds them: sums of squared projected axes,
88 // plus one. That "+1" is what makes the form strictly positive definite.
89 let ca: i64 = e1*e1 + e2*e2 + 1
90 let cc: i64 = e2*e2 + e1*e1 + 1
91 // Cauchy-Schwarz gives cb^2 <= (ca-1)(cc-1); sweep cb across that range, both signs.
92 let cbmax: i64 = fd_isqrt((ca-1)*(cc-1))
93 var icb: i64 = 0
94 while icb < FD_NCB {
95 var cb: i64 = 0
96 if cbmax > 0 { cb = 0 - cbmax + 2*cbmax*icb/(FD_NCB-1) }
97 let det: i64 = ca*cc - cb*cb
98 if det >= 1 {
99 var imul: i64 = 0
100 while imul < FD_NMUL {
101 let mul: i64 = fd_mul_at(imul)
102 let rx: i64 = FD_SIGMA * fd_isqrt(ca)
103 let ry: i64 = FD_SIGMA * fd_isqrt(cc)
104 var idy: i64 = 0
105 while idy < FD_NDY {
106 var dy: i64 = 0
107 if FD_NDY > 1 { dy = 0 - ry + 2*ry*idy/(FD_NDY-1) }
108 cases = cases + 1
109
110 let dx0: i64 = 0 - rx
111 // forward-difference state, exactly as the rasterizer initialises it
112 var p: i64 = (cc*dx0*dx0 - 2*cb*dx0*dy + ca*dy*dy) * mul
113 var d: i64 = (cc*(2*dx0 + 1) - 2*cb*dy) * mul
114 let dd: i64 = 2 * cc * mul
115 var k: i64 = 0
116 var kc: i64 = 0
117 var kn: i64 = det
118 var kfirst: i64 = 0 - 1
119 var kvaried: i64 = 0
120
121 var dx: i64 = dx0
122 while dx <= rx {
123 // IDENTITY 1: the incremental P against the direct quadratic
124 let pdirect: i64 = (cc*dx*dx - 2*cb*dx*dy + ca*dy*dy) * mul
125 if p != pdirect { bad_p = bad_p + 1 }
126 if pdirect < 0 { negp = negp + 1 }
127
128 // IDENTITY 2: the bracket against the division it replaces
129 if p >= 0 {
130 while kc > p { kn = kc; kc = kc - det; k = k - 1; kmoves = kmoves + 1 }
131 while kn <= p { kc = kn; kn = kn + det; k = k + 1; kmoves = kmoves + 1 }
132 let kdiv: i64 = p / det
133 if k != kdiv { bad_k = bad_k + 1 }
134 if kfirst < 0 { kfirst = k }
135 if k != kfirst { kvaried = 1 }
136 }
137
138 pixels = pixels + 1
139 p = p + d
140 d = d + dd
141 dx = dx + 1
142 }
143 if kvaried == 1 { kdistinct = kdistinct + 1 }
144 idy = idy + 1
145 }
146 imul = imul + 1
147 }
148 }
149 icb = icb + 1
150 }
151 je = je + 1
152 }
153 ie = ie + 1
154 }
155
156 gv_puts(" swept cases=" as *u8); gv_num(cases)
157 gv_puts(" pixels=" as *u8); gv_num(pixels)
158 gv_puts(" bracket_steps=" as *u8); gv_num(kmoves)
159 gv_puts(" scanlines_with_varying_k=" as *u8); gv_num(kdistinct)
160 gv_puts("\n forward_diff_mismatches=" as *u8); gv_num(bad_p)
161 gv_puts(" bracket_vs_division_mismatches=" as *u8); gv_num(bad_k)
162 gv_puts(" pixels_with_negative_P=" as *u8); gv_num(negp)
163 gv_puts("\n\n" as *u8)
164
165 // The fixture must have DONE something before any verdict about it means anything.
166 gv_check("fixture-visited-pixels (count bound in condition: an empty sweep cannot pass)" as *u8, pixels > 0, ctr)
167 gv_check("fixture-swept-more-than-one-case" as *u8, cases > 1, ctr)
168
169 // ANTI-VACUITY: if the bracket never stepped, or k never varied along a scanline, then the
170 // identity was tested only where it is trivially true. A sweep that never moves k proves nothing
171 // about a walk whose entire purpose is to move k.
172 gv_check("anti-vacuity-bracket-actually-stepped" as *u8, kmoves > pixels / 2, ctr)
173 gv_check("anti-vacuity-k-varied-within-scanlines" as *u8, kdistinct > cases / 2, ctr)
174
175 // The two identities.
176 gv_check("forward-difference-reproduces-the-direct-quadratic-exactly" as *u8, bad_p == 0, ctr)
177 gv_check("bracket-equals-integer-division-exactly" as *u8, bad_k == 0, ctr)
178
179 // The positive-definiteness claim the rewrite rests on, MEASURED rather than assumed: over a
180 // sweep built the way the projector builds its terms, P is never negative. The renderer keeps a
181 // division fallback for P<0 anyway, so this is corroboration of an argument, not a load-bearing
182 // assumption -- and it is reported as a number so a future change that breaks it is visible.
183 gv_check("positive-definite-form-never-produced-a-negative-P" as *u8, negp == 0, ctr)
184
185 // NEG-CONTROL: a deliberately corrupted increment must be CAUGHT by identity 1. If a wrong
186 // second difference still passed, the tooth above would be decoration.
187 var badcase_detected: i64 = 0
188 var goodcase_detected: i64 = 0
189 let tca: i64 = 137*137 + 40*40 + 1
190 let tcc: i64 = 40*40 + 137*137 + 1
191 let tcb: i64 = fd_isqrt((tca-1)*(tcc-1)) / 2
192 let tmul: i64 = 16
193 let trx: i64 = FD_SIGMA * fd_isqrt(tca)
194 let tdy: i64 = trx / 3
195 var q: i64 = ((tcc*trx*trx) + 2*tcb*trx*tdy + tca*tdy*tdy) * tmul
196 var qd: i64 = (tcc*(0 - 2*trx + 1) - 2*tcb*tdy) * tmul
197 let qdd_good: i64 = 2 * tcc * tmul
198 let qdd_bad: i64 = qdd_good + 1 // ONE off: the smallest possible corruption
199 var qbad: i64 = q
200 var qdbad: i64 = qd
201 var dxx: i64 = 0 - trx
202 while dxx <= trx {
203 let direct: i64 = (tcc*dxx*dxx - 2*tcb*dxx*tdy + tca*tdy*tdy) * tmul
204 if q != direct { goodcase_detected = goodcase_detected + 1 }
205 if qbad != direct { badcase_detected = badcase_detected + 1 }
206 q = q + qd; qd = qd + qdd_good
207 qbad = qbad + qdbad; qdbad = qdbad + qdd_bad
208 dxx = dxx + 1
209 }
210 gv_puts(" neg-control: corrupt-second-difference mismatches=" as *u8); gv_num(badcase_detected)
211 gv_puts(" correct-second-difference mismatches=" as *u8); gv_num(goodcase_detected)
212 gv_puts("\n\n" as *u8)
213 gv_bite("neg-control-corrupted-second-difference-is-caught" as *u8, badcase_detected > 0, goodcase_detected > 0, ctr)
214
215 return gv_verdict("NX-FWDDIFF-EQUIV" as *u8, ctr, "the incremental walk returns the divisions integer, so the rewrite cannot move a pixel" as *u8)
216}