nx_f32_log.nx source

1// nx_f32_log.nx -- IEEE 754 binary32 natural-log log(x) bits-up. 2// 3// L6 of the bits-up numeric tower. Composes L4 mul/add/sub/div/cvt 4// to produce real-valued log(x) without libm. 5// 6// Algorithm (clean-room from cited refs): 7// 8// 1. Special cases per IEEE 754: 9// log(NaN) = NaN 10// log(+0) = -inf 11// log(-0) = -inf 12// log(neg) = NaN 13// log(+inf) = +inf 14// log(1.0) = +0 (exact) 15// 16// 2. Range reduction: x = 2^k * m, m in [1, 2). 17// k = f32_exp_field(x) - 127 18// m = same bits with exp field set to 127 (interpret as 1.frac) 19// Then log(x) = k * ln(2) + log(m). 20// 21// 3. Polynomial for log(m) on m in [1, 2): 22// u = (m - 1) / (m + 1), u in [0, 1/3] 23// log(m) = 2 * atanh(u) 24// ~= 2 * u * (1 + u^2/3 + u^4/5 + u^6/7 + u^8/9) 25// The atanh form converges rapidly on this range (much better 26// than direct Taylor of log(1+t) which is slow near t=1). 27// 28// 4. Combine: log(x) = k * ln(2) + 2*u*P(u^2) 29// where P(u^2) = 1 + u^2/3 + u^4/5 + u^6/7 + u^8/9 (Horner). 30// 31// References absorbed clean-room (no code borrowed): 32// Hart 1968, "Computer Approximations" 33// Cody+Waite 1980, "Software Manual for Elementary Functions" 34// Muller 2016, "Elementary Functions" 35// 36// v1 accuracy: ~50-500 ULPs on f32 (atanh polynomial truncation + 37// accumulator rounding). Adequate for ML loss + entropy math; 38// v2 Remez tightening for sub-100 ULPs. 39// 40// genealogy_id: standard_log_range_reduction + atanh_form_polynomial 41// lineage_id: substrate_f32_log_v1_atanh 42 43import "nx_syscalls.nx" 44import "nx_tier.nx" 45import "nx_f32.nx" 46import "nx_f32_div.nx" 47import "nx_f32_cvt.nx" 48 49// Constants as f32 bit patterns. 50const NX_F32_LOG_ONE: i64 = 0x3F800000 // 1.0 51const NX_F32_LOG_TWO: i64 = 0x40000000 // 2.0 52const NX_F32_LOG_LN2: i64 = 0x3F317218 // ln(2) = 0.6931472 53const NX_F32_LOG_NEG_INF: i64 = 0xFF800000 // -inf 54// Polynomial coefficients (rational 1/N for atanh series). 55const NX_F32_LOG_INV_3: i64 = 0x3EAAAAAB // 1/3 56const NX_F32_LOG_INV_5: i64 = 0x3E4CCCCD // 1/5 57const NX_F32_LOG_INV_7: i64 = 0x3E124925 // 1/7 ~= 0.14285715 58const NX_F32_LOG_INV_9: i64 = 0x3DE38E39 // 1/9 ~= 0.11111111 59 60func nx_f32_log(x: i64) -> i64 { 61 let cls: nx_int = nx_f32_classify(x) 62 63 if cls == NX_F32_CLS_NAN { return NX_F32_NAN_RAW } 64 if cls == NX_F32_CLS_ZERO { return NX_F32_LOG_NEG_INF } 65 66 // Negative finite or -inf -> NaN. 67 if nx_f32_sign(x) == 1 { return NX_F32_NAN_RAW } 68 69 // +inf 70 if cls == NX_F32_CLS_INF { return 0x7F800000 } 71 72 // Exact log(1.0) = 0 73 if x == NX_F32_LOG_ONE { return 0 } 74 75 // Subnormal v1: conservative, return -inf. 76 if cls == NX_F32_CLS_SUBNORMAL { return NX_F32_LOG_NEG_INF } 77 78 // Range reduction: x = 2^k * m, m in [1, 2). 79 let exp_field: i64 = nx_f32_exp_field(x) 80 let k_int: i64 = exp_field - NX_F32_EXP_BIAS 81 let mant: i64 = nx_f32_mant_field(x) 82 // m_raw: same mantissa, exp_field = 127 (= 1.frac form) 83 let m_raw: i64 = (mant) | (NX_F32_EXP_BIAS << NX_F32_EXP_SHIFT) 84 85 // u = (m - 1) / (m + 1). 86 let m_minus_1: i64 = nx_f32_sub(m_raw, NX_F32_LOG_ONE) 87 let m_plus_1: i64 = nx_f32_add(m_raw, NX_F32_LOG_ONE) 88 let u: i64 = nx_f32_div(m_minus_1, m_plus_1) 89 let u2: i64 = nx_f32_mul(u, u) 90 91 // P(u^2) = 1 + u^2/3 + u^4/5 + u^6/7 + u^8/9 via Horner: 92 // = 1 + u^2 * (1/3 + u^2 * (1/5 + u^2 * (1/7 + u^2 * 1/9))) 93 let t1: i64 = nx_f32_mul(u2, NX_F32_LOG_INV_9) 94 let t2: i64 = nx_f32_add(NX_F32_LOG_INV_7, t1) 95 let t3: i64 = nx_f32_mul(u2, t2) 96 let t4: i64 = nx_f32_add(NX_F32_LOG_INV_5, t3) 97 let t5: i64 = nx_f32_mul(u2, t4) 98 let t6: i64 = nx_f32_add(NX_F32_LOG_INV_3, t5) 99 let t7: i64 = nx_f32_mul(u2, t6) 100 let P: i64 = nx_f32_add(NX_F32_LOG_ONE, t7) 101 102 // log(m) = 2 * u * P 103 let two_u: i64 = nx_f32_mul(NX_F32_LOG_TWO, u) 104 let log_m: i64 = nx_f32_mul(two_u, P) 105 106 // log(x) = k * ln(2) + log(m) 107 let k_f32: i64 = nx_i32_to_f32(k_int) 108 let k_ln2: i64 = nx_f32_mul(k_f32, NX_F32_LOG_LN2) 109 return nx_f32_add(k_ln2, log_m) 110}