Merge branch 'grow-libm' into growlibm-newbench

Author: JonasRegehr

.gitignore

@@ -26,3 +26,5 @@ herbie-compiled/
 
 # Python
 .env/
+
+.vscode/

bench/omerc.fpcore

@@ -1,110 +1,131 @@
-; Hard benchmarks extracted from growlibm/omerc.cpp Oblique Mercator setup, focusing on tight domain guards and coupled trig/hyperbolic transforms.
+(FPCore (W)
+ :name "omerc-forward-half-diff"
+ :pre (not (= W 0.0))
+ (* 0.5 (- W (/ 1.0 W))))
 
-(FPCore (S V singam cosgam T)
-  :name "omerc-hard-u-mix"
-  :pre (and (> (fabs T) 1e-12)
-            (< (fabs singam) 1)
-            (< (fabs cosgam) 1)
-            (< (fabs V) 1))
-  (/ (- (* S singam) (* V cosgam)) T))
+(FPCore (W)
+ :name "omerc-forward-half-sum"
+ :pre (not (= W 0.0))
+ (* 0.5 (+ W (/ 1.0 W))))
+
+(FPCore (S T singam cosgam B lam)
+ :name "omerc-forward-U"
+ :pre (not (= T 0.0))
+ (let* ((V (sin (* B lam))))
+   (/ (- (* S singam) (* V cosgam)) T)))
 
 (FPCore (U ArB)
-  :name "omerc-hard-log-ratio"
-  :pre (and (< (fabs U) 0.999999)
-            (> ArB 0))
-  (* 0.5
-     (* ArB
-        (log (/ (- 1 U)
-                 (+ 1 U))))))
-
-(FPCore (ArB S V singam cosgam temp)
-  :name "omerc-hard-atan2-rot"
-  :pre (and (>= (fabs temp) 1e-7)
-            (< (fabs singam) 1)
-            (< (fabs cosgam) 1)
-            (< (fabs V) 1))
-  (* ArB
-     (atan2 (+ (* S cosgam)
-                (* V singam))
-            temp)))
+ :name "omerc-forward-v"
+ :pre (< (fabs U) 1.0)
+ (* 0.5 (* ArB (log (/ (- 1.0 U) (+ 1.0 U))))))
+
+(FPCore (S B lam cosgam singam ArB)
+ :name "omerc-forward-u-atan"
+ :pre (>= (fabs (cos (* B lam))) 1e-7)
+ (let* ((V (sin (* B lam)))
+        (temp (cos (* B lam)))
+        (numer (+ (* S cosgam) (* V singam))))
+   (* ArB (atan2 numer temp))))
+
+(FPCore (A lam)
+ :name "omerc-forward-u-linear"
+ (* A lam))
+
+(FPCore (u v cosrot sinrot u0)
+ :name "omerc-forward-rot-x"
+ (let* ((udiff (- u u0)))
+   (+ (* v cosrot) (* udiff sinrot))))
+
+(FPCore (u v cosrot sinrot u0)
+ :name "omerc-forward-rot-y"
+ (let* ((udiff (- u u0)))
+   (- (* udiff cosrot) (* v sinrot))))
+
+(FPCore (BrA v)
+ :name "omerc-inverse-exp"
+ (exp (- (* BrA v))))
 
 (FPCore (Qp)
-  :name "omerc-hard-sinh-sp"
-  :pre (> Qp 0)
-  (* 0.5
-     (- Qp
-        (/ 1 Qp))))
+ :name "omerc-inverse-Sp"
+ :pre (not (= Qp 0.0))
+ (* 0.5 (- Qp (/ 1.0 Qp))))
 
 (FPCore (Qp)
-  :name "omerc-hard-cosh-tp"
-  :pre (> Qp 0)
-  (* 0.5
-     (+ Qp
-        (/ 1 Qp))))
-
-(FPCore (Vp Sp singam cosgam Tp)
-  :name "omerc-hard-up-frac"
-  :pre (and (> (fabs Tp) 1e-12)
-            (< (fabs singam) 1)
-            (< (fabs cosgam) 1)
-            (< (fabs Vp) 1))
-  (/ (+ (* Vp cosgam)
-        (* Sp singam))
-     Tp))
-
-(FPCore (Esc Up)
-  :name "omerc-hard-phi-scale"
-  :pre (and (> Esc 0)
-            (< (fabs Up) 0.999999))
-  (/ Esc
-     (sqrt (/ (+ 1 Up)
-              (- 1 Up)))))
-
-(FPCore (rB Sp Vp singam cosgam BrAu)
-  :name "omerc-hard-lambda-atan2"
-  :pre (and (< (fabs singam) 1)
-            (< (fabs cosgam) 1)
-            (< (fabs Vp) 1))
-  (* -1
-     (* rB
-        (atan2 (- (* Sp cosgam)
-                   (* Vp singam))
-                (cos BrAu)))))
-
-(FPCore (B lamDiff denom)
-  :name "omerc-hard-gamma0-atan"
-  :pre (and (> (fabs denom) 1e-12)
-            (< (fabs lamDiff) 10))
-  (atan (/ (* 2 (sin (* B lamDiff)))
-          denom)))
-
-(FPCore (lamc F gamma0 B)
-  :name "omerc-hard-lam0-shift"
-  :pre (and (> (fabs B) 1e-12)
-            (> (fabs F) 1e-12)
-            (< (fabs (* (* 0.5 (- F (/ 1 F)))
-                        (tan gamma0)))
-               0.999999))
-  (let ((asinArg (* (* 0.5 (- F (/ 1 F)))
-                    (tan gamma0))))
-    (- lamc
-       (/ (asin asinArg)
-          B))))
-
-(FPCore (ArB D alpha_c)
-  :name "omerc-hard-u0-atan"
-  :pre (and (> ArB 0)
-            (>= (- (* D D) 1) 0)
-            (> (fabs (cos alpha_c)) 1e-12))
-  (fabs (* ArB
-           (atan (/ (sqrt (- (* D D) 1))
-                    (cos alpha_c))))))
-
-(FPCore (ArB gamma0 fortpi)
-  :name "omerc-hard-v-pole"
-  :pre (and (> ArB 0)
-            (< (fabs gamma0) 3)
-            (> (tan (- fortpi (* 0.5 gamma0))) 0))
-  (* ArB
-     (log (tan (- fortpi
-                   (* 0.5 gamma0))))))
+ :name "omerc-inverse-Tp"
+ :pre (not (= Qp 0.0))
+ (* 0.5 (+ Qp (/ 1.0 Qp))))
+
+(FPCore (Sp Tp Vp singam cosgam)
+ :name "omerc-inverse-Up"
+ :pre (not (= Tp 0.0))
+ (/ (+ (* Vp cosgam) (* Sp singam)) Tp))
+
+(FPCore (Up Esc)
+ :name "omerc-inverse-prephi"
+ :pre (< (fabs Up) 1.0)
+ (/ Esc (sqrt (/ (+ 1.0 Up) (- 1.0 Up)))))
+
+(FPCore (Sp Vp cosgam singam BrA u rB)
+ :name "omerc-inverse-lambda"
+ (let* ((temp (cos (* BrA u)))
+        (numer (- (* Sp cosgam) (* Vp singam))))
+   (- (* rB (atan2 numer temp)))))
+
+(FPCore (phi0 es)
+ :name "omerc-init-con"
+ :pre (and (>= es 0.0)
+           (< es 1.0)
+           (< (fabs phi0) (/ PI 2.0)))
+ (let* ((sinphi0 (sin phi0)))
+   (- 1.0 (* es (* sinphi0 sinphi0)))))
+
+(FPCore (phi0 es one_es)
+ :name "omerc-init-B"
+ :pre (and (>= es 0.0)
+           (< es 1.0)
+           (> one_es 0.0)
+           (< (fabs phi0) (/ PI 2.0)))
+ (let* ((cosphi0 (cos phi0))
+        (cos2 (* cosphi0 cosphi0))
+        (cos4 (* cos2 cos2)))
+   (sqrt (+ 1.0 (/ (* es cos4) one_es)))))
+
+(FPCore (phi0 es one_es com)
+ :name "omerc-init-D"
+ :pre (and (>= es 0.0)
+           (< es 1.0)
+           (> one_es 0.0)
+           (> com 0.0)
+           (< (fabs phi0) (/ PI 2.0))
+           (> (fabs (cos phi0)) 0.0)
+           (> (- 1.0 (* es (* (sin phi0) (sin phi0)))) 0.0))
+ (let* ((sinphi0 (sin phi0))
+        (cosphi0 (cos phi0))
+        (cos2 (* cosphi0 cosphi0))
+        (cos4 (* cos2 cos2))
+        (con (- 1.0 (* es (* sinphi0 sinphi0))))
+        (B (sqrt (+ 1.0 (/ (* es cos4) one_es))))
+        (root (sqrt con)))
+   (/ (* B com) (* cosphi0 root))))
+
+(FPCore (ArB D alpha_c phi0)
+ :name "omerc-init-u0"
+ :pre (and (>= (- (* D D) 1.0) 0.0)
+           (not (= (cos alpha_c) 0.0)))
+ (let* ((ratio (/ (sqrt (- (* D D) 1.0)) (cos alpha_c)))
+        (angle (atan ratio))
+        (scaled (* ArB angle))
+        (mag (fabs scaled)))
+   (if (< phi0 0.0)
+       (- mag)
+       mag)))
+
+(FPCore (ArB gamma0)
+ :name "omerc-init-v-pole-n"
+ :pre (> (tan (- (/ PI 4.0) (* 0.5 gamma0))) 0.0)
+ (* ArB (log (tan (- (/ PI 4.0) (* 0.5 gamma0))))))
+
+(FPCore (ArB gamma0)
+ :name "omerc-init-v-pole-s"
+ :pre (> (tan (+ (/ PI 4.0) (* 0.5 gamma0))) 0.0)
+ (* ArB (log (tan (+ (/ PI 4.0) (* 0.5 gamma0))))))

bench/tmerc.fpcore

@@ -0,0 +1,140 @@
+(FPCore (phi)
+ :name "approx-t-guarded"
+ (let* ((sinphi (sin phi))
+        (cosphi (cos phi)))
+   (if (> (fabs cosphi) 1e-10)
+       (/ sinphi cosphi)
+       0.0)))
+
+(FPCore (phi lam es)
+ :name "approx-scaled-lambda"
+ :pre (and (>= es 0.0)
+           (< es 1.0)
+           (> (- 1.0 (* es (* (sin phi) (sin phi)))) 0.0))
+ (let* ((sinphi (sin phi))
+        (cosphi (cos phi))
+        (den (- 1.0 (* es (* sinphi sinphi)))))
+   (/ (* cosphi lam)
+      (sqrt den))))
+
+(FPCore (phi es)
+ :name "approx-inv-meridional-con"
+ :pre (and (>= es 0.0)
+           (< es 1.0))
+ (let* ((sinphi (sin phi)))
+   (- 1.0 (* es (* sinphi sinphi)))))
+
+(FPCore (phi x es k0)
+ :name "approx-inv-d-term"
+ :pre (and (>= es 0.0)
+           (< es 1.0)
+           (> k0 0.0)
+           (> (- 1.0 (* es (* (sin phi) (sin phi)))) 0.0))
+ (let* ((sinphi (sin phi))
+        (con (- 1.0 (* es (* sinphi sinphi)))))
+   (/ (* x (sqrt con)) k0)))
+
+(FPCore (x esp)
+ :name "spherical-inv-g-from-exp"
+ :pre (> esp 0.0)
+ (let* ((h (exp (/ x esp))))
+   (* 0.5 (- h (/ 1.0 h)))))
+
+(FPCore (y esp phi0 g)
+ :name "spherical-inv-latitude-angle"
+ :pre (and (> esp 0.0)
+           (> (+ 1.0 (* g g)) 0.0)
+           (>= (- 1.0 (* (cos (+ phi0 (/ y esp)))
+                         (cos (+ phi0 (/ y esp))))) 0.0))
+ (let* ((D (+ phi0 (/ y esp)))
+        (cosD (cos D))
+        (numer (- 1.0 (* cosD cosD)))
+        (denom (+ 1.0 (* g g))))
+   (asin (sqrt (/ numer denom)))))
+
+(FPCore (y esp phi0 g)
+ :name "spherical-inv-longitude"
+ :pre (and (> esp 0.0)
+           (or (> (fabs g) 0.0)
+               (> (fabs (cos (+ phi0 (/ y esp)))) 0.0)))
+ (atan2 g (cos (+ phi0 (/ y esp)))))
+
+(FPCore (angle esp phi0)
+ :name "spherical-y-scale-shift"
+ :pre (> esp 0.0)
+ (* esp (- angle phi0)))
+
+(FPCore (Cn lam)
+ :name "exact-forward-conformal-lat"
+ (let* ((sC (sin Cn))
+        (cC (cos Cn))
+        (cE (cos lam)))
+   (atan2 sC (* cC cE))))
+
+(FPCore (Cn lam)
+ :name "exact-forward-tanCe"
+ :pre (> (+ (* (sin Cn) (sin Cn))
+            (* (cos Cn) (cos Cn) (* (cos lam) (cos lam)))) 0.0)
+ (let* ((sC (sin Cn))
+        (cC (cos Cn))
+        (sE (sin lam))
+        (cE (cos lam))
+        (den (sqrt (+ (* sC sC)
+                      (* cC cC (* cE cE))))))
+   (/ (* sE cC) den)))
+
+(FPCore (Cn lam)
+ :name "exact-forward-two-inv-denom"
+ :pre (> (+ (* (sin Cn) (sin Cn))
+            (* (cos Cn) (cos Cn) (* (cos lam) (cos lam)))) 0.0)
+ (let* ((sC (sin Cn))
+        (cC (cos Cn))
+        (cE (cos lam))
+        (den (sqrt (+ (* sC sC)
+                      (* cC cC (* cE cE))))))
+   (/ 2.0 den)))
+
+(FPCore (Ce)
+ :name "exact-inv-exp2Ce"
+ (exp (* 2.0 Ce)))
+
+(FPCore (Ce)
+ :name "exact-inv-half-inv-exp"
+ (let* ((exp2 (exp (* 2.0 Ce))))
+   (/ 0.5 exp2)))
+
+(FPCore (Ce)
+ :name "exact-inv-cosh-arg"
+ (let* ((exp2 (exp (* 2.0 Ce)))
+        (half (/ 0.5 exp2)))
+   (+ (* 0.5 exp2) half)))
+
+(FPCore (Cn Ce)
+ :name "exact-inv-rr"
+ (let* ((sinhCe (sinh Ce))
+        (cosC (cos Cn))
+        (sinC (sin Cn))
+        (mod (sqrt (+ (* sinhCe sinhCe)
+                      (* cosC cosC)))))
+   (sqrt (+ (* sinC sinC)
+            (* mod mod)))))
+
+(FPCore (cos_arg_r cosh_arg_i)
+ :name "clenshaw-initial-r"
+ (* 2.0 cos_arg_r cosh_arg_i))
+
+(FPCore (sin_arg_r sinh_arg_i)
+ :name "clenshaw-initial-i"
+ (* -2.0 sin_arg_r sinh_arg_i))
+
+(FPCore (sin_arg_r cosh_arg_i cos_arg_r sinh_arg_i hr hi)
+ :name "clenshaw-final-real"
+ (let* ((r (* sin_arg_r cosh_arg_i))
+        (i (* cos_arg_r sinh_arg_i)))
+   (- (* r hr) (* i hi))))
+
+(FPCore (sin_arg_r cosh_arg_i cos_arg_r sinh_arg_i hr hi)
+ :name "clenshaw-final-imag"
+ (let* ((r (* sin_arg_r cosh_arg_i))
+        (i (* cos_arg_r sinh_arg_i)))
+   (+ (* r hi) (* i hr))))