`nx_f64_quad.nx` source

1// nx_f64_quad.nx -- GAMS-H: sovereign numerical QUADRATURE (the integration class).
2// Composite Simpson's rule (O(h^4)) over a sign... over a function-id dispatch (no fn
3// pointers in the dialect); qf_eval composes our PROVEN kernels (x^2 / sin / exp / 1/x).
4//   integral_a^b f dx ~= h/3 [ f0 + 4(f1+f3+..) + 2(f2+f4+..) + fn ],  h=(b-a)/n, n even.
5// f64 software-IEEE; node x_i = a + i*h built from an exact integer-as-f64 (no drift).
6// (Gauss-Legendre / adaptive Gauss-Kronrod = the accuracy/perf upgrade -- its nodes are
7// Legendre roots, computable by composing the new nx_f64_rootfind.)
8//
9// module: nishi-core.math.f64_quad
10// depends: nishi-core.math.f64, nishi-core.math.f64_exp, nishi-core.math.f64_sincos
11// capability: F64_QUADRATURE_SIMPSON
12// license_tier: ORIGINAL
13import "nx_syscalls.nx"
14import "nx_f64.nx"
15import "nx_f64_div.nx"
16import "nx_f64_exp.nx"
17import "_pe_f64sincos.nx"
18
19const Q_ONE:   i64 = 0x3FF0000000000000    // 1.0
20const Q_TWO:   i64 = 0x4000000000000000    // 2.0
21const Q_THREE: i64 = 0x4008000000000000    // 3.0
22const Q_FOUR:  i64 = 0x4010000000000000    // 4.0
23
24func qf_eval(fid: i64, x: i64) -> i64 {
25    if fid == 0 { return nx_f64_mul(x, x) }          // x^2
26    if fid == 1 { return nx_f64_sin(x) }             // sin(x)
27    if fid == 2 { return nx_f64_exp(x) }             // e^x
28    return nx_f64_div(Q_ONE, x)                       // 1/x
29}
30
31// integer -> f64 by repeated +1.0 (exact for k < 2^53; k small here).
32func qf_i2f(k: i64) -> i64 {
33    var r: i64 = 0
34    var i: i64 = 0
35    while i < k { r = nx_f64_add(r, Q_ONE); i = i + 1 }
36    return r
37}
38
39// composite Simpson over [a,b] with n (even) panels.
40func nx_f64_quad_simpson(fid: i64, a: i64, b: i64, n: i64) -> i64 {
41    let nf: i64 = qf_i2f(n)
42    let h: i64 = nx_f64_div(nx_f64_sub(b, a), nf)
43    var sum: i64 = nx_f64_add(qf_eval(fid, a), qf_eval(fid, b))
44    var if64: i64 = Q_ONE
45    var i: i64 = 1
46    while i < n {
47        let x: i64 = nx_f64_add(a, nx_f64_mul(if64, h))
48        let fx: i64 = qf_eval(fid, x)
49        if (i & 1) == 1 { sum = nx_f64_add(sum, nx_f64_mul(Q_FOUR, fx)) }
50        else { sum = nx_f64_add(sum, nx_f64_mul(Q_TWO, fx)) }
51        if64 = nx_f64_add(if64, Q_ONE)
52        i = i + 1
53    }
54    return nx_f64_mul(sum, nx_f64_div(h, Q_THREE))
55}
nx_f64_quad.nx source

`nx_f64_quad.nx` source
