nx_quadrature.nx source
↩ module page · 114 lines · 5900 B
1// nx_quadrature.nx -- ADAPTIVE QUADRATURE WITH A MEASURED ERROR BOUND (computational CM4, 2026-09-03).
2//
3// ni_quad_adaptive integrates f over [a, b] to a requested absolute tolerance by RECURSIVE SUBDIVISION driven by
4// a MEASURED local error: on each interval it compares Simpson over the whole interval with Simpson over the two
5// halves; the Richardson estimate |S_left + S_right - S_whole| / 15 is the local error of the refined value
6// (Simpson's rule is exact through cubics, so its error scales as h^5 and halving the step divides it by 16).
7// An interval whose estimate meets its share of the tolerance is accepted at the refined value plus the
8// Richardson correction; otherwise it is split and each half gets half the tolerance.
9//
10// THE CONTRACT THE MATRIX ROW NAMES: this organ REFUSES rather than silently returning a wrong answer. When the
11// requested tolerance cannot be met on the interval -- the recursion reaches its depth bound or the evaluation
12// budget with the estimate still above tolerance, or the integrand returns a non-finite value -- the status
13// says so and the caller receives the best value AND the achieved error bound, never a bare number that looks
14// converged. A tolerance below what f64 rounding can resolve on the result's scale is the textbook case: the
15// estimate floors at rounding noise and the organ reports REFUSED-TOLERANCE instead of recursing to a lie.
16//
17// Bounds are DERIVED, not tuned (rule 11): NQ_MAX_DEPTH = 50 subdivisions halves an interval 2^50 times, past
18// which the midpoint is not representable apart from its endpoints in binary64 (53-bit significand); the
19// evaluation budget follows from it. Every call reports its evaluation count so a caller can price it.
20//
21// API (one ruler, every caller composes it):
22// ni_quad_adaptive(f, a, b, tol, out) -> status; out[0] = value, out[1] = error bound, out[2] = evaluations
23// status NQ_OK 0 | NQ_REFUSED_TOL 1 | NQ_BAD_INPUT 2 (a > b, tol <= 0, or out null) | NQ_NONFINITE 3
24// Pure f64 arithmetic on the sovereign compiler (LN16). license_tier: ORIGINAL lib (no main)
25import "nx_syscalls.nx"
26
27const NQ_OK: i64 = 0
28const NQ_REFUSED_TOL: i64 = 1
29const NQ_BAD_INPUT: i64 = 2
30const NQ_NONFINITE: i64 = 3
31const NQ_MAX_DEPTH: i64 = 50 // 2^50 halvings: past this a binary64 midpoint collapses onto an endpoint
32const NQ_MAX_EVALS: i64 = 4194304 // 2^22 integrand evaluations: the budget a depth-50 pathological tree may not exceed
33const NQ_OUT_VALUE: i64 = 0
34const NQ_OUT_ERR: i64 = 1
35const NQ_OUT_EVALS: i64 = 2
36const NQ_OUT_SLOTS: i64 = 4
37const NQ_OUT_BYTES: i64 = 32 // NQ_OUT_SLOTS * 8
38
39// finite iff x == x (NaN fails) and x - x is itself a number: for an infinity, inf - inf is NaN. No largest-finite
40// literal is needed (and the lexer has no exponent form yet -- see lang LN34).
41func nq_finite(x: f64) -> i64 {
42 if x != x { return 0 }
43 let d: f64 = x - x
44 if d != d { return 0 }
45 return 1
46}
47func nq_abs(x: f64) -> f64 { if x < 0.0 { return 0.0 - x } return x }
48
49// Simpson over [a, b] given f(a), f(mid), f(b): (b - a) / 6 * (fa + 4 fm + fb)
50func nq_simpson(a: f64, b: f64, fa: f64, fm: f64, fb: f64) -> f64 {
51 return (b - a) / 6.0 * (fa + 4.0 * fm + fb)
52}
53
54// One node of the recursion. st: [0]=evals [1]=nonfinite-seen [2]=refused (depth or budget exhausted with the
55// estimate above tolerance) [3]=accumulated error bound. Returns the interval's refined integral.
56func nq_node(f: func(f64) -> f64, a: f64, b: f64, fa: f64, fm: f64, fb: f64, whole: f64, tol: f64, depth: i64, st: *f64, cnt: *i64) -> f64 {
57 let m: f64 = (a + b) / 2.0
58 let lm: f64 = (a + m) / 2.0
59 let rm: f64 = (m + b) / 2.0
60 let flm: f64 = f(lm)
61 let frm: f64 = f(rm)
62 cnt[0] = cnt[0] + 2
63 if nq_finite(flm) == 0 { cnt[1] = 1; return whole }
64 if nq_finite(frm) == 0 { cnt[1] = 1; return whole }
65 let left: f64 = nq_simpson(a, m, fa, flm, fm)
66 let right: f64 = nq_simpson(m, b, fm, frm, fb)
67 let refined: f64 = left + right
68 let est: f64 = nq_abs(refined - whole) / 15.0
69 if est <= tol {
70 st[3] = st[3] + est
71 return refined + (refined - whole) / 15.0
72 }
73 if depth >= NQ_MAX_DEPTH { cnt[2] = 1; st[3] = st[3] + est; return refined }
74 if cnt[0] >= NQ_MAX_EVALS { cnt[2] = 1; st[3] = st[3] + est; return refined }
75 let half: f64 = tol / 2.0
76 let vl: f64 = nq_node(f, a, m, fa, flm, fm, left, half, depth + 1, st, cnt)
77 let vr: f64 = nq_node(f, m, b, fm, frm, fb, right, half, depth + 1, st, cnt)
78 return vl + vr
79}
80
81// THE RULER. See the header for the contract; the status is the answer's honesty and must be read first.
82func ni_quad_adaptive(f: func(f64) -> f64, a: f64, b: f64, tol: f64, out: *f64) -> i64 {
83 if (out as i64) == 0 { return NQ_BAD_INPUT }
84 out[NQ_OUT_VALUE] = 0.0
85 out[NQ_OUT_ERR] = 0.0
86 out[NQ_OUT_EVALS] = 0.0
87 if tol <= 0.0 { return NQ_BAD_INPUT }
88 if nq_finite(a) == 0 { return NQ_BAD_INPUT }
89 if nq_finite(b) == 0 { return NQ_BAD_INPUT }
90 if nq_finite(tol) == 0 { return NQ_BAD_INPUT }
91 if a > b { return NQ_BAD_INPUT }
92 if a == b { return NQ_OK }
93 let cnt: *i64 = sys_mmap(NQ_OUT_BYTES) as *i64
94 let st: *f64 = sys_mmap(NQ_OUT_BYTES) as *f64
95 cnt[0] = 0; cnt[1] = 0; cnt[2] = 0
96 st[3] = 0.0
97 let m: f64 = (a + b) / 2.0
98 let fa: f64 = f(a)
99 let fm: f64 = f(m)
100 let fb: f64 = f(b)
101 cnt[0] = 3
102 if nq_finite(fa) == 0 { return NQ_NONFINITE }
103 if nq_finite(fm) == 0 { return NQ_NONFINITE }
104 if nq_finite(fb) == 0 { return NQ_NONFINITE }
105 let whole: f64 = nq_simpson(a, b, fa, fm, fb)
106 let v: f64 = nq_node(f, a, b, fa, fm, fb, whole, tol, 0, st, cnt)
107 out[NQ_OUT_VALUE] = v
108 out[NQ_OUT_ERR] = st[3]
109 out[NQ_OUT_EVALS] = cnt[0] as f64
110 if cnt[1] == 1 { return NQ_NONFINITE }
111 if cnt[2] == 1 { return NQ_REFUSED_TOL }
112 if nq_finite(v) == 0 { return NQ_NONFINITE }
113 return NQ_OK
114}