nx_pattern_emit6_gamma.nx source

1// nx_pattern_emit6_gamma.nx -- PATTERN EMITTER: MATH_KERNEL / GAMMA_STIRLING (DLMF ch5, 2// ME2-GAMMA-DD4). Transcribes the VALIDATED gd_gamma (dd lnGamma, max_ulp=2 vs bf 3// oracle) into an emitted f64 kernel -- the novel design (probe->dd-verify->root-cause) 4// is DONE; this is the mechanical emission. Composes the authored log + sinhcosh 5// kernels (kernel composition: the emitted module imports them). Every constant from 6// the sovereign spec table c[] (no magic numbers). Rung-1 domain x>=0.5 (x<0.5 7// reflection = LOUD NAN, ME2-GAMMA-reflection debt -- the erf gap-contract gene). 8// license_tier: ORIGINAL 9import "nx_pattern_emit.nx" 10import "nx_syscalls.nx" 11func pe6g_wlit(fd: i64, v: i64) -> i64 { 12 if v < 0 { pe_w(fd, "0 - " as *u8); pe_wn(fd, 0 - v); return 0 } 13 pe_wn(fd, v) 14 return 0 15} 16func pe6_author_gamma(modfile: *u8, modpath: *u8, testpath: *u8, c: *i64) -> i64 { 17 let mf: i64 = sys_openat_wr(modpath, 0x1a4); if mf < 0 { return 0 } 18 pe_w(mf, "// AUTHORED BY THE NISHI BUILDER (pattern: MATH_KERNEL / GAMMA_STIRLING) -- no Claude core logic.\n" as *u8) 19 pe_w(mf, "// dd lnGamma (Stirling + pull-up), composed from authored log+sinhcosh kernels; spec=sovereign bigfloat.\n" as *u8) 20 pe_w(mf, "// Validated max_ulp=2 vs bf oracle. Rung-1 x>=0.5; x<0.5 reflection = LOUD NAN (named debt).\n" as *u8) 21 pe_w(mf, "import \"nx_syscalls.nx\"\nimport \"nx_tier.nx\"\nimport \"nx_f64.nx\"\nimport \"nx_f64_div.nx\"\nimport \"nx_f64_cvt.nx\"\n" as *u8) 22 pe_w(mf, "import \"_pe_f64log.nx\"\nimport \"_pe_f64sinhcosh.nx\"\n" as *u8) 23 // exp via proven cosh+sinh; ln_dd via Newton correction 24 pe_w(mf, "func _gm_expf(v: i64) -> i64 { return nx_f64_add(nx_f64_cosh(v), nx_f64_sinh(v)) }\n" as *u8) 25 pe_w(mf, "func _gm_lndd(w: i64, out: *i64) -> i64 {\n" as *u8) 26 pe_w(mf, " let hi: i64 = nx_f64_log(w)\n" as *u8) 27 pe_w(mf, " let eh: i64 = _gm_expf(hi)\n" as *u8) 28 pe_w(mf, " out[0] = hi\n out[1] = nx_f64_div(nx_f64_sub(w, eh), w)\n return 0\n}\n" as *u8) 29 // the kernel 30 pe_w(mf, "func nx_f64_gamma(x: i64) -> i64 {\n" as *u8) 31 pe_w(mf, " let cls: nx_int = nx_f64_classify(x)\n" as *u8) 32 pe_w(mf, " if cls == NX_F64_CLS_NAN { return NX_F64_NAN_RAW }\n" as *u8) 33 pe_w(mf, " if nx_f64_sign(x) == 1 { return NX_F64_NAN_RAW }\n" as *u8) 34 pe_w(mf, " if cls == NX_F64_CLS_INF { return x }\n" as *u8) 35 // x < 0.5 (lo_cut c[13]) -> reflection debt = LOUD NAN 36 pe_w(mf, " if (x & 0x7FFFFFFFFFFFFFFF) < " as *u8); pe6g_wlit(mf, c[13]); pe_w(mf, " { return NX_F64_NAN_RAW }\n" as *u8) 37 // overflow cut c[4] 38 pe_w(mf, " if nx_f64_sub(x, " as *u8); pe6g_wlit(mf, c[4]); pe_w(mf, ") >> 63 == 0 { return 0x7FF0000000000000 }\n" as *u8) 39 // pull-up to w>=8 (c[3]); prod = x*(x+1)*... 40 pe_w(mf, " var w: i64 = x\n var prod: i64 = " as *u8); pe6g_wlit(mf, c[0]); pe_w(mf, "\n var guard: i64 = 0\n" as *u8) 41 pe_w(mf, " while guard < 64 {\n" as *u8) 42 pe_w(mf, " if nx_f64_sub(w, " as *u8); pe6g_wlit(mf, c[3]); pe_w(mf, ") >> 63 == 0 { guard = 64 } else {\n" as *u8) 43 pe_w(mf, " prod = nx_f64_mul(prod, w)\n w = nx_f64_add(w, " as *u8); pe6g_wlit(mf, c[0]); pe_w(mf, ")\n guard = guard + 1\n }\n }\n" as *u8) 44 pe_w(mf, " let ln: *i64 = sys_mmap(32) as *i64\n _gm_lndd(w, ln)\n" as *u8) 45 pe_w(mf, " let wm: i64 = nx_f64_sub(w, " as *u8); pe6g_wlit(mf, c[1]); pe_w(mf, ")\n" as *u8) 46 pe_w(mf, " let mp: *i64 = sys_mmap(32) as *i64\n _sh_mul2(wm, ln[0], mp)\n" as *u8) 47 pe_w(mf, " var a_hi: i64 = mp[0]\n var a_lo: i64 = nx_f64_add(mp[1], nx_f64_mul(wm, ln[1]))\n" as *u8) 48 pe_w(mf, " let acc: *i64 = sys_mmap(32) as *i64\n" as *u8) 49 pe_w(mf, " _sh_2sum(a_hi, nx_f64_neg(w), acc)\n var g_hi: i64 = acc[0]\n var g_lo: i64 = nx_f64_add(acc[1], a_lo)\n" as *u8) 50 pe_w(mf, " _sh_2sum(g_hi, " as *u8); pe6g_wlit(mf, c[2]); pe_w(mf, ", acc)\n g_hi = acc[0]\n g_lo = nx_f64_add(g_lo, acc[1])\n" as *u8) 51 pe_w(mf, " let winv: i64 = nx_f64_div(" as *u8); pe6g_wlit(mf, c[0]); pe_w(mf, ", w)\n let w2inv: i64 = nx_f64_mul(winv, winv)\n" as *u8) 52 pe_w(mf, " var term: i64 = winv\n var ser: i64 = nx_f64_mul(" as *u8); pe6g_wlit(mf, c[12]); pe_w(mf, ", term)\n" as *u8) 53 var k: i64 = 2 54 while k <= 8 { 55 pe_w(mf, " term = nx_f64_mul(term, w2inv)\n ser = nx_f64_add(ser, nx_f64_mul(" as *u8) 56 pe6g_wlit(mf, c[12 - (k - 1)]) 57 pe_w(mf, ", term))\n" as *u8) 58 k = k + 1 59 } 60 // series into g_hi (the validated fix: significant lnGamma term -> exponent) 61 pe_w(mf, " _sh_2sum(g_hi, ser, acc)\n g_hi = acc[0]\n g_lo = nx_f64_add(g_lo, acc[1])\n" as *u8) 62 pe_w(mf, " let eg: i64 = _gm_expf(g_hi)\n let gw: i64 = nx_f64_add(eg, nx_f64_mul(eg, g_lo))\n" as *u8) 63 pe_w(mf, " return nx_f64_div(gw, prod)\n}\n" as *u8) 64 sys_close(mf) 65 // emitted invariant test: Gamma(integer) = (integer-1)! exactly-ish + specials + monotone-ish 66 let tf: i64 = sys_openat_wr(testpath, 0x1a4); if tf < 0 { return 0 } 67 pe_w(tf, "// GENERATED invariant KATs for the Builder-authored gamma (GAMMA_STIRLING).\n" as *u8) 68 pe_w(tf, "import \"" as *u8); pe_w(tf, modfile); pe_w(tf, "\"\n" as *u8) 69 pe_w(tf, "func gt_ulp(a: i64, b: i64) -> i64 { if a >= b { return a - b } return b - a }\n" as *u8) 70 pe_w(tf, "func main() -> i64 {\n var bad: i64 = 0\n" as *u8) 71 pe_w(tf, " if nx_f64_gamma(0x7FF8000000000000) != NX_F64_NAN_RAW { bad = bad + 1 }\n" as *u8) // NaN 72 pe_w(tf, " if nx_f64_gamma(0xBFF0000000000000) != NX_F64_NAN_RAW { bad = bad + 1 }\n" as *u8) // -1 -> NaN 73 pe_w(tf, " if nx_f64_gamma(0x3FD0000000000000) != NX_F64_NAN_RAW { bad = bad + 1 }\n" as *u8) // 0.25 < 0.5 -> NaN 74 // Gamma(2)=1, Gamma(3)=2, Gamma(4)=6, Gamma(5)=24, Gamma(6)=120 (<=2 ulp) 75 pe_w(tf, " if gt_ulp(nx_f64_gamma(0x4000000000000000), 0x3FF0000000000000) > 2 { bad = bad + 1 }\n" as *u8) // G(2)=1 76 pe_w(tf, " if gt_ulp(nx_f64_gamma(0x4008000000000000), 0x4000000000000000) > 2 { bad = bad + 1 }\n" as *u8) // G(3)=2 77 pe_w(tf, " if gt_ulp(nx_f64_gamma(0x4010000000000000), 0x4018000000000000) > 2 { bad = bad + 1 }\n" as *u8) // G(4)=6 78 pe_w(tf, " if gt_ulp(nx_f64_gamma(0x4014000000000000), 0x4038000000000000) > 2 { bad = bad + 1 }\n" as *u8) // G(5)=24 79 pe_w(tf, " if gt_ulp(nx_f64_gamma(0x4018000000000000), 0x405E000000000000) > 2 { bad = bad + 1 }\n" as *u8) // G(6)=120 80 pe_w(tf, " sys_exit(bad)\n return bad\n}\n" as *u8) 81 sys_close(tf) 82 return 1 83}