Forked from
cmucl / cmucl
1663 commits behind the upstream repository.
-
Raymond Toy authored
Also replae WITH-INXACT-EXCEPTION-ENABLED with WITH-FLOAT-TRAPS-ENABLED. All tests still pass, as expected.
Raymond Toy authoredAlso replae WITH-INXACT-EXCEPTION-ENABLED with WITH-FLOAT-TRAPS-ENABLED. All tests still pass, as expected.
fdlibm.lisp 24.89 KiB
;; Test that fdlibm routines signals exceptions as expected.
(defpackage :fdlibm-tests
(:use :cl :lisp-unit))
(in-package "FDLIBM-TESTS")
(defparameter *qnan*
(ext:with-float-traps-masked (:invalid)
(* 0 ext:double-float-positive-infinity))
"Some randon quiet MaN value")
(defparameter *snan*
(kernel:make-double-float #x7ff00000 1)
"A randon signaling MaN value")
(defmacro with-inexact-exception-enabled (&body body)
`(ext:with-float-traps-enabled (:inexact)
,@body))
(define-test %cosh.exceptions
(:tag :fdlibm)
(assert-error 'floating-point-overflow
(kernel:%cosh 1000d0))
(assert-error 'floating-point-overflow
(kernel:%cosh -1000d0))
(assert-error 'floating-point-invalid-operation
(kernel:%cosh *snan*))
(assert-true (ext:float-nan-p (kernel:%cosh *qnan*)))
;; Same, but with overflow's masked
(ext:with-float-traps-masked (:overflow)
(assert-equal ext:double-float-positive-infinity
(kernel:%cosh 1000d0))
(assert-equal ext:double-float-positive-infinity
(kernel:%cosh -1000d0))
(assert-equal ext:double-float-positive-infinity
(kernel:%cosh ext:double-float-positive-infinity))
(assert-equal ext:double-float-positive-infinity
(kernel:%cosh ext:double-float-negative-infinity)))
;; Test NaN
(ext:with-float-traps-masked (:invalid)
(assert-true (ext:float-nan-p (kernel:%cosh *snan*)))))
(define-test %sinh.exceptions
(:tag :fdlibm)
(assert-error 'floating-point-overflow
(kernel:%sinh 1000d0))
(assert-error 'floating-point-overflow
(kernel:%sinh -1000d0))
(assert-error 'floating-point-invalid-operation
(kernel:%sinh *snan*))
(assert-true (ext:float-nan-p (kernel:%sinh *qnan*)))
;; Same, but with overflow's masked
(ext:with-float-traps-masked (:overflow)
(assert-equal ext:double-float-positive-infinity
(kernel:%sinh 1000d0))
(assert-equal ext:double-float-negative-infinity
(kernel:%sinh -1000d0))
(assert-equal ext:double-float-positive-infinity
(kernel:%sinh ext:double-float-positive-infinity))
(assert-equal ext:double-float-negative-infinity
(kernel:%sinh ext:double-float-negative-infinity)))
;; Test NaN
(ext:with-float-traps-masked (:invalid)
(assert-true (ext:float-nan-p (kernel:%sinh *qnan*))))
;; sinh(x) = x for |x| < 2^-28. Should signal inexact unless x = 0.
(let ((x (scale-float 1d0 -29))
(x0 0d0))
(with-inexact-exception-enabled
;; This must not throw an inexact exception because the result
;; is exact when the arg is 0.
(assert-eql 0d0 (kernel:%sinh x0)))
(with-inexact-exception-enabled
;; This must throw an inexact exception for non-zero x even
;; though the result is exactly x.
(assert-error 'floating-point-inexact
(kernel:%sinh x)))))
(define-test %tanh.exceptions
(:tag :fdlibm)
(assert-true (ext:float-nan-p (kernel:%tanh *qnan*)))
(assert-error 'floating-point-invalid-operation
(kernel:%tanh *snan*))
(ext:with-float-traps-masked (:invalid)
(assert-true (ext:float-nan-p (kernel:%tanh *snan*))))
;; tanh(x) = +/- 1 for |x| > 22, raising inexact, always.
(let ((x 22.1d0))
(with-inexact-exception-enabled
;; This must throw an inexact exception for non-zero x even
;; though the result is exactly x.
(assert-error 'floating-point-inexact
(kernel:%tanh x)))))
(define-test %acosh.exceptions
(:tag :fdlibm)
(assert-error 'floating-point-overflow
(kernel:%acosh ext:double-float-positive-infinity))
(assert-error 'floating-point-invalid-operation
(kernel:%acosh 0d0))
(ext:with-float-traps-masked (:overflow)
(assert-equal ext:double-float-positive-infinity
(kernel:%acosh ext:double-float-positive-infinity)))
(ext:with-float-traps-masked (:invalid)
(assert-true (ext:float-nan-p (kernel:%acosh 0d0)))))
(define-test %asinh.exceptions
(:tag :fdlibm)
(assert-error 'floating-point-invalid-operation
(kernel:%asinh *snan*))
(assert-error 'floating-point-overflow
(kernel:%asinh ext:double-float-positive-infinity))
(assert-error 'floating-point-overflow
(kernel:%asinh ext:double-float-negative-infinity))
(assert-true (ext:float-nan-p (kernel:%asinh *qnan*)))
(ext:with-float-traps-masked (:overflow)
(assert-equal ext:double-float-positive-infinity
(kernel:%asinh ext:double-float-positive-infinity))
(assert-error ext:double-float-negative-infinity
(kernel:%asinh ext:double-float-negative-infinity)))
(ext:with-float-traps-masked (:invalid)
(assert-true (ext:float-nan-p (kernel:%asinh *snan*))))
(let ((x (scale-float 1d0 -29))
(x0 0d0))
(with-inexact-exception-enabled
;; This must not throw an inexact exception because the result
;; is exact when the arg is 0.
(assert-eql 0d0 (asinh x0)))
(with-inexact-exception-enabled
;; This must throw an inexact exception for non-zero x even
;; though the result is exactly x.
(assert-error 'floating-point-inexact
(asinh x)))))
(define-test %atanh.exceptions
(:tag :fdlibm)
(assert-error 'floating-point-invalid-operation
(kernel:%atanh 2d0))
(assert-error 'floating-point-invalid-operation
(kernel:%atanh -2d0))
(assert-error 'division-by-zero
(kernel:%atanh 1d0))
(assert-error 'division-by-zero
(kernel:%atanh -1d0))
(ext:with-float-traps-masked (:invalid)
(assert-true (ext:float-nan-p (kernel:%atanh 2d0)))
(assert-true (ext:float-nan-p (kernel:%atanh -2d0))))
(ext:with-float-traps-masked (:divide-by-zero)
(assert-equal ext:double-float-positive-infinity
(kernel:%atanh 1d0))
(assert-equal ext:double-float-negative-infinity
(kernel:%atanh -1d0))))
(define-test %expm1.exceptions
(:tag :fdlibm)
(assert-error 'floating-point-overflow
(kernel:%expm1 709.8d0))
(assert-equal ext:double-float-positive-infinity
(kernel:%expm1 ext:double-float-positive-infinity))
(assert-error 'floating-point-invalid-operation
(kernel:%expm1 *snan*))
(assert-true (ext:float-nan-p (kernel:%expm1 *qnan*)))
(ext:with-float-traps-masked (:overflow)
(assert-equal ext:double-float-positive-infinity
(kernel:%expm1 709.8d0))
)
(ext:with-float-traps-masked (:invalid)
(assert-true (ext::float-nan-p (kernel:%expm1 *snan*))))
;; expm1(x) = -1 for x < -56*log(2), signaling inexact
(let ((x (* -57 (log 2d0))))
(with-inexact-exception-enabled
(assert-error 'floating-point-inexact
(kernel:%expm1 x)))))
(define-test %log1p.exceptions
(:tag :fdlibm)
(assert-error 'floating-point-invalid-operation
(kernel:%log1p -2d0))
(assert-error 'floating-point-overflow
(kernel:%log1p -1d0))
(assert-true (ext:float-nan-p (kernel:%log1p *qnan*)))
(ext:with-float-traps-masked (:overflow)
(assert-equal ext:double-float-negative-infinity
(kernel:%log1p -1d0)))
(ext:with-float-traps-masked (:invalid)
(assert-true (ext:float-nan-p (kernel:%log1p *snan*))))
;; log1p(x) = x for |x| < 2^-54, signaling inexact except for x = 0.
(let ((x (scale-float 1d0 -55))
(x0 0d0))
(with-inexact-exception-enabled
;; This must not throw an inexact exception because the result
;; is exact when the arg is 0.
(assert-eql 0d0 (kernel:%log1p x0)))
(with-inexact-exception-enabled
;; This must throw an inexact exception for non-zero x even
;; though the result is exactly x.
(assert-error 'floating-point-inexact
(kernel:%log1p x)))))
(define-test %exp.exceptions
(:tag :fdlibm)
(assert-error 'floating-point-overflow
(kernel:%exp 710d0))
(assert-true (ext:float-nan-p (kernel:%exp *qnan*)))
(assert-error 'floating-point-invalid-operation
(kernel:%exp *snan*))
(assert-equal ext:double-float-positive-infinity
(kernel:%exp ext:double-float-positive-infinity))
(assert-equal 0d0
(kernel:%exp -1000d0))
(ext:with-float-traps-masked (:overflow)
(assert-equal ext:double-float-positive-infinity
(kernel:%exp 710d0)))
(let ((modes (ext:get-floating-point-modes)))
(unwind-protect
(progn
(ext:set-floating-point-modes :traps '(:underflow))
(assert-error 'floating-point-underflow
(kernel:%exp -1000d0)))
(apply #'ext:set-floating-point-modes modes)))
(let ((x (scale-float 1d0 -29))
(x0 0d0))
;; exp(x) = x, |x| < 2^-28, with inexact exception unlees x = 0
(with-inexact-exception-enabled
;; This must not throw an inexact exception because the result
;; is exact when the arg is 0.
(assert-eql 1d0 (kernel:%exp x0)))
(with-inexact-exception-enabled
;; This must throw an inexact exception for non-zero x even
;; though the result is exactly x.
(assert-error 'floating-point-inexact
(kernel:%exp x)))))
(define-test %log.exception
(:tag :fdlibm)
(assert-error 'division-by-zero
(kernel:%log 0d0))
(assert-error 'division-by-zero
(kernel:%log -0d0))
(assert-error 'floating-point-invalid-operation
(kernel:%log -1d0))
(assert-error 'floating-point-invalid-operation
(kernel:%log *snan*))
(assert-true (ext:float-nan-p (kernel:%log *qnan*)))
(ext:with-float-traps-masked (:divide-by-zero)
(assert-equal ext:double-float-negative-infinity
(kernel:%log 0d0))
(assert-equal ext:double-float-negative-infinity
(kernel:%log -0d0)))
(ext:with-float-traps-masked (:invalid)
(assert-true (ext:float-nan-p (kernel:%log -1d0)))
(assert-true (ext:float-nan-p (kernel:%log *snan*)))))
(define-test %acos.exceptions
(:tag :fdlibm)
(assert-error 'floating-point-invalid-operation
(kernel:%acos 2d0))
(assert-error 'floating-point-invalid-operation
(kernel:%acos -2d0))
(ext:with-float-traps-masked (:invalid)
(assert-true (ext:float-nan-p (kernel:%acos 2d0)))
(assert-true (ext:float-nan-p (kernel:%acos -2d0)))))
(define-test %asin.exceptions
(:tag :fdlibm)
(assert-error 'floating-point-invalid-operation
(kernel:%asin 2d0))
(assert-error 'floating-point-invalid-operation
(kernel:%asin -2d0))
(ext:with-float-traps-masked (:invalid)
(assert-true (ext:float-nan-p (kernel:%asin 2d0)))
(assert-true (ext:float-nan-p (kernel:%asin -2d0)))))
(define-test %atan.exceptions
(:tag :fdlibm)
(assert-error 'floating-point-invalid-operation
(kernel:%atan *snan*))
(assert-true (ext:float-nan-p (kernel:%atan *qnan*)))
(ext:with-float-traps-masked (:invalid)
(assert-true (ext:float-nan-p (kernel:%atan *snan*))))
;; atan(x) = x for |x| < 2^-29, signaling inexact.
(let ((x (scale-float 1d0 -30))
(x0 0d0))
(with-inexact-exception-enabled
;; This must not throw an inexact exception because the result
;; is exact when the arg is 0.
(assert-eql 0d0 (kernel:%atan x0)))
(with-inexact-exception-enabled
;; This must throw an inexact exception for non-zero x even
;; though the result is exactly x.
(assert-error 'floating-point-inexact
(kernel:%atan x)))))
(define-test %log10.exceptions
(:tag :fdlibm)
;; %log10(2^k) = k
(dotimes (k 23)
(assert-equalp k
(kernel:%log10 (float (expt 10 k) 1d0))))
(assert-error 'division-by-zero
(kernel:%log10 0d0))
(assert-error 'floating-point-invalid-operation
(kernel:%log10 -1d0))
(assert-true (ext:float-nan-p (kernel:%log10 *qnan*)))
(assert-equal ext:double-float-positive-infinity
(kernel:%log10 ext:double-float-positive-infinity))
(ext:with-float-traps-masked (:divide-by-zero)
(assert-equal ext:double-float-negative-infinity
(kernel:%log10 0d0))
(assert-equal ext:double-float-negative-infinity
(kernel:%log10 -0d0)))
(ext:with-float-traps-masked (:invalid)
(assert-true (ext:float-nan-p (kernel:%log10 -1d0)))))
(define-test %scalbn.exceptions
(:tag :fdlibm)
(let ((modes (ext:get-floating-point-modes)))
(unwind-protect
(progn
(ext:set-floating-point-modes :traps '(:underflow))
(assert-error 'floating-point-underflow
(kernel:%scalbn 1d0 -51000)))
(apply #'ext:set-floating-point-modes modes)))
(assert-true 0d0
(kernel:%scalbn 1d0 -51000))
(assert-true -0d0
(kernel:%scalbn -1d0 -51000))
(assert-error 'floating-point-overflow
(kernel:%scalbn ext:double-float-positive-infinity 1))
(assert-error 'floating-point-invalid-operation
(kernel:%scalbn *snan* 1))
(assert-error 'floating-point-overflow
(kernel:%scalbn most-positive-double-float 2))
(assert-error 'floating-point-overflow
(kernel:%scalbn most-negative-double-float 2))
(ext:with-float-traps-masked (:overflow)
(assert-equal ext:double-float-positive-infinity
(kernel:%scalbn ext:double-float-positive-infinity 1))
(assert-equal ext:double-float-positive-infinity
(kernel:%scalbn most-positive-double-float 2))
(assert-equal ext:double-float-negative-infinity
(kernel:%scalbn most-negative-double-float 2))))
;;; These tests taken from github.com/rtoy/fdlibm-js
(define-test acosh-basic-tests
(:tag :fdlibm)
;; acosh(1) = 0
(assert-eql 0d0 (acosh 1d0))
;; acosh(1.5) = log((sqrt(5)+3)/2, case 1 < x < 2
(assert-eql 0.9624236501192069d0 (acosh 1.5d0))
;; acosh(4) = log(sqrt(15)+4), case 2 < x < 2^28
(assert-eql 2.0634370688955608d0 (acosh 4d0))
;; acosh(2^50), case 2^28 < x
(assert-eql 35.35050620855721d0 (acosh (scale-float 1d0 50)))
;; No overflow for most positive
(assert-eql 710.4758600739439d0 (acosh most-positive-double-float)))
(define-test asinh-basic-tests
(:tag :fdlibm)
(assert-eql -0d0 (asinh -0d0))
(assert-eql 0d0 (asinh 0d0))
(let ((x (scale-float 1d0 -29))
(x0 0d0))
;; asinh(x) = x for x < 2^-28
(assert-eql x (asinh x))
(assert-eql (- x) (asinh (- x))))
(let ((x (scale-float 1d0 -28)))
;; Case 2 > |x| >= 2^-28
(assert-eql 3.725290298461914d-9 (asinh x))
(assert-eql -3.725290298461914d-9 (asinh (- x))))
(let ((x 1d0))
;; Case 2 > |x| >= 2^-28
(assert-eql 0.881373587019543d0 (asinh x))
(assert-eql -0.881373587019543d0 (asinh (- x))))
(let ((x 5d0))
;; Case 2^28 > |x| > 2
(assert-eql 2.3124383412727525d0 (asinh x))
(assert-eql -2.3124383412727525d0 (asinh (- x))))
(let ((x (scale-float 1d0 28)))
;; Case 2^28 > |x|
(assert-eql 20.101268236238415d0 (asinh x))
(assert-eql -20.101268236238415d0 (asinh (- x))))
(let ((x most-positive-double-float))
;; No overflow for most-positive-double-float
(assert-eql 710.4758600739439d0 (asinh x))
(assert-eql -710.4758600739439d0 (asinh (- x)))))
(define-test atanh-basic-tests
(:tag :fdlibm)
(assert-eql +0d0 (atanh +0d0))
(assert-eql -0d0 (atanh -0d0))
;; atanh(x) = x, |x| < 2^-28
(let ((x (scale-float 1d0 -29)))
(assert-eql x (atanh x))
(assert-eql (- x) (atanh (- x))))
;; atanh(0.25) = log(5/3)/2, |x| < 0.5
(let ((x 0.25d0))
(assert-eql 0.25541281188299536d0 (atanh x))
(assert-eql -0.25541281188299536d0 (atanh (- x)))
;; There's no guarantee that atanh(1/4) = log(5/3)2 in floating
;; point, but it's true in this case with fdlibm
(assert-eql (/ (log (float 5/3 1d0)) 2) (atanh x)))
;; atanh(0.75) = log(7)/2, 0.5 < |x| < 1
(let ((x 0.75d0))
(assert-eql 0.9729550745276566d0 (atanh x))
(assert-eql -0.9729550745276566d0 (atanh (- x)))
;; There's no guarantee that atanh(3/4) = log(7)2 in floating
;; point, but it's true in this case with fdlibm
(assert-eql (/ (log 7d0) 2) (atanh x))))
(define-test cosh-basic-tests
(:tag :fdlibm)
;; cosh(2^-55) = 1, tiny x case
(let ((x (scale-float 1d0 -55)))
(assert-eql 1d0 (cosh x))
(assert-eql 1d0 (cosh (- x))))
;; cosh(2^-55) = 1, tiny x case
(let ((x (scale-float 1d0 -56)))
(assert-eql 1d0 (cosh x))
(assert-eql 1d0 (cosh (- x))))
;; cosh(log(2)/4) = (sqrt(2) + 1)/2^(5/4), case |x| < log(2)/2
(let ((x (/ (log 2d0) 4)))
;; This depends on (/ (log 2d0) 4) producing the value we really
;; want as the arg.
(assert-eql 1.0150517651282178d0 (cosh x))
(assert-eql 1.0150517651282178d0 (cosh (- x))))
;; cosh(10*log(2)) = 1048577/2048, case log(2)/2 < |x| < 22
(let ((x (* 10 (log 2d0)))
(y (float 1048577/2048 1d0)))
(assert-eql y (cosh x))
(assert-eql y (cosh (- x))))
;; cosh(32*log(2)), case 22 <= |x| < log(maxdouble)
(let ((x (* 32 (log 2d0))))
(assert-eql 2.1474836479999983d9 (cosh x))
(assert-eql 2.1474836479999983d9 (cosh (- x))))
;; cosh(710.4758600739439), case log(maxdouble) <= |x| <= overflowthreshold
(let ((x 710.4758600739439d0))
(assert-eql 1.7976931348621744d308 (cosh x))
(assert-eql 1.7976931348621744d308 (cosh (- x)))))
(define-test exp-basic-tests
(:tag :fdlibm)
;; No overflow and no underflow
(let ((x 709.7822265625d0))
(assert-eql 1.7968190737295725d308 (exp x))
(assert-eql 5.565390609552841d-309 (exp (- x))))
;; exp(7.09782712893383973096e+02), no overflow
(assert-eql 1.7976931348622732d308 (exp 7.09782712893383973096d+02))
;; exp(-7.45133219101941108420e+02), no underflow
(assert-eql 4.9406564584124654d-324 (exp -7.45133219101941108420d+02))
;; Overflow
(assert-error 'floating-point-overflow (exp 709.7827128933841d0))
;; Case |x| < 2^-28
(let ((x (scale-float 1d0 -29)))
(assert-eql (+ 1 x) (exp x))
(assert-eql (- 1 x) (exp (- x))))
;; exp(0.5), case log(2)/2 < |x| < 3/2*log(2)
(let ((x 0.5d0))
(assert-eql 1.6487212707001282d0 (exp x))
(assert-eql 0.6065306597126334d0 (exp (- x))))
;; exp(2), case |x| > 3/2*log(2)
(let ((x 2d0))
(assert-eql 7.38905609893065d0 (exp x))
(assert-eql 0.1353352832366127d0 (exp (- x))))
;; exp(2^-1022), case k < -1021
(assert-eql 1d0 (exp (scale-float 1d0 -1022)))
;; exp(2^-1021), case k >= -1021
(assert-eql 1d0 (exp (scale-float 1d0 -1021)))
;; exp(7.09782712893383973096e+02), no overflow
(assert-eql 1.7976931348622732d308 (exp 7.09782712893383973096d+02))
;; overflow
(assert-error 'floating-point-overflow (exp 709.7827128933841d0))
;; exp(-7.45133219101941108420e+02), no underflow
(assert-eql 4.9406564584124654d-324 (exp -745.1332191019411d0))
;; exp(-745.1332191019412), underflows
(assert-eql 0d0 (exp -745.1332191019412d0))
;; exp(1000) overflow
(assert-error 'floating-point-overflow (exp 1000d0))
;; exp(-1000) underflow
(assert-eql 0d0 (exp -1000d0)))
(define-test log-basic-tests
(:tag :fdlibm)
(assert-eql 0d0 (log 1d0))
(assert-eql 1d0 (log (exp 1d0)))
(assert-eql -1d0 (log (exp -1d0)))
(assert-eql 0.5d0 (log (sqrt (exp 1d0))))
(assert-eql -0.5d0 (log (sqrt (exp -1d0))))
;; Test a denormal arg
(assert-eql -709.08956571282410d0 (log (scale-float 1d0 -1023)))
;; Largest double value
(assert-eql 709.7827128933840d0 (log most-positive-double-float))
;; Tests case 0 < f < 2^-20, k = 0
;; log(1+2^-21)
(assert-eql 4.7683704451632344d-7 (log (+ 1 (scale-float 1d0 -21))))
;; Tests case 0 < f < 2^-20, k = 1
;; log(2 + 2^-20)
(assert-eql 0.6931476573969898d0 (log (+ 2(scale-float 1d0 -20))))
(assert-eql 1.3862943611198906d0 (log 4d0))
;; Tests main path, i > 0, k = 0
(assert-eql 0.3220828910287846d0
(log (kernel:make-double-float (+ #x3ff00000 #x6147a) 0)))
;; Tests main path, i > 0, k = 1
(assert-eql 0.35065625373947773d0
(log (kernel:make-double-float (+ #x3ff00000 #x6b851) 0)))
;; Tests main path, i > 0, k = -1
(assert-eql -0.3710642895311607d0
(log (kernel:make-double-float (+ #x3fe00000 #x6147a) 0)))
;; Tests main path, i < 0, k = 0
(assert-eql 0.3220821999597803d0
(log (kernel:make-double-float (+ #x3ff00000 #x61479) 0)))
;; Tests main path, i < 0, k = 1
(assert-eql 1.0152293805197257d0
(log (kernel:make-double-float (+ #x40000000 #x61479) 0)))
;; Tests main path, i < 0, k = -1
(assert-eql -0.37106498060016496d0
(log (kernel:make-double-float (+ #x3fe00000 #x61479) 0))))
(define-test log-consistency
(:tag :fdlibm)
;; |log(x) + log(1/x)| < 1.77635684e-15, x = 1.2^k, 0 <= k < 2000
;; The threshold is experimentally determined
(let ((x 1d0)
(max-value -1d0))
(declare (double-float max-value)
(type (double-float 1d0) x))
(dotimes (k 2000)
(let ((y (abs (+ (log x) (log (/ x))))))
(setf max-value (max max-value y))
(setf x (* x 1.4d0))))
(assert-true (< max-value 1.77635684d-15)))
;; |exp(log(x)) - x|/x < 5.6766649d-14, x = 1.4^k, 0 <= k < 2000
(let ((x 1d0)
(max-error 0d0))
(declare (double-float max-error)
(type (double-float 1d0) x))
(dotimes (k 2000)
(let ((y (abs (/ (- (exp (log x)) x) x))))
(setf max-error (max max-error y))
(setf x (* x 1.4d0))))
(assert-true (< max-error 5.6766649d-14)))
;; |exp(log(x)) - x|/x < 5.68410245d-14, x = 1.4^(-k), 0 <= k < 2000
(let ((x 1d0)
(max-error 0d0))
(declare (double-float max-error)
(type (double-float (0d0)) x))
(dotimes (k 2000)
(let ((y (abs (/ (- (exp (log x)) x) x))))
(setf max-error (max max-error y))
(setf x (/ x 1.4d0))))
(assert-true (< max-error 5.68410245d-14))))
(define-test sinh-basic-tests
(:tag :fdlibm)
(assert-eql +0d0 (sinh 0d0))
(assert-eql -0d0 (sinh -0d0))
;; sinh(x) = x, |x| < 2^-28
(let ((x (scale-float 1d0 -29)))
(assert-eql x (sinh x))
(assert-eql (- x) (sinh (- x))))
;; case |x| < 1
(assert-eql 0.5210953054937474d0 (sinh 0.5d0))
(assert-eql -0.5210953054937474d0 (sinh -0.5d0))
;; sinh(10*log(2)) = 1048575/2048, case |x| < 22
(let ((x (* 10 (log 2d0)))
(y (float 1048575/2048 1d0)))
(assert-eql y (sinh x))
(assert-eql (- y) (sinh (- x))))
;; sinh(10), case |x| < 22
(let ((y 11013.232874703393d0))
(assert-eql y (sinh 10d0))
(assert-eql (- y) (sinh -10d0)))
;; sinh(32*log(2)), case |x| in [22, log(maxdouble)]
(let ((x (* 32 (log 2d0)))
(y 2.1474836479999983d9))
(assert-eql y (sinh x))
(assert-eql (- y) (sinh (- x))))
;; sinh(100), case |x| in [22, log(maxdouble)]
(let ((y 1.3440585709080678d43))
(assert-eql y (sinh 100d0))
(assert-eql (- y) (sinh -100d0)))
;; sinh(710....), no overflow, case |x| in [log(maxdouble), overflowthreshold]
(let ((x 710.4758600739439d0)
(y 1.7976931348621744d308))
(assert-eql y (sinh x))
(assert-eql (- y) (sinh (- x))))
;; sinh(710.475860073944), overflow, case |x| > ovfthreshold]
(let ((x 710.475860073944d0))
(assert-error 'floating-point-overflow (sinh x))
(assert-error 'floating-point-overflow (sinh (- x))))
(assert-error 'floating-point-overflow (sinh 1000d0))
(assert-error 'floating-point-overflow (sinh -1000d0)))
(define-test tanh-basic-tests
(:tag :fdlibm)
;; case |x| < 2^-55
(let ((x (scale-float 1d0 -56)))
(assert-eql x (tanh x))
(assert-eql (- x) (tanh (- x))))
;; tanh(log(2)) = 3/5, case |x| < 1
(let ((x (log 2d0))
(y (float 3/5 1d0)))
(assert-eql y (tanh x))
(assert-eql (- y) (tanh (- x))))
;; tanh(2*log(2)) = 15/17, case |x| < 22
(let ((x (* 2 (log 2d0)))
(y (float 15/17 1d0)))
(assert-eql y (tanh x))
(assert-eql (- y) (tanh (- x))))
;; tanh(100) = 1, case |x| > 22
(assert-eql 1d0 (tanh 100d0))
(assert-eql -1d0 (tanh -100d0))
;; tanh(1d300), no overflow
(assert-eql 1d0 (tanh most-positive-double-float))
(assert-eql -1d0 (tanh (- most-positive-double-float))))
(define-test %asin-basic-tests
(:tag :fdlibm)
(let ((x (scale-float 1d0 -28))
(x0 0d0))
;; asin(x) = x for |x| < 2^-27, with inexact exception if x is not 0.
(assert-eql x (kernel:%asin x))
(assert-eql (- x) (kernel:%asin (- x)))))
(define-test %asin-exception
(:tag :fdlibm)
(let ((x (scale-float 1d0 -28))
(x0 0d0))
;; asin(x) = x for |x| < 2^-27, with inexact exception if x is not 0.
(assert-eql x (kernel:%asin x))
(assert-eql (- x) (kernel:%asin (- x)))
(with-inexact-exception-enabled
;; This must not throw an inexact exception because the result
;; is exact when the arg is 0.
(assert-eql 0d0 (kernel:%asin x0)))
(with-inexact-exception-enabled
;; This must throw an inexact exception for non-zero x even
;; though the result is exactly x.
(assert-error 'floating-point-inexact
(kernel:%asin x)))))
(define-test %cos.exceptions
(:tag :fdlibm)
;; cos(x) = 1 for |x| < 2^-27. Signal inexact unless x = 0
(let ((x (scale-float 1d0 -28))
(x0 0d0))
(with-inexact-exception-enabled
;; This must not throw an inexact exception because the result
;; is exact when the arg is 0.
(assert-eql 1d0 (kernel:%cos x0)))
(with-inexact-exception-enabled
;; This must throw an inexact exception for non-zero x even
;; though the result is exactly x.
(assert-error 'floating-point-inexact
(kernel:%cos x)))))
(define-test %sin.exceptions
(:tag :fdlibm)
;; sin(x) = x for |x| < 2^-27. Signal inexact unless x = 0
(let ((x (scale-float 1d0 -28))
(x0 0d0))
(with-inexact-exception-enabled
;; This must not throw an inexact exception because the result
;; is exact when the arg is 0.
(assert-eql 0d0 (kernel:%sin x0)))
(with-inexact-exception-enabled
;; This must throw an inexact exception for non-zero x even
;; though the result is exactly x.
(assert-error 'floating-point-inexact
(kernel:%sin x)))))
(define-test %tan.exceptions
(:tag :fdlibm)
;; tan(x) = x for |x| < 2^-28. Signal inexact unless x = 0
(let ((x (scale-float 1d0 -29))
(x0 0d0))
(with-inexact-exception-enabled
;; This must not throw an inexact exception because the result
;; is exact when the arg is 0.
(assert-eql 0d0 (kernel:%tan x0)))
(with-inexact-exception-enabled
;; This must throw an inexact exception for non-zero x even
;; though the result is exactly x.
(assert-error 'floating-point-inexact
(kernel:%tan x)))))