code wiki / _hdl_build / _pe_f64erf.nx
_pe_f64erf.nx source
↩ module page · 72 lines · 3496 B
1// AUTHORED BY THE NISHI BUILDER (pattern: MATH_KERNEL / ERF_TAYLOR) -- no Claude core logic.
2// All constants from the SOVEREIGN bigfloat spec (nx_mathspec_erf; rails R1+R2 held).
3// Rung-1 domain: |x|<=1.75 series; |x|>=6 exact +-1; gap = LOUD NAN (ME2-ERF-002 erfc).
4import "nx_syscalls.nx"
5import "nx_tier.nx"
6import "nx_f64.nx"
7import "nx_f64_div.nx"
8import "nx_f64_cvt.nx"
9const ER_SMASK: i64 = 0xFFFFFFFFF8000000
10func _er_mul2(a: i64, b: i64, out: *i64) -> i64 {
11 let ah: i64 = a & ER_SMASK
12 let al: i64 = nx_f64_sub(a, ah)
13 let bh: i64 = b & ER_SMASK
14 let bl: i64 = nx_f64_sub(b, bh)
15 let p: i64 = nx_f64_mul(a, b)
16 var e: i64 = nx_f64_sub(nx_f64_mul(ah, bh), p)
17 e = nx_f64_add(e, nx_f64_mul(ah, bl))
18 e = nx_f64_add(e, nx_f64_mul(al, bh))
19 e = nx_f64_add(e, nx_f64_mul(al, bl))
20 out[0] = p
21 out[1] = e
22 return 0
23}
24func nx_f64_erf(x: i64) -> i64 {
25 let cls: nx_int = nx_f64_classify(x)
26 if cls == NX_F64_CLS_NAN { return NX_F64_NAN_RAW }
27 if cls == NX_F64_CLS_ZERO { return x }
28 let sgn: i64 = x & (1 << 63)
29 let ax: i64 = x & 0x7FFFFFFFFFFFFFFF
30 if cls == NX_F64_CLS_INF { return sgn | 4607182418800017408 }
31 if ax >= 4618441417868443648 { return sgn | 4607182418800017408 }
32 if ax > 4610560118520545280 { return NX_F64_NAN_RAW }
33 let z: i64 = nx_f64_mul(x, x)
34 var p: i64 = 0 - 5150843714325961130
35 p = nx_f64_add(nx_f64_mul(p, z), 4095049280634423949)
36 p = nx_f64_add(nx_f64_mul(p, z), 0 - 5105982346802226491)
37 p = nx_f64_add(nx_f64_mul(p, z), 4139560558244362668)
38 p = nx_f64_add(nx_f64_mul(p, z), 0 - 5061783253457760980)
39 p = nx_f64_add(nx_f64_mul(p, z), 4183183548076058075)
40 p = nx_f64_add(nx_f64_mul(p, z), 0 - 5018969547106274875)
41 p = nx_f64_add(nx_f64_mul(p, z), 4225603550606982745)
42 p = nx_f64_add(nx_f64_mul(p, z), 0 - 4976522948511461254)
43 p = nx_f64_add(nx_f64_mul(p, z), 4267132850232544578)
44 p = nx_f64_add(nx_f64_mul(p, z), 0 - 4935608437509841955)
45 p = nx_f64_add(nx_f64_mul(p, z), 4307600509085743982)
46 p = nx_f64_add(nx_f64_mul(p, z), 0 - 4895664306295453080)
47 p = nx_f64_add(nx_f64_mul(p, z), 4346949740736225411)
48 p = nx_f64_add(nx_f64_mul(p, z), 0 - 4857370668689759252)
49 p = nx_f64_add(nx_f64_mul(p, z), 4384842725451610701)
50 p = nx_f64_add(nx_f64_mul(p, z), 0 - 4820041113690619343)
51 p = nx_f64_add(nx_f64_mul(p, z), 4421362170136678432)
52 p = nx_f64_add(nx_f64_mul(p, z), 0 - 4784466991167303820)
53 p = nx_f64_add(nx_f64_mul(p, z), 4456017475183478241)
54 p = nx_f64_add(nx_f64_mul(p, z), 0 - 4750534051605340358)
55 p = nx_f64_add(nx_f64_mul(p, z), 4489013713692986168)
56 p = nx_f64_add(nx_f64_mul(p, z), 0 - 4718796608662762291)
57 p = nx_f64_add(nx_f64_mul(p, z), 4519496366895530008)
58 p = nx_f64_add(nx_f64_mul(p, z), 0 - 4689446265707761909)
59 p = nx_f64_add(nx_f64_mul(p, z), 4547511648352583708)
60 p = nx_f64_add(nx_f64_mul(p, z), 0 - 4663245410524981037)
61 p = nx_f64_add(nx_f64_mul(p, z), 4571987621712047976)
62 p = nx_f64_add(nx_f64_mul(p, z), 0 - 4640852187442739688)
63 p = nx_f64_add(nx_f64_mul(p, z), 4591870180066957722)
64 p = nx_f64_add(nx_f64_mul(p, z), 0 - 4623695617433709227)
65 p = nx_f64_add(nx_f64_mul(p, z), 4607182418800017408)
66 let sv: *i64 = sys_mmap(16) as *i64
67 _er_mul2(x, p, sv)
68 let pv: *i64 = sys_mmap(16) as *i64
69 _er_mul2(sv[0], 4607760587169110893, pv)
70 let tail: i64 = nx_f64_add(pv[1], nx_f64_mul(sv[1], 4607760587169110893))
71 return nx_f64_add(pv[0], tail)
72}