Merge pull request #138 from herbie-fp/direct-rounding-pow

Direct rounding for `ival-pow` and `ival-pow2`

Author: obround

mpfr.rkt

@@ -138,6 +138,8 @@
 (define mpfr-atan2!
   (get-mpfr-fun 'mpfr_atan2 (_fun _mpfr-pointer _mpfr-pointer _mpfr-pointer _rnd_t -> _int)))
 
+(define mpfr-pow2! (get-mpfr-fun 'mpfr_sqr (_fun _mpfr-pointer _mpfr-pointer _rnd_t -> _int)))
+
 (define mpfr-set! (get-mpfr-fun 'mpfr_set (_fun _mpfr-pointer _mpfr-pointer _rnd_t -> _void)))
 
 (define mpfr-set-prec! set-mpfr-prec!)
@@ -346,5 +348,6 @@
          mpfr-log!
          mpfr-const-pi!
          mpfr-atan2!
+         mpfr-pow2!
          mpfr-set-prec!
          mpfr-set!)

ops/core.rkt

@@ -32,6 +32,12 @@
          ival-union
          ival-split
          (rename-out [monotonic monotonic->ival] [comonotonic comonotonic->ival])
+         define*
+         immutable
+         monotonic-mpfr!
+         comonotonic-mpfr!
+         monotonic-mpfr
+         comonotonic-mpfr
          ival-illegal
          ival-pi
          ival-e

ops/pow.rkt

@@ -13,15 +13,15 @@
     [(and (< (mpfr-exp (ival-hi-val x)) 1) (not (bfinfinite? (ival-hi-val x)))) -1]
     [else 0]))
 
-(define (eppow a a! b b! a-class b-class)
-  (when (bfzero? a)
-    (set! a 0.bf)) ; Handle (-0)^(-1)
-  (define-values (val exact?) (bf-return-exact? bfexpt (list a b)))
-  (values val
+(define (eppow! out a a! b b! a-class b-class rnd)
+  (define base (if (bfzero? a) 0.bf a)) ; Handle (-0)^(-1)
+  (mpfr-set-prec! out (bf-precision))
+  (define exact? (= 0 (mpfr-pow! out base b rnd)))
+  (values out
           (or (and a! b! exact?)
-              (and a! (bf=? a 1.bf))
-              (and a! (bfzero? a) (not (= b-class 0)))
-              (and a! (bfinfinite? a) (not (= b-class 0)))
+              (and a! (bf=? base 1.bf))
+              (and a! (bfzero? base) (not (= b-class 0)))
+              (and a! (bfinfinite? base) (not (= b-class 0)))
               (and b! (bfzero? b))
               (and b! (bfinfinite? b) (not (= a-class 0))))))
 
@@ -39,24 +39,25 @@
   (define y-class (classify-ival y))
 
   (define (mk-pow a b c d)
+    (define out (new-ival (bf-precision)))
     (define-values (lo lo!)
-      (rnd 'down
-           eppow
-           (endpoint-val a)
-           (endpoint-immovable? a)
-           (endpoint-val b)
-           (endpoint-immovable? b)
-           x-class
-           y-class))
+      (eppow! (ival-lo-val out)
+              (endpoint-val a)
+              (endpoint-immovable? a)
+              (endpoint-val b)
+              (endpoint-immovable? b)
+              x-class
+              y-class
+              'down))
     (define-values (hi hi!)
-      (rnd 'up
-           eppow
-           (endpoint-val c)
-           (endpoint-immovable? c)
-           (endpoint-val d)
-           (endpoint-immovable? d)
-           x-class
-           y-class))
+      (eppow! (ival-hi-val out)
+              (endpoint-val c)
+              (endpoint-immovable? c)
+              (endpoint-val d)
+              (endpoint-immovable? d)
+              x-class
+              y-class
+              'up))
 
     (define-values (real-lo! real-hi!)
       (cond
@@ -69,17 +70,16 @@
          ;; Can't use >=, even though exp2-overflow-threshold is a
          ;; power of 2, because mpfr-exp is offset by 1 from the real
          ;; exponent, which matters when we add them.
+         (define log2-base (bf 0))
+         (define (log2-sum-exceeds-threshold? exponent-value base-value)
+           (mpfr-log2! log2-base base-value 'zero)
+           (> (+ (mpfr-exp exponent-value) (mpfr-exp log2-base)) (mpfr-exp exp2-overflow-threshold)))
+
          (define must-overflow
-           (and (bfinfinite? hi)
-                (= (* x-class y-class) 1)
-                (> (+ (mpfr-exp bval) (mpfr-exp (rnd 'zero bflog2 aval)))
-                   (mpfr-exp exp2-overflow-threshold))))
+           (and (bfinfinite? hi) (= (* x-class y-class) 1) (log2-sum-exceeds-threshold? bval aval)))
 
          (define must-underflow
-           (and (bfzero? lo)
-                (= (* x-class y-class) -1)
-                (> (+ (mpfr-exp dval) (mpfr-exp (rnd 'zero bflog2 cval)))
-                   (mpfr-exp exp2-overflow-threshold))))
+           (and (bfzero? lo) (= (* x-class y-class) -1) (log2-sum-exceeds-threshold? dval cval)))
 
          (define real-lo! (or lo! must-underflow (and (bfzero? lo) a! b!)))
          (define real-hi! (or hi! must-underflow must-overflow (and (bfinfinite? hi) c! d!)))
@@ -152,8 +152,8 @@
       (let ([pospow (ival-pow-pos (ival-exact-fabs x) y)])
         (ival-then (ival-assert ival-maybe) (ival-union (ival-neg pospow) pospow)))))
 
-(define (ival-pow2 x)
-  ((monotonic->ival (lambda (x) (bfmul x x))) (ival-pre-fabs x)))
+(define* ival-pow2! (lambda (out x) ((monotonic-mpfr! mpfr-pow2!) out (ival-pre-fabs x))))
+(define* ival-pow2 (immutable ival-pow2!))
 
 (define (ival-pow x y)
   (cond