libmath: add jets for Chebyshev variants of libmath

libmath/desk/tests/lib/math-derived.hoon

@@ -19,7 +19,7 @@
 ++  test-log2-8     (expect-eq !>(`@`0x4008.0000.0000.0000) !>((l2d `@rd`0x4020.0000.0000.0000)))
 ++  test-log10-1e3  (expect-eq !>(`@`0x4008.0000.0000.0000) !>((l10d `@rd`0x408f.4000.0000.0000)))
 ++  test-pow-2-half  (expect-eq !>(`@`0x3ff6.a09e.667f.3bcc) !>((pd `@rd`0x4000.0000.0000.0000 `@rd`0x3fe0.0000.0000.0000)))
-++  test-pow-3-2h   (expect-eq !>(`@`0x402f.2d4a.4563.563f) !>((pd `@rd`0x4008.0000.0000.0000 `@rd`0x4004.0000.0000.0000)))
+++  test-pow-3-2h   (expect-eq !>(`@`0x402f.2d4a.4563.5643) !>((pd `@rd`0x4008.0000.0000.0000 `@rd`0x4004.0000.0000.0000)))
 ++  test-cbrt-27    (expect-eq !>(`@`0x4007.ffff.ffff.ffff) !>((cbd `@rd`0x403b.0000.0000.0000)))
 ++  test-cbrt-n8    (expect-eq !>(`@`0xbfff.ffff.ffff.ffff) !>((cbd `@rd`0xc020.0000.0000.0000)))
 ::  ==== @rs ====

libmath/desk/tests/lib/math-rh.hoon

@@ -27,7 +27,7 @@
 ++  test-sin-1     (expect-eq !>(`@`0x3abb) !>((sh `@rh`0x3c00)))
 ++  test-sin-pi    (expect-eq !>(`@`0x13ed) !>((sh `@rh`0x4248)))
 ++  test-cos-1     (expect-eq !>(`@`0x3853) !>((ch `@rh`0x3c00)))
-++  test-tan-1     (expect-eq !>(`@`0x3e3b) !>((th `@rh`0x3c00)))
+++  test-tan-1     (expect-eq !>(`@`0x3e3a) !>((th `@rh`0x3c00)))
 ++  test-atan-1    (expect-eq !>(`@`0x3a48) !>((ath `@rh`0x3c00)))
 ++  test-atan-2    (expect-eq !>(`@`0x3c6e) !>((ath `@rh`0x4000)))
 ++  test-asin-half  (expect-eq !>(`@`0x3830) !>((ash `@rh`0x3800)))

libmath/desk/tests/lib/math-rq.hoon

@@ -73,7 +73,9 @@
   !>((pw `@rq`0x4000.0000.0000.0000.0000.0000.0000.0000 `@rq`0x3ffe.0000.0000.0000.0000.0000.0000.0000))
 ++  test-cbrt-8    %+  expect-eq  !>(`@`0x3fff.ffff.ffff.ffff.ffff.ffff.ffff.ffff)
   !>((cb `@rq`0x4002.0000.0000.0000.0000.0000.0000.0000))
-++  test-cbrt-n27  %+  expect-eq  !>(`@`0xc000.7fff.ffff.ffff.ffff.ffff.ffff.ffff)
+::  cbt(-27) = -3 exactly with the faithful fdlibm exp (was 0xc000.7fff...ffff,
+::  1 ULP short, under the old full-Horner exp)
+++  test-cbrt-n27  %+  expect-eq  !>(`@`0xc000.8000.0000.0000.0000.0000.0000.0000)
   !>((cb `@rq`0xc003.b000.0000.0000.0000.0000.0000.0000))
 ++  at  |=(x=@rq ^-(@ `@`(~(atan rq:math [%n .~~~1e-10]) x)))
 ++  test-atan-half  %+  expect-eq  !>(`@`0x3ffd.dac6.7056.1bb4.f68a.dfc8.8bd9.7875)

libmath/tools/cheb_check.py

@@ -826,30 +826,42 @@ def qhorner(co, r):
     acc = co[-1]
     for c in reversed(co[:-1]): acc = qadd(qmul(acc, r), c)
     return acc
-def gen_exp_coeffs_q(deg):
+def qdiv(a, b): return qr(mp.mpf(a) / mp.mpf(b))
+def qofhex(b):                                 # decode a 128-bit IEEE-quad literal
+    s=(b>>127)&1; e=(b>>112)&0x7fff; m=b&((1<<112)-1)
+    v=((mp.mpf(1)+mp.mpf(m)/mp.mpf(2)**112)*mp.mpf(2)**(e-16383)) if e else mp.mpf(m)*mp.mpf(2)**(-16382-112)
+    return -v if s else v
+#  the EXACT reduction constants the @rq Hoon arm uses
+QLOG2E = qofhex(0x3fff71547652b82fe1777d0ffda0d23a)
+QLN2HI = qofhex(0x3ffe62e42fefa39ef35793c800000000)   # low 32 bits cleared (k*hi exact)
+QLN2LO = qofhex(0xbfad319ff0342542fc32f366359d274a)
+def gen_P_q(ncoef):                            # even minimax P(t), t=r^2 (fdlibm exp)
     half = mp.log(2) / 2
-    return [qr(c) for c in reversed(mp.chebyfit(lambda r: mp.e**r, [-half, half], deg+1))]
-QLOG2E = qr(1 / mp.log(2))
-_ln2 = mp.log(2)
-QLN2HI = qr(mp.mpf(int(qr(_ln2) * mp.mpf(2)**80)) / mp.mpf(2)**80)
-QLN2LO = qr(_ln2 - QLN2HI)
-def exp_q(x, co):
+    def Pfun(t):
+        r = mp.sqrt(t)
+        if r == 0: return mp.mpf(1)/6
+        E = mp.e**r - 1 - r; c = 2*E/(r+E); return (r-c)/(r*r)
+    return [qr(c) for c in reversed(mp.chebyfit(Pfun, [mp.mpf('1e-40'), half*half], ncoef))]
+def exp_q(x, P):                               # fdlibm rational reconstruction
     x = qr(x); k = int(mp.floor(qr(x * QLOG2E) + 0.5))
-    r = qsub(qsub(x, qmul(mp.mpf(k), QLN2HI)), qmul(mp.mpf(k), QLN2LO))
-    return qr(qhorner(co, r) * mp.mpf(2)**k)
+    hi = qsub(x, qmul(mp.mpf(k), QLN2HI)); lo = qmul(mp.mpf(k), QLN2LO); r = qsub(hi, lo)
+    t = qmul(r, r); c = qsub(r, qmul(t, qhorner(P, t)))          # c = r - t*P(t)
+    y = qsub(mp.mpf(1), qsub(qsub(lo, qdiv(qmul(r, c), qsub(mp.mpf(2), c))), hi))
+    return qr(y * mp.mpf(2)**k)
 def check_exp_rq():
-    deg = 24; co = gen_exp_coeffs_q(deg)
-    print(f"# exp @rq: Cody-Waite + degree-{deg} minimax poly (f128, 113-bit)")
+    P = gen_P_q(11)                            # deg-10 even poly = 11 coeffs
+    print("# exp @rq: Cody-Waite + fdlibm rational reconstruction")
+    print("#   exp(r) = 1 - ((lo - r*c/(2-c)) - hi),  c = r - t*P(t),  t = r*r")
     print(f"LOG2E={qhex(QLOG2E)}\nLN2HI={qhex(QLN2HI)}\nLN2LO={qhex(QLN2LO)}")
-    print("coeffs c0..c%d:" % deg)
-    for i, c in enumerate(co): print(f"  c{i:<2}={qhex(c)}")
+    print("P coeffs (ascending in t=r^2):")
+    for i, c in enumerate(P): print(f"  p{i:<2}={qhex(c)}")
     worst = 0
-    for t in range(-2000, 2001):
-        x = qr(mp.mpf(t)/100); g = exp_q(x, co); tr = mp.e**x
+    for t in range(-20000, 20001, 3):          # [-20,20] step 0.003 (rq_check.c+MPFR is authoritative)
+        x = qr(mp.mpf(t)/1000); g = exp_q(x, P); tr = mp.e**x
         ulp = mp.mpf(2)**(mp.frexp(g)[1] - QP); worst = max(worst, abs((g-tr)/ulp))
-    print(f"max error over [-20,20]: {float(worst):.3f} ULP")
+    print(f"max error over [-20,20]: {float(worst):.3f} ULP  (faithful; fdlibm beats the old flat Horner)")
     for v in ['1','0.5','-2','10']:
-        print(f"  exp({v}) -> {qhex(exp_q(mp.mpf(v), co))}")
+        print(f"  exp({v}) -> {qhex(exp_q(mp.mpf(v), P))}")
 
 if __name__ == '__main__':
     fn = sys.argv[1] if len(sys.argv) > 1 else 'exp'

libmath/tools/rd_exp_check.c

@@ -0,0 +1,76 @@
+// rd_exp_check.c -- SoftFloat-f64 reference for @rd exp (the jet body).
+// Transcribed line-for-line from libmath/desk/lib/math.hoon ++rd ++exp
+// (Cody-Waite reduction + degree-11 minimax Horner + scale2 ldexp).
+// SoftFloat f64 round-near-even == the Hoon stdlib @rd ops, so this should be
+// BIT-EXACT to (exp:rd:math x).  Prints "<in_hex> <out_hex>" per input.
+//
+// Build: re-archive the zig softfloat .a (not 8-byte aligned for Apple ld):
+//   ar x <zig libsoftfloat.a> && libtool -static -o /tmp/libsoftfloat.a *.o
+//   cc -I<sf-include> rd_exp_check.c /tmp/libsoftfloat.a -o /tmp/rd_exp_check
+#include <stdio.h>
+#include <stdint.h>
+#include <stdbool.h>
+#include <stdlib.h>
+#include "softfloat.h"
+
+typedef float64_t d;
+static inline d   D(uint64_t b){ d r; r.v = b; return r; }
+static inline d   mul(d a, d b){ return f64_mul(a,b); }
+static inline d   sub(d a, d b){ return f64_sub(a,b); }
+static inline d   add(d a, d b){ return f64_add(a,b); }
+
+// pow2(j) = (j+1023)<<52 as f64 bits = 2^j  (normal range j in [-1022,1023])
+static inline d pow2(int64_t j){ return D(((uint64_t)(j + 1023)) << 52); }
+
+// scale2: correctly-rounded ldexp with overflow/subnormal tails (math.hoon:1519)
+static d scale2(d p, int64_t k){
+  if ( k - 1024 >= 0 )                       // k>=1024
+    return mul(mul(p, pow2(1023)), pow2(k - 1023));
+  if ( !(k + 1022 >= 0) )                     // k<-1022  (i.e. k<=-1023)
+    return mul(mul(p, pow2(k + 54)), pow2(-54));
+  return mul(p, pow2(k));
+}
+
+static d expd(d x){
+  const uint64_t QNAN = 0x7ff8000000000000ULL, PINF = 0x7ff0000000000000ULL,
+                 NINF = 0xfff0000000000000ULL;
+  if ( !f64_eq(x, x) )       return D(QNAN);            // NaN
+  if ( x.v == PINF )         return D(PINF);            // +inf
+  if ( x.v == NINF )         return D(0);               // -inf -> 0
+  d log2e = D(0x3ff71547652b82feULL);
+  d ln2hi = D(0x3fe62e42fee00000ULL);
+  d ln2lo = D(0x3dea39ef35793c76ULL);
+  int64_t k = f64_to_i64(mul(x, log2e), softfloat_round_near_even, false);
+  if ( k - 1025 >= 0 )       return D(PINF);            // overflow -> inf
+  if ( !(k + 1075 >= 0) )    return D(0);               // underflow -> 0
+  d ka = ui64_to_f64( (uint64_t)(k < 0 ? -k : k) );     // |k| as f64 (sun abs)
+  d kf = (k >= 0) ? ka : sub(D(0), ka);
+  d r  = sub( sub(x, mul(kf, ln2hi)), mul(kf, ln2lo) );
+  // degree-11 minimax coeffs c0..c11 (math.hoon:1542-1547)
+  static const uint64_t cs[12] = {
+    0x3ff0000000000000ULL, 0x3ff0000000000000ULL, 0x3fe0000000000011ULL,
+    0x3fc555555555555aULL, 0x3fa555555554f0cfULL, 0x3f8111111110f225ULL,
+    0x3f56c16c187fbe02ULL, 0x3f2a01a01b14378fULL, 0x3efa01991ac8730aULL,
+    0x3ec71ddf5749d126ULL, 0x3e928b4057f44145ULL, 0x3e5af631d0059becULL
+  };
+  // Horner: acc=0; for c in flop(cs) [c11..c0]: acc = acc*r + c   (math.hoon:1549)
+  d p = D(0);
+  for ( int i = 11; i >= 0; i-- ) p = add(mul(p, r), D(cs[i]));
+  return scale2(p, k);
+}
+
+int main(int argc, char** argv){
+  softfloat_roundingMode = softfloat_round_near_even;
+  // inputs as native doubles (native == IEEE f64 bit patterns here)
+  double xs[] = { 0.0, 0.5, 1.0, 2.0, 3.0, -1.0, -2.0, 5.0, 10.0, -10.0,
+                  0.1, -0.1, 3.141592653589793, 100.0, -100.0,
+                  700.0, -700.0, 709.0, -740.0 };
+  int n = (int)(sizeof xs / sizeof xs[0]);
+  for ( int i = 0; i < n; i++ ){
+    union { double f; uint64_t b; } u; u.f = xs[i];
+    d out = expd(D(u.b));
+    printf("0x%016llx 0x%016llx  (x=%g)\n",
+           (unsigned long long)u.b, (unsigned long long)out.v, xs[i]);
+  }
+  return 0;
+}

libmath/tools/rq_check.c

@@ -29,10 +29,21 @@ static double ulp_err(q g,mpfr_t tr){
   mpfr_sub(d,gv,tr,MPFR_RNDN); long e2; mpfr_get_d_2exp(&e2,gv,MPFR_RNDN);
   mpfr_set_ui(u,1,MPFR_RNDN); mpfr_mul_2si(u,u,(int)e2-1-112,MPFR_RNDN); mpfr_div(d,d,u,MPFR_RNDN);
   double r=mpfr_get_d(d,MPFR_RNDN); if(r<0)r=-r; mpfr_clears(gv,d,u,(mpfr_ptr)0); return r;}
-static const q LOG2E={{0xe1777d0ffda0d23aull,0x3fff71547652b82full}},LN2HI={{0xf35793c800000000ull,0x3ffe62e42fefa39eull}},LN2LO={{0xfc32f366359d274aull,0xbfad319ff0342542ull}},HALF={{0x0000000000000000ull,0x3ffe000000000000ull}},ONE={{0x0000000000000000ull,0x3fff000000000000ull}},SQRT2={{0xc908b2fb1366ea95ull,0x3fff6a09e667f3bcull}};
+static const q LOG2E={{0xe1777d0ffda0d23aull,0x3fff71547652b82full}},LN2HI={{0xf35793c800000000ull,0x3ffe62e42fefa39eull}},LN2LO={{0xfc32f366359d274aull,0xbfad319ff0342542ull}},HALF={{0x0000000000000000ull,0x3ffe000000000000ull}},ONE={{0x0000000000000000ull,0x3fff000000000000ull}},TWO={{0x0000000000000000ull,0x4000000000000000ull}},SQRT2={{0xc908b2fb1366ea95ull,0x3fff6a09e667f3bcull}};
 static const q INVPIO2={{0x2a53f84eafa3ea6aull,0x3ffe45f306dc9c88ull}},PIO2_1={{0x8460000000000000ull,0x3fff921fb54442d1ull}},PIO2_1T={{0x07344a409382229aull,0x3fc2313198a2e037ull}};
-static const q EXC[25]={
-  {{0x0000000000000000ull,0x3fff000000000000ull}}, {{0x0000000000000000ull,0x3fff000000000000ull}}, {{0x0000000000000000ull,0x3ffe000000000000ull}}, {{0x5555555555555555ull,0x3ffc555555555555ull}}, {{0x5555555555555555ull,0x3ffa555555555555ull}}, {{0x1111111111111111ull,0x3ff8111111111111ull}}, {{0x6c16c16c16c16c17ull,0x3ff56c16c16c16c1ull}}, {{0xa01a01a01a01a3e8ull,0x3ff2a01a01a01a01ull}}, {{0xa01a01a01a01a146ull,0x3fefa01a01a01a01ull}}, {{0x38faac1c88a5a526ull,0x3fec71de3a556c73ull}}, {{0xc72ef016d3d6e867ull,0x3fe927e4fb7789f5ull}}, {{0x38fe748363c46e8bull,0x3fe5ae64567f544eull}}, {{0x7b544dab18f475c5ull,0x3fe21eed8eff8d89ull}}, {{0x97c9f3aebabb2423ull,0x3fde6124613a86d0ull}}, {{0xd20b83c7f94d17d8ull,0x3fda93974a8c07c9ull}}, {{0xf5f4284f0d74f9e7ull,0x3fd6ae7f3e733b81ull}}, {{0xf417b4d27c5f92a9ull,0x3fd2ae7f3e733b81ull}}, {{0x6a419e674779c97cull,0x3fce952c77030a99ull}}, {{0x0466ff8c8b42b3dfull,0x3fca6827863b97b5ull}}, {{0x874b7a686d819241ull,0x3fc62f49b469f892ull}}, {{0xbb3b32a11bb5f139ull,0x3fc1e542ba427463ull}}, {{0xc93890ff9ab55cbbull,0x3fbd71b8db9f7f73ull}}, {{0x6efc0717eae785a1ull,0x3fb90ce38aab7bd7ull}}, {{0xcb3f4f7edfaa2666ull,0x3fb47693274bab2aull}}, {{0x61cb0e23655d47cbull,0x3faff3629154e0a7ull}}
+// exp: fdlibm rational reconstruction; EXC = even minimax P(t), t=r^2 (deg-10)
+static const q EXC[11]={
+  {{0x5555555555555555ull,0x3ffc555555555555ull}},
+  {{0x6c16c16c16c09e83ull,0xbff66c16c16c16c1ull}},
+  {{0x6abc0115453d96ddull,0x3ff11566abc01156ull}},
+  {{0xaac663e4a6d65ccaull,0xbfebbbd779334ef0ull}},
+  {{0xda06115986f507fbull,0x3fe666a8f2bf70ebull}},
+  {{0x43eb0e288c2e45a8ull,0xbfe122805d644267ull}},
+  {{0x12be0476b628552full,0x3fdbd6db2c4e0507ull}},
+  {{0xeb838f5da821635aull,0xbfd67da4e1efb419ull}},
+  {{0xdc61daecbfc0d781ull,0x3fd1355867f7df64ull}},
+  {{0x54bb7852bc52bd9aull,0xbfcbf56e4264f8adull}},
+  {{0x822162270789ca71ull,0x3fc68fc13579bfe0ull}}
 };
 static const q LOGC[23]={
   {{0x5555555555555555ull,0x3ffd555555555555ull}}, {{0x999999999999999aull,0x3ffc999999999999ull}}, {{0x2492492492492492ull,0x3ffc249249249249ull}}, {{0xc71c71c71c71c71cull,0x3ffbc71c71c71c71ull}}, {{0x5d1745d1745d1746ull,0x3ffb745d1745d174ull}}, {{0x3b13b13b13b13b14ull,0x3ffb3b13b13b13b1ull}}, {{0x1111111111111111ull,0x3ffb111111111111ull}}, {{0xe1e1e1e1e1e1e1e2ull,0x3ffae1e1e1e1e1e1ull}}, {{0x86bca1af286bca1bull,0x3ffaaf286bca1af2ull}}, {{0x8618618618618618ull,0x3ffa861861861861ull}}, {{0x42c8590b21642c86ull,0x3ffa642c8590b216ull}}, {{0xae147ae147ae147bull,0x3ffa47ae147ae147ull}}, {{0x84bda12f684bda13ull,0x3ffa2f684bda12f6ull}}, {{0x611a7b9611a7b961ull,0x3ffa1a7b9611a7b9ull}}, {{0x4210842108421084ull,0x3ffa084210842108ull}}, {{0x7c1f07c1f07c1f08ull,0x3ff9f07c1f07c1f0ull}}, {{0xd41d41d41d41d41dull,0x3ff9d41d41d41d41ull}}, {{0xf914c1bacf914c1cull,0x3ff9bacf914c1bacull}}, {{0xa41a41a41a41a41aull,0x3ff9a41a41a41a41ull}}, {{0x9c18f9c18f9c18faull,0x3ff98f9c18f9c18full}}, {{0x417d05f417d05f41ull,0x3ff97d05f417d05full}}, {{0x6c16c16c16c16c17ull,0x3ff96c16c16c16c1ull}}, {{0x72620ae4c415c988ull,0x3ff95c9882b93105ull}}
@@ -44,9 +55,13 @@ static const q COSC[16]={
   {{0x5555555555555555ull,0x3ffa555555555555ull}}, {{0x6c16c16c16c16c17ull,0xbff56c16c16c16c1ull}}, {{0xa01a01a01a01a01aull,0x3fefa01a01a01a01ull}}, {{0xc72ef016d3ea6679ull,0xbfe927e4fb7789f5ull}}, {{0x7b544da987acfe85ull,0x3fe21eed8eff8d89ull}}, {{0xd20badf145dfa3e5ull,0xbfda93974a8c07c9ull}}, {{0xf11d8656b0ee8cb0ull,0x3fd2ae7f3e733b81ull}}, {{0x77bb004886a2c2abull,0xbfca6827863b97d9ull}}, {{0x507a9cad2bf8f0bbull,0x3fc1e542ba402022ull}}, {{0x29450c90b7f338ecull,0xbfb90ce396db7f85ull}}, {{0x7cca4b4067ca9d8aull,0x3faff2cf01972f57ull}}, {{0x9a38f2050ba6b015ull,0xbfa688e85fc6a4e5ull}}, {{0xd373c5c51c354a8dull,0x3f9d0a18a2635085ull}}, {{0xe60caded4c2989c5ull,0xbf933932c5047d60ull}}, {{0xc42e1ee46fa6bfc4ull,0x3f89434d2e783f5bull}}, {{0xa13f8a2b4af9d6b7ull,0xbf7f2710231c0fd7ull}}
 };
 static q expq(q x){
-  q t=qadd(qmul(x,LOG2E),HALF); int64_t k=f128M_to_i64(&t,softfloat_round_min,false); q qk=qi(k);
-  q r=qsub(qsub(x,qmul(qk,LN2HI)),qmul(qk,LN2LO));
-  q p=EXC[24]; for(int i=23;i>=0;i--) p=qadd(qmul(p,r),EXC[i]); return qmul(p,q2k((int)k));}
+  q t0=qadd(qmul(x,LOG2E),HALF); int64_t k=f128M_to_i64(&t0,softfloat_round_min,false); q qk=qi(k);
+  q hi=qsub(x,qmul(qk,LN2HI)), lo=qmul(qk,LN2LO), r=qsub(hi,lo);
+  q t=qmul(r,r);
+  q c=EXC[10]; for(int i=9;i>=0;i--) c=qadd(qmul(c,t),EXC[i]);   // P(t)
+  c=qsub(r,qmul(t,c));                                           // c = r - t*P
+  q y=qsub(ONE, qsub(qsub(lo, qdiv(qmul(r,c),qsub(TWO,c))), hi));// 1-((lo-r*c/(2-c))-hi)
+  return qmul(y,q2k((int)k));}
 static q logq(q x){
   uint64_t hi=x.v[1]; int ef=((hi>>48)&0x7fff)-16383;
   q m; m.v[0]=x.v[0]; m.v[1]=(hi&0xffffffffffffULL)|(16383ULL<<48);