code wiki / _hdl_build / nx_f64_rootfind.nx
nx_f64_rootfind.nx source
↩ module page · 48 lines · 2169 B
1// nx_f64_rootfind.nx -- GAMS-F: sovereign nonlinear root-finder (the nonlinear-solve
2// capability class). Robust BISECTION on a sign-changing bracket [a,b] -- guaranteed
3// convergence (1 bit/iter), no derivative needed. The function under test is selected by
4// a function-id (the dialect has no function pointers); rf_eval composes our PROVEN
5// kernels (exp / cos / polynomial), so the root-finder rests only on gated primitives.
6// fid 0: x^2 - 2 (root sqrt2)
7// fid 1: cos(x) - x (Dottie number)
8// fid 2: exp(x) - 2 (root ln2)
9// fid 3: x^3 - x - 2 (root ~1.5214)
10// f64 software-IEEE throughout. (Brent superlinear convergence = the perf-exceed upgrade.)
11//
12// module: nishi-core.math.f64_rootfind
13// depends: nishi-core.math.f64, nishi-core.math.f64_exp, nishi-core.math.f64_sincos
14// capability: F64_ROOTFIND_BISECTION
15// license_tier: ORIGINAL
16import "nx_syscalls.nx"
17import "nx_f64.nx"
18import "nx_f64_exp.nx"
19import "_pe_f64sincos.nx"
20
21const RF_TWO: i64 = 0x4000000000000000 // 2.0
22const RF_HALF: i64 = 0x3FE0000000000000 // 0.5
23
24func rf_eval(fid: i64, x: i64) -> i64 {
25 if fid == 0 { return nx_f64_sub(nx_f64_mul(x, x), RF_TWO) } // x^2 - 2
26 if fid == 1 { return nx_f64_sub(nx_f64_cos(x), x) } // cos(x) - x
27 if fid == 2 { return nx_f64_sub(nx_f64_exp(x), RF_TWO) } // exp(x) - 2
28 // x^3 - x - 2
29 let x3: i64 = nx_f64_mul(nx_f64_mul(x, x), x)
30 return nx_f64_sub(nx_f64_sub(x3, x), RF_TWO)
31}
32
33// bisection root of f_fid in [a,b] (requires sign change); 70 halvings -> below f64 ulp.
34func nx_f64_root_bisect(fid: i64, a: i64, b: i64) -> i64 {
35 var lo: i64 = a
36 var hi: i64 = b
37 var flo: i64 = rf_eval(fid, lo)
38 var i: i64 = 0
39 while i < 70 {
40 let m: i64 = nx_f64_mul(nx_f64_add(lo, hi), RF_HALF)
41 let fm: i64 = rf_eval(fid, m)
42 // same sign as flo? sign bits equal -> root is in [m, hi]; else [lo, m]
43 let diff: i64 = ((flo ^ fm) >> 63) & 1
44 if diff == 0 { lo = m; flo = fm } else { hi = m }
45 i = i + 1
46 }
47 return nx_f64_mul(nx_f64_add(lo, hi), RF_HALF)
48}