nx_quadrature_gate.nx source
↩ module page · 71 lines · 5788 B
1// nx_quadrature_gate.nx -- computational CM4: adaptive quadrature with a MEASURED error bound that REFUSES rather
2// than lie. In-process over nx_quadrature (the one ruler), analytic KATs whose exact values are closed forms:
3// x^2 on [0,1] = 1/3 · x^3 on [0,2] = 4 (Simpson-exact: must cost the minimum evaluations) · 1/(1+x^2) on [0,1] =
4// pi/4 · a sharp peak 1/(0.01+x^2) on [-1,1] = 20 atan(10) (forces deep subdivision near 0) · |x| on [-1,1] = 1
5// (a kink the adaptive split must resolve).
6// Honesty teeth: |value - exact| <= the reported error bound on EVERY KAT; a tolerance below binary64 resolution
7// on the result's scale returns NQ_REFUSED_TOL with the best value still carried; a > b and tol <= 0 are refused as
8// bad input; an integrand that returns NaN is refused as NONFINITE. The refusal teeth are the neg-controls: an
9// implementation that always returns OK cannot pass them. license_tier: ORIGINAL
10import "nx_syscalls.nx"
11import "nx_gate_verdict.nx"
12import "nx_quadrature.nx"
13
14// The KAT constants are FUNCTIONS, not module consts: the const-expression parser accepts only integer literals
15// today (lang LN35, found by this gate) and the lexer has no exponent form (LN34), so 1e-9 is spelled in decimals.
16func qg_pi_4() -> f64 { return 0.7853981633974483 }
17func qg_peak_exact() -> f64 { return 29.42255348607467 } // 20 * atan(10)
18func qg_tol() -> f64 { return 0.000000001 } // 1e-9 absolute, the KAT tolerance
19// 1e-20 built by multiplication: a 20-digit fraction literal overflows the exact-rational packer (10^20 > i64) and the
20// compiler LOOPS on it -- lang LN36, found by this line; the exponent form 1e-20 is LN34.
21func qg_tol_impossible() -> f64 { return 0.000000001 * 0.000000001 * 0.01 }
22func qg_slack() -> f64 { return 0.000000000001 } // 1e-12 rounding slack allowed on top of the reported bound
23
24func qg_sq(x: f64) -> f64 { return x * x }
25func qg_cube(x: f64) -> f64 { return x * x * x }
26func qg_lorentz(x: f64) -> f64 { return 1.0 / (1.0 + x * x) }
27func qg_peak(x: f64) -> f64 { return 1.0 / (0.01 + x * x) }
28func qg_absx(x: f64) -> f64 { if x < 0.0 { return 0.0 - x } return x }
29func qg_nan(x: f64) -> f64 { let z: f64 = 0.0; return z / z }
30
31func qg_abs(x: f64) -> f64 { if x < 0.0 { return 0.0 - x } return x }
32// 1 iff |value - exact| <= tol AND |value - exact| <= reported bound + slack (the bound is honest)
33func qg_kat(name: *u8, f: func(f64) -> f64, a: f64, b: f64, exact: f64, c: *i64) -> i64 {
34 let out: *f64 = sys_mmap(NQ_OUT_BYTES) as *f64
35 let rc: i64 = ni_quad_adaptive(f, a, b, qg_tol(), out)
36 let err: f64 = qg_abs(out[NQ_OUT_VALUE] - exact)
37 var ok: i64 = 0
38 if rc == NQ_OK { if err <= qg_tol() { if err <= out[NQ_OUT_ERR] + qg_slack() { ok = 1 } } }
39 gv_puts(" " as *u8); gv_puts(name); gv_puts(": rc=" as *u8); gv_num(rc)
40 gv_puts(" evals=" as *u8); gv_num(out[NQ_OUT_EVALS] as i64)
41 gv_puts(" err_x1e12=" as *u8); gv_num((err * 1000000000000.0) as i64)
42 gv_puts(" bound_x1e12=" as *u8); gv_num((out[NQ_OUT_ERR] * 1000000000000.0) as i64); gv_puts("\n" as *u8)
43 gv_check(name, ok, c)
44 return out[NQ_OUT_EVALS] as i64
45}
46
47func main() -> i64 {
48 let c: *i64 = gv_ctr()
49 gv_puts("=== NX-QUADRATURE-GATE: adaptive Simpson with a Richardson error bound; refuses rather than lies ===\n" as *u8)
50 qg_kat("kat-x2-on-0-1-equals-one-third-within-1e-9-and-inside-its-own-bound" as *u8, qg_sq, 0.0, 1.0, 0.3333333333333333, c)
51 let ev3: i64 = qg_kat("kat-x3-on-0-2-equals-4-simpson-exact" as *u8, qg_cube, 0.0, 2.0, 4.0, c)
52 gv_check("simpson-exact-cubic-costs-the-minimum-5-evaluations (no needless subdivision)" as *u8, ev3 == 5, c)
53 qg_kat("kat-lorentzian-on-0-1-equals-pi-over-4" as *u8, qg_lorentz, 0.0, 1.0, qg_pi_4(), c)
54 let evp: i64 = qg_kat("kat-sharp-peak-1-over-0.01-plus-x2-on-minus1-1-equals-20-atan-10" as *u8, qg_peak, 0.0 - 1.0, 1.0, qg_peak_exact(), c)
55 gv_check("sharp-peak-forced-real-subdivision (more evaluations than the cubic)" as *u8, evp > 5, c)
56 qg_kat("kat-abs-x-on-minus1-1-equals-1 (a kink is resolved by splitting)" as *u8, qg_absx, 0.0 - 1.0, 1.0, 1.0, c)
57
58 // neg-control: a tolerance no binary64 result of order 1 can meet -> REFUSED-TOLERANCE, value still carried
59 let out: *f64 = sys_mmap(NQ_OUT_BYTES) as *f64
60 let rr: i64 = ni_quad_adaptive(qg_lorentz, 0.0, 1.0, qg_tol_impossible(), out)
61 gv_puts(" impossible-tol rc=" as *u8); gv_num(rr); gv_puts(" evals=" as *u8); gv_num(out[NQ_OUT_EVALS] as i64); gv_puts("\n" as *u8)
62 gv_check("neg-control-impossible-tolerance-is-REFUSED-not-returned-as-converged" as *u8, rr == NQ_REFUSED_TOL, c)
63 gv_check("neg-control-refusal-still-carries-the-best-value-within-1e-9-of-pi-over-4" as *u8, qg_abs(out[NQ_OUT_VALUE] - qg_pi_4()) <= qg_tol(), c)
64 gv_check("neg-control-refusal-reports-a-nonzero-error-bound" as *u8, out[NQ_OUT_ERR] > 0.0, c)
65 // neg-controls: bad input and a non-finite integrand are named, never integrated
66 gv_check("neg-control-reversed-interval-is-BAD-INPUT" as *u8, ni_quad_adaptive(qg_sq, 1.0, 0.0, qg_tol(), out) == NQ_BAD_INPUT, c)
67 gv_check("neg-control-zero-tolerance-is-BAD-INPUT" as *u8, ni_quad_adaptive(qg_sq, 0.0, 1.0, 0.0, out) == NQ_BAD_INPUT, c)
68 gv_check("neg-control-nan-integrand-is-NONFINITE" as *u8, ni_quad_adaptive(qg_nan, 0.0, 1.0, qg_tol(), out) == NQ_NONFINITE, c)
69 gv_check("empty-interval-is-OK-with-value-zero" as *u8, (ni_quad_adaptive(qg_sq, 2.0, 2.0, qg_tol(), out) == NQ_OK) & (out[NQ_OUT_VALUE] == 0.0), c)
70 return gv_verdict("QUADRATURE-GATE" as *u8, c, "CM4: adaptive Simpson with a measured Richardson bound; every KAT lands inside its own reported error, an impossible tolerance is refused with the best value carried, bad input and NaN are named" as *u8)
71}