nx_f32_sincos.nx source

1// nx_f32_sincos.nx -- IEEE 754 binary32 sin(x) + cos(x) bits-up. 2// 3// L6 transcendentals. Composes L4 mul/add/sub/cvt. No libm. 4// 5// Algorithm (clean-room from cited refs): 6// 7// 1. Range reduction: 8// k = round(x * 2/pi) 9// r = x - k * pi/2, r in [-pi/4, pi/4] 10// Quadrant q = k mod 4. 11// 12// 2. Reduced-range polynomials (Taylor, 4 odd / 4 even terms): 13// sin(r) ~= r - r^3/6 + r^5/120 - r^7/5040 14// ~= r * (1 - r^2/6 + r^4/120 - r^6/5040) 15// cos(r) ~= 1 - r^2/2 + r^4/24 - r^6/720 16// 17// 3. Quadrant dispatch: 18// sin: q=0 -> sin(r) ; q=1 -> cos(r) ; q=2 -> -sin(r) ; q=3 -> -cos(r) 19// cos: q=0 -> cos(r) ; q=1 -> -sin(r) ; q=2 -> -cos(r) ; q=3 -> sin(r) 20// 21// References absorbed clean-room: 22// Hart 1968, Cody+Waite 1980, Muller 2016 23// Standard quadrant-based range-reduction canon 24// 25// v1 accuracy: ~1000-10000 ULPs (Taylor truncation + range reduction 26// loss; classic single-word pi/2 issue). Adequate for ML positional 27// encoding (RoPE). v2 Remez minimax + Cody-Waite split pi/2 for 28// tighter bounds queued. 29// 30// genealogy_id: standard_quadrant_reduction + taylor_horner_canon 31// lineage_id: substrate_f32_sincos_v1_taylor4 32 33import "nx_syscalls.nx" 34import "nx_tier.nx" 35import "nx_f32.nx" 36import "nx_f32_cvt.nx" 37const NX_MAGIC_4294967295: i64 = 4294967295 38 39// f32 constants. 40const NX_F32_SC_ONE: i64 = 0x3F800000 // 1.0 41const NX_F32_SC_PI_2: i64 = 0x3FC90FDB // pi/2 = 1.5707963 42const NX_F32_SC_2_PI: i64 = 0x3F22F983 // 2/pi = 0.6366198 43const NX_F32_SC_INV_6: i64 = 0x3E2AAAAB // 1/6 44const NX_F32_SC_INV_120: i64 = 0x3C088889 // 1/120 45const NX_F32_SC_INV_5040: i64 = 0x39500D01 // 1/5040 46const NX_F32_SC_INV_2: i64 = 0x3F000000 // 1/2 47const NX_F32_SC_INV_24: i64 = 0x3D2AAAAB // 1/24 48const NX_F32_SC_INV_720: i64 = 0x3AB60B61 // 1/720 49 50// f32 -> i32 round-to-nearest-even (private helper; duplicated from 51// nx_f32_exp.nx to keep this brick standalone per Task #10 discipline). 52 53func _sc_f32_to_i32_rne(value: i64) -> i64 { 54 let cls: nx_int = nx_f32_classify(value) 55 if cls == NX_F32_CLS_ZERO { return 0 } 56 if cls == NX_F32_CLS_NAN { return 0 } 57 if cls == NX_F32_CLS_INF { return 0 } 58 59 let sign: i64 = nx_f32_sign(value) 60 let exp_field: i64 = nx_f32_exp_field(value) 61 let mant: i64 = nx_f32_mant_field(value) 62 let sig: i64 = mant | NX_F32_IMPLICIT_1 63 let real_e: i64 = exp_field - NX_F32_EXP_BIAS 64 65 var result: i64 = 0 66 if real_e < 0 - 1 { return 0 } 67 if real_e >= 23 { 68 let shift_up: i64 = real_e - 23 69 if shift_up > 30 { return 0 } 70 result = sig << shift_up 71 } else { 72 let shift_down: i64 = 23 - real_e 73 let lost_mask: i64 = (1 << shift_down) - 1 74 let lost: i64 = sig & lost_mask 75 result = sig >> shift_down 76 let halfway: i64 = 1 << (shift_down - 1) 77 if lost > halfway { result = result + 1 } 78 if lost == halfway { 79 if (result & 1) == 1 { result = result + 1 } 80 } 81 } 82 if sign == 1 { return 0 - result } 83 return result 84} 85 86// sin polynomial on reduced range r in [-pi/4, pi/4]: 87// sin(r) ~= r * (1 - r^2/6 + r^4/120 - r^6/5040) 88// = r * (1 + r^2 * (-1/6 + r^2 * (1/120 + r^2 * (-1/5040)))) 89// (Horner form; signs absorbed into coefficient signs.) 90 91func _sc_sin_poly(r: i64) -> i64 { 92 let r2: i64 = nx_f32_mul(r, r) 93 // Build Horner: 94 // inner3 = -1/5040 95 // inner2 = 1/120 + r^2 * inner3 96 // inner1 = -1/6 + r^2 * inner2 97 // shell = 1 + r^2 * inner1 98 // result = r * shell 99 let neg_inv_5040: i64 = nx_f32_neg(NX_F32_SC_INV_5040) 100 let neg_inv_6: i64 = nx_f32_neg(NX_F32_SC_INV_6) 101 let t1: i64 = nx_f32_mul(r2, neg_inv_5040) 102 let t2: i64 = nx_f32_add(NX_F32_SC_INV_120, t1) 103 let t3: i64 = nx_f32_mul(r2, t2) 104 let t4: i64 = nx_f32_add(neg_inv_6, t3) 105 let t5: i64 = nx_f32_mul(r2, t4) 106 let shell: i64 = nx_f32_add(NX_F32_SC_ONE, t5) 107 return nx_f32_mul(r, shell) 108} 109 110// cos polynomial on reduced range r in [-pi/4, pi/4]: 111// cos(r) ~= 1 - r^2/2 + r^4/24 - r^6/720 112// = 1 + r^2 * (-1/2 + r^2 * (1/24 + r^2 * (-1/720))) 113 114func _sc_cos_poly(r: i64) -> i64 { 115 let r2: i64 = nx_f32_mul(r, r) 116 let neg_inv_720: i64 = nx_f32_neg(NX_F32_SC_INV_720) 117 let neg_inv_2: i64 = nx_f32_neg(NX_F32_SC_INV_2) 118 let t1: i64 = nx_f32_mul(r2, neg_inv_720) 119 let t2: i64 = nx_f32_add(NX_F32_SC_INV_24, t1) 120 let t3: i64 = nx_f32_mul(r2, t2) 121 let t4: i64 = nx_f32_add(neg_inv_2, t3) 122 let t5: i64 = nx_f32_mul(r2, t4) 123 return nx_f32_add(NX_F32_SC_ONE, t5) 124} 125 126// ===== Range reduction shared by sin and cos ======================= 127// 128// Returns (k_quad, r) where x = k*pi/2 + r. 129// k_quad is k mod 4 (0..3), r in [-pi/4, pi/4]. 130// Caller passes out_r and gets k_quad via return value. 131 132func _sc_range_reduce(x: i64, out_r: *i64) -> i64 { 133 let y: i64 = nx_f32_mul(x, NX_F32_SC_2_PI) 134 let k: i64 = _sc_f32_to_i32_rne(y) 135 let kf: i64 = nx_i32_to_f32(k) 136 let r: i64 = nx_f32_sub(x, nx_f32_mul(kf, NX_F32_SC_PI_2)) 137 out_r[0] = r 138 let q: i64 = k & 3 // k mod 4 (sign-safe for our bounded usage) 139 if q < 0 { return q + 4 } 140 return q 141} 142 143// PACKED range reduce: (q << 32) | (r_raw & 0xFFFFFFFF). An f32 raw fits in 32 bits, so ONE return carries 144// both -- no out-param slot, NO sys_mmap. The out-param variant above cost one mmap PER TRIG CALL, which in a 145// training loop = millions of leaked 4KB pages -> OOM (measured: the R3e f32 reader died run-exit=137 from 146// per-element RoPE sin/cos before the fix). Debt eaten 2026-07-10: sin/cos now allocation-free. 147func _sc_range_reduce_pk(x: i64) -> i64 { 148 let y: i64 = nx_f32_mul(x, NX_F32_SC_2_PI) 149 let k: i64 = _sc_f32_to_i32_rne(y) 150 let kf: i64 = nx_i32_to_f32(k) 151 let r: i64 = nx_f32_sub(x, nx_f32_mul(kf, NX_F32_SC_PI_2)) 152 var q: i64 = k & 3 153 if q < 0 { q = q + 4 } 154 return (q << 32) | (r & NX_MAGIC_4294967295) 155} 156 157// ===== sin(x) ====================================================== 158 159func nx_f32_sin(x: i64) -> i64 { 160 let cls: nx_int = nx_f32_classify(x) 161 if cls == NX_F32_CLS_NAN { return NX_F32_NAN_RAW } 162 if cls == NX_F32_CLS_ZERO { return x } // sin(+-0) = +-0 163 if cls == NX_F32_CLS_INF { return NX_F32_NAN_RAW } 164 165 let pk: i64 = _sc_range_reduce_pk(x) // allocation-free (was mmap-per-call) 166 let q: i64 = pk >> 32 167 let r: i64 = pk & NX_MAGIC_4294967295 168 169 if q == 0 { return _sc_sin_poly(r) } 170 if q == 1 { return _sc_cos_poly(r) } 171 if q == 2 { return nx_f32_neg(_sc_sin_poly(r)) } 172 return nx_f32_neg(_sc_cos_poly(r)) 173} 174 175// ===== cos(x) ====================================================== 176 177func nx_f32_cos(x: i64) -> i64 { 178 let cls: nx_int = nx_f32_classify(x) 179 if cls == NX_F32_CLS_NAN { return NX_F32_NAN_RAW } 180 if cls == NX_F32_CLS_ZERO { return NX_F32_SC_ONE } // cos(+-0) = 1 181 if cls == NX_F32_CLS_INF { return NX_F32_NAN_RAW } 182 183 let pk: i64 = _sc_range_reduce_pk(x) // allocation-free (was mmap-per-call) 184 let q: i64 = pk >> 32 185 let r: i64 = pk & NX_MAGIC_4294967295 186 187 if q == 0 { return _sc_cos_poly(r) } 188 if q == 1 { return nx_f32_neg(_sc_sin_poly(r)) } 189 if q == 2 { return nx_f32_neg(_sc_cos_poly(r)) } 190 return _sc_sin_poly(r) 191}