nx_f64_exp.nx source
↩ module page · 78 lines · 3296 B
1// nx_f64_exp.nx -- binary64 exponential, bits-up (ME1 rung 1, DLMF ch.4).
2//
3// THE TEMPLATE for the ME1 transcendental family (and the future Builder
4// MATH_KERNEL pattern emitter): range-reduce -> polynomial -> rescale, with
5// every constant GENERATED from a high-precision runtime oracle
6// (_f64_exp_consts.nx, mpmath @60dps) and a ULP-bounded KAT gate
7// (_f64_exp_gate_authored, vectors correctly rounded by the same oracle).
8//
9// Algorithm (Cody & Waite 1980 range reduction; Hart 1968 poly approach):
10// 1. specials: NaN->NaN, +inf->+inf, -inf->0, +-0->1; |x| beyond the
11// over/underflow cutoffs short-circuits to inf / 0.
12// 2. k = nearest-int(x / ln2) (trunc of x*INV_LN2 +- 0.5)
13// 3. r = (x - k*LN2_HI) - k*LN2_LO -- k*LN2_HI exact (25-bit hi split,
14// |k| <= 1075 = 11 bits), so r carries ~full f64 precision; |r| <= ln2/2.
15// 4. exp(r) = Horner sum of r^i/i!, i = 0..15 (truncation < 1e-21 rel,
16// far below half-ulp; observed error is rounding accumulation only).
17// 5. result = ldexp(exp(r), k) -- correct gradual underflow via round_pack.
18//
19// ULP bound: GATED, not asserted -- see math_engine.log EXP-ULP line.
20// One transcendental per file (family law; see nx_f32_div.nx quirk note).
21//
22// license_tier: ORIGINAL
23
24import "nx_syscalls.nx"
25import "nx_tier.nx"
26import "nx_f64.nx"
27import "nx_f64_cvt.nx"
28import "_f64_exp_consts.nx"
29
30func nx_f64_exp(x: i64) -> i64 {
31 let cls: nx_int = nx_f64_classify(x)
32 if cls == NX_F64_CLS_NAN { return NX_F64_NAN_RAW }
33 if cls == NX_F64_CLS_ZERO { return EXPC_ONE }
34 if cls == NX_F64_CLS_INF {
35 if nx_f64_sign(x) == 1 { return 0 }
36 return NX_F64_INF_RAW
37 }
38 // Overflow / underflow cutoffs (exact thresholds handled by ldexp; these
39 // only guard the f64->i64 conversion from huge k).
40 if nx_f64_lt(EXPC_HI_CUT, x) == 1 { return NX_F64_INF_RAW }
41 if nx_f64_lt(x, EXPC_LO_CUT) == 1 { return 0 }
42
43 // k = nearest integer to x/ln2 (half-away-from-zero; a half-ulp slop in
44 // the tie direction only moves r within the poly's convergence margin).
45 let t: i64 = nx_f64_mul(x, EXPC_INV_LN2)
46 var tk: i64 = 0
47 if nx_f64_sign(t) == 1 {
48 tk = nx_f64_sub(t, EXPC_HALF)
49 } else {
50 tk = nx_f64_add(t, EXPC_HALF)
51 }
52 let k: i64 = nx_f64_to_i64(tk)
53 let kf: i64 = nx_i64_to_f64(k)
54
55 // r = (x - k*ln2_hi) - k*ln2_lo
56 var r: i64 = nx_f64_sub(x, nx_f64_mul(kf, EXPC_LN2_HI))
57 r = nx_f64_sub(r, nx_f64_mul(kf, EXPC_LN2_LO))
58
59 // Horner: p = C15; p = p*r + C(i) down to C0 = 1.
60 var p: i64 = EXPC_C15
61 p = nx_f64_add(nx_f64_mul(p, r), EXPC_C14)
62 p = nx_f64_add(nx_f64_mul(p, r), EXPC_C13)
63 p = nx_f64_add(nx_f64_mul(p, r), EXPC_C12)
64 p = nx_f64_add(nx_f64_mul(p, r), EXPC_C11)
65 p = nx_f64_add(nx_f64_mul(p, r), EXPC_C10)
66 p = nx_f64_add(nx_f64_mul(p, r), EXPC_C9)
67 p = nx_f64_add(nx_f64_mul(p, r), EXPC_C8)
68 p = nx_f64_add(nx_f64_mul(p, r), EXPC_C7)
69 p = nx_f64_add(nx_f64_mul(p, r), EXPC_C6)
70 p = nx_f64_add(nx_f64_mul(p, r), EXPC_C5)
71 p = nx_f64_add(nx_f64_mul(p, r), EXPC_C4)
72 p = nx_f64_add(nx_f64_mul(p, r), EXPC_C3)
73 p = nx_f64_add(nx_f64_mul(p, r), EXPC_C2)
74 p = nx_f64_add(nx_f64_mul(p, r), EXPC_C1)
75 p = nx_f64_add(nx_f64_mul(p, r), EXPC_C0)
76
77 return nx_f64_ldexp(p, k)
78}