Simple-minded criterion to compute x^-y by 1/x^y was too simple.

For any |x| > 1, we can't do this without comparing the size of
log(x) and y well enough to discern what might yield subnormal, so
might as well just compute the general case then.

diff --git a/src/runtime/arith.scm b/src/runtime/arith.scm
index 134e9990851be813d8e62bea1a929744ae0b34f2..71c66b39967a4a1cb01c6e6612955fdd331f31bb 100644 (file)
--- a/src/runtime/arith.scm
+++ b/src/runtime/arith.scm
@@ -1265,12 +1265,15 @@ USA.
                                      (loop x y answer)))))))))
                 (cond ((int:positive? y) (exact-method y))
                       ((int:negative? y)
-                       (if (int:< y microcode-id/floating-exponent-min)
-                           ;; The exact method cannot generate
-                           ;; subnormals because the negated exponent
-                           ;; overflows, so use the general case.
-                           (general-case x (int:->flonum y))
-                           (flo:/ flo:1 (exact-method (int:negate y)))))
+                       (if (flo:< (flo:abs x) flo:1)
+                           (flo:/ flo:1 (exact-method (int:negate y)))
+                           ;; For any base above 1 and sufficiently large
+                           ;; exponents, or for any negative exponent and
+                           ;; sufficiently large base, the division above
+                           ;; would overflow into infinite and then yield
+                           ;; 0 when it should yield a subnormal.  So use
+                           ;; the general case.
+                           (general-case x (int:->flonum y))))
                       (else flo:1))))
              (else
               (general-case x (rat:->inexact y))))

diff --git a/tests/runtime/test-arith.scm b/tests/runtime/test-arith.scm
index 7017c7e0e76e8ea302925666eb6a072fe721c2c4..9384402e94c4450c20815560b75eedccae712e6e 100644 (file)
--- a/tests/runtime/test-arith.scm
+++ b/tests/runtime/test-arith.scm
@@ -245,4 +245,8 @@ USA.
         (let ((x (vector-ref v 0))
               (y (vector-ref v 1))
               (x^y (string->number (vector-ref v 2))))
-          (assert-eqv (expt x y) x^y))))))
\ No newline at end of file
+          (assert-eqv (expt x y) x^y)
+          ;; For all the inputs, reciprocal is exact.
+          (assert-eqv (expt (/ 1 x) (- y)) x^y)
+          (assert-eqv (expt (* 2 x) (/ y 2)) x^y)
+          (assert-eqv (expt (/ 1 (* 2 x)) (- (/ y 2))) x^y))))))
\ No newline at end of file