nx_henyey_greenstein.nx source
↩ module page · 265 lines · 13356 B
1// nx_henyey_greenstein.nx -- anisotropic single-scatter phase function in Q10.
2//
3// Epoch 2 FIRMAMENT primitive per nxc2/docs/NISHI_GAME_ENGINE_ROADMAP.md.
4//
5// The Henyey-Greenstein phase function approximates light scattering
6// off small particles (atmospheric aerosols, water droplets, fog,
7// pigment). Published in Henyey & Greenstein 1941 "Diffuse Radiation
8// in the Galaxy" (Astrophysical Journal 93, 70-83) -- pre-1970
9// astrophysics, public-domain math.
10//
11// Formula:
12// P(cos θ, g) = (1 - g^2) / (4π · (1 + g^2 - 2g·cos θ)^1.5)
13//
14// where:
15// cos θ in [-1, 1]: scattering angle cosine (1 = forward, -1 = backward)
16// g in [-1, 1]: asymmetry parameter
17// g > 0 = forward-scattering (atmospheric haze ~ 0.8)
18// g = 0 = isotropic
19// g < 0 = backward-scattering (some pigments)
20//
21// Output: probability density. Normalized so ∫P dΩ = 1 over sphere.
22// Range can exceed 1.0 when |g| near 1 and cos θ near sign(g)
23// (forward/backward peaking). Capped at NX_HG_SATURATED to keep
24// downstream arithmetic well-behaved.
25//
26// Fixed-point strategy: keep numerator and denominator in raw
27// integer space and defer all divisions to a single final step.
28// Avoids the Q10 precision starvation that would otherwise hit the
29// (1 + g^2 - 2gc) term as it approaches 0.
30//
31// Loss audit: the only information-discarding ops are (a) integer
32// sqrt (_isqrt) which truncates the fractional part, introducing
33// ~1% error at typical g values; (b) clamp at |g| >= NX_HG_G_MAX
34// to avoid singularity; (c) final integer division. All three are
35// inherent to fixed-point evaluation of a singular function and are
36// named here per cardinal rule 25.
37//
38// genealogy_id: henyey_greenstein_1941_phase
39// lineage_id: nx_henyey_greenstein_q10
40
41// nx_safety_envelope:
42// intended_use: AUTO_APPLIED -- primitive-specific tuning queued
43// sil_target: SIL1
44// evidence: [bulk_applied_2026-05-16, see-file-comment-for-detail]
45// verdict: NOT_YET_EVALUATED
46
47import "nx_syscalls.nx"
48import "nx_tier.nx"
49const NX_MAGIC_1048576: i64 = 1048576
50const NX_MAGIC_1024: i64 = 1024
51const NX_MAGIC_1300: i64 = 1300
52const NX_MAGIC_1700: i64 = 1700
53const NX_MAGIC_9999: i64 = 9999
54
55// ===== Q10 constants ================================================
56// Q-scale (1.0 in Q10).
57const NX_HG_Q: nx_int = 1024
58
59// 4π in Q10: 4 * 3.14159265 * 1024 = 12867.96 -> 12868
60// Rounding error ~5e-6 -- well below other approximations.
61const NX_HG_FOUR_PI_Q10: nx_int = 12868
62
63// Maximum |g| supported. At |g| = Q (1.0 exactly), the function is
64// a delta -- infinity at cos θ = sign(g), zero elsewhere -- and
65// purely-integer evaluation collapses to a 0/0 form. We clamp at
66// 1023 to keep numerical evaluation bounded; callers needing a true
67// delta should use a different primitive (planned: nx_phase_dirac).
68const NX_HG_G_MAX: nx_int = 1023
69
70// Saturation value returned when the denominator collapses to <= 0
71// due to the singularity at |g| -> 1. Chosen large enough to be
72// distinguishable from any well-defined output without overflowing
73// downstream Q10 arithmetic that composes phase function values.
74const NX_HG_SATURATED: nx_int = 100000
75
76// ===== Sealed-enum verdict bands ====================================
77// Per dual-reading cardinal -- every primitive ships Q10 + qualitative.
78// 7-band classifier on g (the asymmetry parameter) -- the science
79// uses sealed bins because "moderate forward" vs "strong forward" is
80// the level at which atmospheric haze vs cloud droplet differs.
81const NX_HG_BAND_STRONG_BACKWARD: nx_int = 0 // g <= -700
82const NX_HG_BAND_MODERATE_BACKWARD: nx_int = 1 // -700 < g <= -300
83const NX_HG_BAND_WEAK_BACKWARD: nx_int = 2 // -300 < g <= -50
84const NX_HG_BAND_ISOTROPIC: nx_int = 3 // -50 < g < 50
85const NX_HG_BAND_WEAK_FORWARD: nx_int = 4 // 50 <= g < 300
86const NX_HG_BAND_MODERATE_FORWARD: nx_int = 5 // 300 <= g < 700
87const NX_HG_BAND_STRONG_FORWARD: nx_int = 6 // g >= 700
88
89const NX_HG_THRESH_STRONG_BACK: nx_int = -700
90const NX_HG_THRESH_MOD_BACK: nx_int = -300
91const NX_HG_THRESH_WEAK_BACK: nx_int = -50
92const NX_HG_THRESH_WEAK_FWD: nx_int = 50
93const NX_HG_THRESH_MOD_FWD: nx_int = 300
94const NX_HG_THRESH_STRONG_FWD: nx_int = 700
95
96// ===== Validity predicate ==========================================
97func nx_henyey_greenstein_band_is_valid(b: nx_int) -> nx_int {
98 if b == NX_HG_BAND_STRONG_BACKWARD { return 1 }
99 if b == NX_HG_BAND_MODERATE_BACKWARD { return 1 }
100 if b == NX_HG_BAND_WEAK_BACKWARD { return 1 }
101 if b == NX_HG_BAND_ISOTROPIC { return 1 }
102 if b == NX_HG_BAND_WEAK_FORWARD { return 1 }
103 if b == NX_HG_BAND_MODERATE_FORWARD { return 1 }
104 if b == NX_HG_BAND_STRONG_FORWARD { return 1 }
105 return 0
106}
107
108// ===== Sealed-enum classifier on g =================================
109func nx_henyey_greenstein_classify(g_q10: nx_int) -> nx_int {
110 if g_q10 <= NX_HG_THRESH_STRONG_BACK { return NX_HG_BAND_STRONG_BACKWARD }
111 if g_q10 <= NX_HG_THRESH_MOD_BACK { return NX_HG_BAND_MODERATE_BACKWARD }
112 if g_q10 <= NX_HG_THRESH_WEAK_BACK { return NX_HG_BAND_WEAK_BACKWARD }
113 if g_q10 < NX_HG_THRESH_WEAK_FWD { return NX_HG_BAND_ISOTROPIC }
114 if g_q10 < NX_HG_THRESH_MOD_FWD { return NX_HG_BAND_WEAK_FORWARD }
115 if g_q10 < NX_HG_THRESH_STRONG_FWD { return NX_HG_BAND_MODERATE_FORWARD }
116 return NX_HG_BAND_STRONG_FORWARD
117}
118
119// ===== Integer square root ==========================================
120// Digit-by-digit binary sqrt; returns floor(sqrt(n)) for n >= 0.
121// No division, no multiplication except shifts. O(log n) iterations.
122// Reference algorithm: Wikipedia "Methods of computing square roots",
123// the "binary numeral system (base 2)" section -- public-domain math.
124func _isqrt(n: nx_int) -> nx_int {
125 if n <= 0 { return 0 }
126 var x: nx_int = n
127 var c: nx_int = 0
128 var d: nx_int = 1
129
130 // Find the highest power-of-4 not exceeding n.
131 while d <= n / 4 {
132 d = d * 4
133 }
134
135 while d != 0 {
136 if x >= c + d {
137 x = x - c - d
138 c = c / 2 + d
139 } else {
140 c = c / 2
141 }
142 d = d / 4
143 }
144 return c
145}
146
147// ===== Phase function ===============================================
148// cos_theta_q10: Q10, in [-1024, 1024] (clamped if out of range)
149// g_q10: Q10, in [-1023, 1023] (clamped if out of range)
150// returns: Q10 probability density; ranges from very small
151// to NX_HG_SATURATED at the singularity-near case.
152//
153// Strategy: evaluate (Q² - G²) · Q² / (4π · D · sqrt(D))
154// where D = Q² + G² - 2GC is the inner denominator in "Q20" form
155// (i.e. multiplied by Q² so all terms stay integer-clean).
156//
157// This avoids the Q10-truncation problem that hits when 1+g^2-2gc
158// approaches 0 -- by deferring the divide-by-Q² to the very end,
159// we keep all the precision available in 64-bit integers.
160func nx_henyey_greenstein_phase(cos_theta_q10: nx_int, g_q10: nx_int) -> nx_int {
161 var c: nx_int = cos_theta_q10
162 if c > NX_HG_Q { c = NX_HG_Q }
163 if c < (0 - NX_HG_Q) { c = 0 - NX_HG_Q }
164
165 var g: nx_int = g_q10
166 if g > NX_HG_G_MAX { g = NX_HG_G_MAX }
167 if g < (0 - NX_HG_G_MAX) { g = 0 - NX_HG_G_MAX }
168
169 let q2: nx_int = NX_HG_Q * NX_HG_Q // NX_MAGIC_1048576
170 let g_sq: nx_int = g * g // max ~ 1023^2 ~ 1.05e6
171 let two_gc: nx_int = 2 * g * c // max ~ 2 * 1023 * NX_MAGIC_1024 ~ 2.1e6
172
173 let num: nx_int = q2 - g_sq // (Q^2 - G^2), always >= 0 since |g| < Q
174 let d: nx_int = q2 + g_sq - two_gc // (Q^2 + G^2 - 2GC), >= 0 ideally
175
176 if d <= 0 { return NX_HG_SATURATED }
177
178 let sqrt_d: nx_int = _isqrt(d) // approx sqrt(D) in same scale as D
179 let denom_15: nx_int = d * sqrt_d // D * sqrt(D) ~ D^1.5 (raw)
180
181 if denom_15 <= 0 { return NX_HG_SATURATED }
182
183 // P_q10 = (Q^2 - G^2) * Q^2 / (4π · D · sqrt(D))
184 let scaled_num: nx_int = num * q2
185 let scaled_denom: nx_int = NX_HG_FOUR_PI_Q10 * denom_15
186
187 if scaled_denom <= 0 { return NX_HG_SATURATED }
188
189 let p_q10: nx_int = scaled_num / scaled_denom
190 if p_q10 > NX_HG_SATURATED { return NX_HG_SATURATED }
191 return p_q10
192}
193
194// ===== Self-test ====================================================
195func main() -> i64 {
196 // T1: Integer sqrt sanity.
197 if _isqrt(0) != 0 { return __syscall(93, 1, 0, 0, 0, 0, 0) }
198 if _isqrt(1) != 1 { return __syscall(93, 2, 0, 0, 0, 0, 0) }
199 if _isqrt(4) != 2 { return __syscall(93, 3, 0, 0, 0, 0, 0) }
200 if _isqrt(NX_MAGIC_1024) != 32 { return __syscall(93, 4, 0, 0, 0, 0, 0) }
201 if _isqrt(NX_MAGIC_1048576) != NX_MAGIC_1024 { return __syscall(93, 5, 0, 0, 0, 0, 0) }
202 // Floor-rounded: sqrt(99) = 9.94, floor -> 9.
203 if _isqrt(99) != 9 { return __syscall(93, 6, 0, 0, 0, 0, 0) }
204
205 // T2: Isotropic (g = 0) returns 1 / (4π) regardless of angle.
206 // 1/(4π) = 0.07958 -> Q10 = 81.5, integer-sqrt route gives ~79.
207 let p_iso_fwd: nx_int = nx_henyey_greenstein_phase(NX_HG_Q, 0)
208 let p_iso_back: nx_int = nx_henyey_greenstein_phase(0 - NX_HG_Q, 0)
209 let p_iso_side: nx_int = nx_henyey_greenstein_phase(0, 0)
210 if p_iso_fwd < 75 { return __syscall(93, 10, 0, 0, 0, 0, 0) }
211 if p_iso_fwd > 85 { return __syscall(93, 11, 0, 0, 0, 0, 0) }
212 if p_iso_back < 75 { return __syscall(93, 12, 0, 0, 0, 0, 0) }
213 if p_iso_back > 85 { return __syscall(93, 13, 0, 0, 0, 0, 0) }
214 if p_iso_side < 75 { return __syscall(93, 14, 0, 0, 0, 0, 0) }
215 if p_iso_side > 85 { return __syscall(93, 15, 0, 0, 0, 0, 0) }
216
217 // T3: Strong forward (g = 0.7). At cos θ = 1, analytic P ~= 1.504.
218 // P_q10 expected ~= 1540; allow ±10% for integer-sqrt approximation.
219 let p_fwd_peak: nx_int = nx_henyey_greenstein_phase(NX_HG_Q, 717)
220 if p_fwd_peak < NX_MAGIC_1300 { return __syscall(93, 16, 0, 0, 0, 0, 0) }
221 if p_fwd_peak > NX_MAGIC_1700 { return __syscall(93, 17, 0, 0, 0, 0, 0) }
222
223 // T4: Strong forward (g = 0.7), at cos θ = -1 (backward direction)
224 // -- minimum value. Analytic P ~= 0.0089 -> Q10 ~= 9.
225 let p_fwd_back: nx_int = nx_henyey_greenstein_phase(0 - NX_HG_Q, 717)
226 if p_fwd_back < 7 { return __syscall(93, 18, 0, 0, 0, 0, 0) }
227 if p_fwd_back > 12 { return __syscall(93, 19, 0, 0, 0, 0, 0) }
228
229 // T5: Symmetry under (g, c) -> (-g, -c). Phase function obeys
230 // P(c, g) == P(-c, -g) exactly (1-g^2 invariant, denom invariant).
231 let p_a: nx_int = nx_henyey_greenstein_phase( 512, 500)
232 let p_b: nx_int = nx_henyey_greenstein_phase(0 - 512, 0 - 500)
233 if p_a != p_b { return __syscall(93, 20, 0, 0, 0, 0, 0) }
234
235 // T6: Clamping at |g| = NX_HG_G_MAX -- input g > Q must return
236 // the same value as g = NX_HG_G_MAX (1023).
237 let p_at_max: nx_int = nx_henyey_greenstein_phase(0, NX_HG_G_MAX)
238 let p_overshoot: nx_int = nx_henyey_greenstein_phase(0, NX_MAGIC_9999)
239 if p_at_max != p_overshoot { return __syscall(93, 21, 0, 0, 0, 0, 0) }
240
241 // T7: Classifier band boundaries.
242 if nx_henyey_greenstein_classify(0 - 800) != NX_HG_BAND_STRONG_BACKWARD { return __syscall(93, 30, 0, 0, 0, 0, 0) }
243 if nx_henyey_greenstein_classify(0 - 700) != NX_HG_BAND_STRONG_BACKWARD { return __syscall(93, 31, 0, 0, 0, 0, 0) }
244 if nx_henyey_greenstein_classify(0 - 500) != NX_HG_BAND_MODERATE_BACKWARD { return __syscall(93, 32, 0, 0, 0, 0, 0) }
245 if nx_henyey_greenstein_classify(0 - 100) != NX_HG_BAND_WEAK_BACKWARD { return __syscall(93, 33, 0, 0, 0, 0, 0) }
246 if nx_henyey_greenstein_classify( 0) != NX_HG_BAND_ISOTROPIC { return __syscall(93, 34, 0, 0, 0, 0, 0) }
247 if nx_henyey_greenstein_classify( 100) != NX_HG_BAND_WEAK_FORWARD { return __syscall(93, 35, 0, 0, 0, 0, 0) }
248 if nx_henyey_greenstein_classify( 500) != NX_HG_BAND_MODERATE_FORWARD { return __syscall(93, 36, 0, 0, 0, 0, 0) }
249 if nx_henyey_greenstein_classify( 800) != NX_HG_BAND_STRONG_FORWARD { return __syscall(93, 37, 0, 0, 0, 0, 0) }
250
251 // T8: Validity predicate accepts every enum value and rejects
252 // out-of-band integers.
253 if nx_henyey_greenstein_band_is_valid(NX_HG_BAND_STRONG_BACKWARD) != 1 { return __syscall(93, 40, 0, 0, 0, 0, 0) }
254 if nx_henyey_greenstein_band_is_valid(NX_HG_BAND_MODERATE_BACKWARD) != 1 { return __syscall(93, 41, 0, 0, 0, 0, 0) }
255 if nx_henyey_greenstein_band_is_valid(NX_HG_BAND_WEAK_BACKWARD) != 1 { return __syscall(93, 42, 0, 0, 0, 0, 0) }
256 if nx_henyey_greenstein_band_is_valid(NX_HG_BAND_ISOTROPIC) != 1 { return __syscall(93, 43, 0, 0, 0, 0, 0) }
257 if nx_henyey_greenstein_band_is_valid(NX_HG_BAND_WEAK_FORWARD) != 1 { return __syscall(93, 44, 0, 0, 0, 0, 0) }
258 if nx_henyey_greenstein_band_is_valid(NX_HG_BAND_MODERATE_FORWARD) != 1 { return __syscall(93, 45, 0, 0, 0, 0, 0) }
259 if nx_henyey_greenstein_band_is_valid(NX_HG_BAND_STRONG_FORWARD) != 1 { return __syscall(93, 46, 0, 0, 0, 0, 0) }
260 if nx_henyey_greenstein_band_is_valid(0 - 1) != 0 { return __syscall(93, 47, 0, 0, 0, 0, 0) }
261 if nx_henyey_greenstein_band_is_valid(7) != 0 { return __syscall(93, 48, 0, 0, 0, 0, 0) }
262 if nx_henyey_greenstein_band_is_valid(99) != 0 { return __syscall(93, 49, 0, 0, 0, 0, 0) }
263
264 return 0
265}