nx_bigfloat120_sinh.nx source
↩ module page · 113 lines · 4000 B
1// nx_bigfloat120_sinh.nx -- SOVEREIGN sinh/cosh oracle on the 120-bit bigfloat.
2// cosh = (E + 1/E)/2 with E = exp(x) from the pow organ's bigfloat exp core
3// (bfp_exp_pm) -- no cancellation anywhere. sinh uses the same form for
4// x >= 1 (cancellation there costs < 1 bit) and the EXACT-RATIONAL Taylor
5// series x + x^3/3! + ... for x < 1, where the exp difference would cancel
6// catastrophically at 120 bits. Positive-only core; sinh odd / cosh even
7// handled at the boundary. Overflow: |x| > 711 decided as inf (true cut
8// ~710.476 lands inside bf_to_f64's natural overflow rounding).
9//
10// Self-anchor (see _bf_sinhcosh_gate_authored): double-angle identities
11// sinh(2u) = 2 sinh u cosh u and cosh(2u) = 1 + 2 sinh^2 u (cancellation-free,
12// valid across the whole domain), the pythagorean strip at small |x|, and an
13// f64 tie sinh+cosh = e^x against the BLESSED bf_exp_f64 organ.
14// license_tier: ORIGINAL
15
16import "nx_syscalls.nx"
17import "nx_tier.nx"
18import "nx_bigfloat120.nx"
19import "nx_bigfloat120_div.nx"
20import "nx_bigfloat120_exp.nx"
21import "nx_bigfloat120_ln.nx"
22import "nx_bigfloat120_pow.nx"
23const K_MAGIC_2047: i64 = 2047
24
25// sinh Taylor for bigfloat 0 <= x < 1: x * (1 + sum x^(2k)/(2k+1)!), 22 terms
26func _bsh_taylor_sinh(out: *i64, x: *i64) -> i64 {
27 if bf_is_zero(x) == 1 { return bf_copy(out, x) }
28 let z: *i64 = bf_new()
29 bf_mul(z, x, x)
30 let term: *i64 = bf_new()
31 bf_copy(term, x)
32 let sum: *i64 = bf_new()
33 bf_copy(sum, x)
34 let tmul: *i64 = bf_new()
35 let tdiv: *i64 = bf_new()
36 let snew: *i64 = bf_new()
37 var k: i64 = 1
38 while k <= 22 {
39 bf_mul(tmul, term, z)
40 bf_div_small(tdiv, tmul, (2 * k) * (2 * k + 1))
41 bf_copy(term, tdiv)
42 bf_add(snew, sum, term)
43 bf_copy(sum, snew)
44 k = k + 1
45 }
46 return bf_copy(out, sum)
47}
48
49// sinh(xb) for positive bigfloat xb (caller bounds xb <= ~711)
50func bf_sinh_xb(out: *i64, xb: *i64) -> i64 {
51 if bf_is_zero(xb) == 1 { return bf_copy(out, xb) }
52 let one: *i64 = bf_new()
53 bf_set_int(one, 1)
54 if bf_cmp(xb, one) < 0 { return _bsh_taylor_sinh(out, xb) }
55 let e: *i64 = bf_new()
56 bfp_exp_pm(e, xb, 0)
57 let ei: *i64 = bf_new()
58 bfp_exp_pm(ei, xb, 1)
59 bf_sub(out, e, ei) // e > 1/e for x >= 1
60 out[0] = out[0] - 1 // / 2
61 return 0
62}
63
64// cosh(xb) for positive-or-zero bigfloat xb (caller bounds xb <= ~711)
65func bf_cosh_xb(out: *i64, xb: *i64) -> i64 {
66 let one: *i64 = bf_new()
67 bf_set_int(one, 1)
68 if bf_is_zero(xb) == 1 { return bf_copy(out, one) }
69 let e: *i64 = bf_new()
70 bfp_exp_pm(e, xb, 0)
71 let ei: *i64 = bf_new()
72 bfp_exp_pm(ei, xb, 1)
73 bf_add(out, e, ei)
74 out[0] = out[0] - 1 // / 2
75 return 0
76}
77
78func bf_sinh_f64(x: i64) -> i64 {
79 let ax: i64 = x & 0x7FFFFFFFFFFFFFFF
80 let ef: i64 = (x >> 52) & 0x7FF
81 let sgn: i64 = (x >> 63) & 1
82 if ef == K_MAGIC_2047 {
83 if (ax & 0x000FFFFFFFFFFFFF) != 0 { return 0x7FF8000000000000 }
84 return x // sinh(+-inf) = +-inf
85 }
86 if ax == 0 { return x }
87 let xb: *i64 = bf_new()
88 bf_set_f64(xb, ax)
89 let cut: *i64 = bf_new()
90 bf_set_int(cut, 711)
91 if bf_cmp(xb, cut) > 0 { return (sgn << 63) | 0x7FF0000000000000 }
92 let m: *i64 = bf_new()
93 bf_sinh_xb(m, xb)
94 return bf_to_f64(m, sgn, 0)
95}
96
97func bf_cosh_f64(x: i64) -> i64 {
98 let ax: i64 = x & 0x7FFFFFFFFFFFFFFF
99 let ef: i64 = (x >> 52) & 0x7FF
100 if ef == K_MAGIC_2047 {
101 if (ax & 0x000FFFFFFFFFFFFF) != 0 { return 0x7FF8000000000000 }
102 return 0x7FF0000000000000 // cosh(+-inf) = +inf
103 }
104 if ax == 0 { return 0x3FF0000000000000 } // cosh(0) = 1
105 let xb: *i64 = bf_new()
106 bf_set_f64(xb, ax)
107 let cut: *i64 = bf_new()
108 bf_set_int(cut, 711)
109 if bf_cmp(xb, cut) > 0 { return 0x7FF0000000000000 }
110 let m: *i64 = bf_new()
111 bf_cosh_xb(m, xb)
112 return bf_to_f64(m, 0, 0)
113}