Raymond Toy pushed to branch master at cmucl / cmucl
Commits: a53d7ef4 by Raymond Toy at 2015-12-23T14:14:13Z Use setexception to raise the inexact exception.
o Update fdlibm.h and setexception.c to support the inexact execption. o Use this in asinh.
Tests pass.
- - - - - 0d53bc7f by Raymond Toy at 2015-12-23T14:24:49Z Merge branch 'master' into rtoy-setexception-inexact
- - - - - e655d017 by Raymond Toy at 2015-12-23T14:30:10Z Use setexception to raise the inexact exception for asin.
o Add tests for this o Use setexception for inexact in e_asin.c.
- - - - - 8b36c06e by Raymond Toy at 2015-12-23T15:48:58Z Group the inexact exception test with the exceptions tests.
- - - - - b4c91767 by Raymond Toy at 2015-12-23T15:55:19Z Use setexception to raise the inexact exception for exp.
o Add tests for this o Use setexception for inexact in e_exp.c.
- - - - - e90e91d4 by Raymond Toy at 2015-12-23T19:43:17Z Use setexception to raise the inexact exception for sinh.
- - - - - d448ca78 by Raymond Toy at 2015-12-23T19:49:21Z Use setexception to raise the inexact exception for cos.
- - - - - 71bdff74 by Raymond Toy at 2015-12-24T08:54:03Z Use setexception to raise the inexact exception for sin.
- - - - - a31150c5 by Raymond Toy at 2015-12-24T08:59:29Z Use setexception to raise the inexact exception for tan.
- - - - - ae70cdd3 by Raymond Toy at 2015-12-24T09:06:54Z Use setexception to raise the inexact exception for atan.
- - - - - 89097cd0 by Raymond Toy at 2015-12-24T09:15:27Z Use setexception to raise the inexact exception for %expm1.
- - - - - 91ff3607 by Raymond Toy at 2015-12-24T09:23:06Z Use setexception to raise the inexact exception for %log1p.
- - - - - b9e3a511 by Raymond Toy at 2015-12-24T09:29:08Z Use setexception to raise the inexact exception for tanh.
Note that original code didn't actually signal inexact probably because the compiler constant-folded one - tiny to one.
- - - - - 97bd0eaa by Raymond Toy at 2015-12-24T10:37:57Z Add WITH-FLOAT-TRAPS-ENABLED to enable specific traps.
This works like WITH-FLOAT-TRAPS-MASKED, except that the specified traps are enabled.
Use this in fdlibm to enable the inexact trap.
- - - - - 5721ddd2 by Raymond Toy at 2015-12-24T11:46:36Z Simplify WITH-FLOAT-TRAPS-MASKED and WITH-FLOAT-TRAPS-ENABLED.
Merge the body of both macros into one since they only differ in how the bits are merged with the actual mode bits.
- - - - - d2a8d5c7 by Raymond Toy at 2015-12-24T22:25:13Z ADD docstrings for WITH-FLOAT-TRAPS-MASKED and WITH-FLOAT-TRAPS-ENABLED.
- - - - - 519d5133 by Raymond Toy at 2015-12-24T22:47:24Z (setf floating-point-modes) wants (unsigned-byte 24)
When enabling traps, need to take just the low 24 bits of the arg because (setf floating-point-modes) wants an (unsigned-byte 24) argument. The logorc2 makes the result negative when enabling traps.
- - - - - 46e43aed by Raymond Toy at 2015-12-24T22:48:49Z Use correct package (EXT) for WITH-FLOAT-TRAPS-MASKED.
Also replae WITH-INXACT-EXCEPTION-ENABLED with WITH-FLOAT-TRAPS-ENABLED.
All tests still pass, as expected.
- - - - - 38e8ce5c by Raymond Toy at 2015-12-26T09:09:56Z Fix bug on sparc and clean up.
On sparc and ppc (setf vm:floating-point-modes) takes an (unsigned-byte 32) arg, so adjust the ldb byte appopriately.
Clean up code by putting the docstring into the macro.
- - - - - 6cc16b9b by Raymond Toy at 2015-12-26T09:11:40Z Regenerated du to float-traps.lisp changes.
- - - - - 55b541e5 by Raymond Toy at 2015-12-26T09:49:25Z Add shell script to run the test suite.
This makes it quite a bit easier to run the test suite instead of trying to remember exactly how to invoke it from the command line.
- - - - - 62acaf64 by Raymond Toy at 2015-12-27T10:36:59Z Merge branch 'master' into rtoy-setexception-inexact
- - - - - 90b9651b by Raymond Toy at 2015-12-27T21:02:14Z Clean up with-float-traps macro.
* Add some comments. * Change x86 (setf floating-point-modes) to accept (unsigned-byte 32). * Remove unneeded x86 conditionalization on the byte size.
- - - - - 80d7ca4e by Raymond Toy at 2015-12-28T10:26:27Z Fix compiler warnings by removing unused vars.
- - - - - 6a6908fb by Raymond Toy at 2015-12-28T10:28:57Z Remove trailing blank line.
- - - - - fc1c9daa by Raymond Toy at 2015-12-28T10:31:55Z Replace with-inexact-exception-enabled with with-float-traps-enabled.
- - - - - 5b83139e by Raymond Toy at 2015-12-28T18:37:04Z Use setexception to set inexact exception
Fix issue #12 by replacing the code with an explicit call to set the inexact exception when needed.
- - - - -
20 changed files:
- + bin/run-tests.sh - src/code/exports.lisp - src/code/float-trap.lisp - src/i18n/locale/cmucl.pot - src/lisp/e_asin.c - src/lisp/e_cosh.c - src/lisp/e_exp.c - src/lisp/e_sinh.c - src/lisp/fdlibm.h - src/lisp/k_cos.c - src/lisp/k_sin.c - src/lisp/k_tan.c - src/lisp/s_asinh.c - src/lisp/s_atan.c - src/lisp/s_expm1.c - src/lisp/s_log1p.c - src/lisp/s_scalbn.c - src/lisp/s_tanh.c - src/lisp/setexception.c - tests/fdlibm.lisp
Changes:
===================================== bin/run-tests.sh ===================================== --- /dev/null +++ b/bin/run-tests.sh @@ -0,0 +1,50 @@ +#! /bin/bash + +# Run the testsuite. +# +# By default, all the tests are run, but if additional args are given, +# then just those tests are run. + +usage() { + echo "run-tests.sh [?] [-l lisp] [tests]" + echo " -l lisp Lisp to use for the tests; defaults to lisp" + echo " -? This help message" + echo "" + echo "Run the test suite" + echo "" + echo "Any remaining args are the names of the tests to run." + echo "These are basically the file names (without extension)" + echo "in the tests/ directory." + echo "" + echo "This script expects to be run from the top level of the" + echo "cmucl source tree. That is, is should be invoked as" + echo "bin/run-tests.sh" + exit 0; +} + +LISP=lisp +while getopts "h?l:" arg +do + case $arg in + l) LISP=$OPTARG ;; + ?) usage ;; + esac +done + +# Shift out the options +shift $[$OPTIND - 1] + +if [ $# -eq 0 ]; then + # No args so run all the tests + $LISP -noinit -load tests/run-tests.lisp -eval '(cmucl-test-runner:run-all-tests)' +else + # Run selected files. Convert each file name to uppercase and append "-TESTS" + result="" + for f in $* + do + new=`echo $f | tr '[a-z]' '[A-Z]'` + result="$result ""$new-TESTS" + done + $LISP -noinit -load tests/run-tests.lisp -eval "(progn (cmucl-test-runner:load-test-files) (cmucl-test-runner:run-test $result))" +fi +
===================================== src/code/exports.lisp ===================================== --- a/src/code/exports.lisp +++ b/src/code/exports.lisp @@ -1583,7 +1583,8 @@ "FLOAT-DENORMALIZED-P" "FLOAT-INFINITY-P" "FLOAT-NAN-P" "FLOAT-TRAPPING-NAN-P" "FLOAT-SIGNALING-NAN-P" - "WITH-FLOAT-TRAPS-MASKED") + "WITH-FLOAT-TRAPS-MASKED" + "WITH-FLOAT-TRAPS-ENABLED") ;; More float extensions #+double-double (:export "LEAST-POSITIVE-NORMALIZED-DOUBLE-DOUBLE-FLOAT"
===================================== src/code/float-trap.lisp ===================================== --- a/src/code/float-trap.lisp +++ b/src/code/float-trap.lisp @@ -23,7 +23,8 @@ ) (in-package "EXTENSIONS") (export '(set-floating-point-modes get-floating-point-modes - with-float-traps-masked)) + with-float-traps-masked + with-float-traps-enabled)) (in-package "VM")
(eval-when (compile load eval) @@ -103,7 +104,7 @@
final-mode)) (defun (setf floating-point-modes) (new-mode) - (declare (type (unsigned-byte 24) new-mode)) + (declare (type (unsigned-byte 32) new-mode)) ;; Set the floating point modes for both X87 and SSE2. This ;; include the rounding control bits. (let* ((rc (ldb float-rounding-mode new-mode)) @@ -116,8 +117,8 @@ ;; is ok and would be the correct setting if we ;; ever support long-floats. (ash 3 8)))) - (setf (vm::sse2-floating-point-modes) new-mode) - (setf (vm::x87-floating-point-modes) x87-modes)) + (setf (vm::sse2-floating-point-modes) (ldb (byte 24 0) new-mode)) + (setf (vm::x87-floating-point-modes) (ldb (byte 24 0) x87-modes))) new-mode) )
@@ -364,45 +365,66 @@ (error _"SIGFPE with no exceptions currently enabled? (si-code = ~D)" code)))))))
-;;; WITH-FLOAT-TRAPS-MASKED -- Public -;;; -(defmacro with-float-traps-masked (traps &body body) - "Execute BODY with the floating point exceptions listed in TRAPS +(macrolet + ((with-float-traps (name logical-op docstring) + ;; Define macros to enable or disable floating-point + ;; exceptions. Masked exceptions and enabled exceptions only + ;; differ whether we AND in the bits or OR them, respectively. + ;; Logical-op is the operation to use. + (let ((macro-name (symbolicate "WITH-FLOAT-TRAPS-" name))) + `(progn + (defmacro ,macro-name (traps &body body) + ,docstring + (let ((traps (dpb (float-trap-mask traps) float-traps-byte 0)) + (exceptions (dpb (float-trap-mask traps) float-sticky-bits 0)) + (trap-mask (dpb (lognot (float-trap-mask traps)) + float-traps-byte #xffffffff)) + (exception-mask (dpb (lognot (vm::float-trap-mask traps)) + float-sticky-bits #xffffffff)) + ;; On ppc if we are masking the invalid trap, we need to make + ;; sure we wipe out the various individual sticky bits + ;; representing the invalid operation. Otherwise, if we + ;; enable the invalid trap later, these sticky bits will cause + ;; an exception. + #+ppc + (invalid-mask (if (member :invalid traps) + (dpb 0 + (byte 1 31) + (dpb 0 vm::float-invalid-op-2-byte + (dpb 0 vm:float-invalid-op-1-byte #xffffffff))) + #xffffffff)) + (orig-modes (gensym))) + `(let ((,orig-modes (floating-point-modes))) + (unwind-protect + (progn + (setf (floating-point-modes) + (ldb (byte 32 0) + (,',logical-op ,orig-modes ,(logand trap-mask exception-mask)))) + ,@body) + ;; Restore the original traps and exceptions. + (setf (floating-point-modes) + (logior (logand ,orig-modes ,(logior traps exceptions)) + (logand (floating-point-modes) + ,(logand trap-mask exception-mask) + #+ppc + ,invalid-mask + #+mips ,(dpb 0 float-exceptions-byte #xffffffff)))))))))))) + + ;; WITH-FLOAT-TRAPS-MASKED -- Public + (with-float-traps masked logand + _N"Execute BODY with the floating point exceptions listed in TRAPS masked (disabled). TRAPS should be a list of possible exceptions which includes :UNDERFLOW, :OVERFLOW, :INEXACT, :INVALID and :DIVIDE-BY-ZERO and on the X86 :DENORMALIZED-OPERAND. The respective accrued exceptions are cleared at the start of the body to support - their testing within, and restored on exit." - (let ((traps (dpb (float-trap-mask traps) float-traps-byte 0)) - (exceptions (dpb (float-trap-mask traps) float-sticky-bits 0)) - (trap-mask (dpb (lognot (float-trap-mask traps)) - float-traps-byte #xffffffff)) - (exception-mask (dpb (lognot (vm::float-trap-mask traps)) - float-sticky-bits #xffffffff)) - ;; On ppc if we are masking the invalid trap, we need to make - ;; sure we wipe out the various individual sticky bits - ;; representing the invalid operation. Otherwise, if we - ;; enable the invalid trap later, these sticky bits will cause - ;; an exception. - #+ppc - (invalid-mask (if (member :invalid traps) - (dpb 0 - (byte 1 31) - (dpb 0 vm::float-invalid-op-2-byte - (dpb 0 vm:float-invalid-op-1-byte #xffffffff))) - #xffffffff)) - (orig-modes (gensym))) - `(let ((,orig-modes (floating-point-modes))) - (unwind-protect - (progn - (setf (floating-point-modes) - (logand ,orig-modes ,(logand trap-mask exception-mask))) - ,@body) - ;; Restore the original traps and exceptions. - (setf (floating-point-modes) - (logior (logand ,orig-modes ,(logior traps exceptions)) - (logand (floating-point-modes) - ,(logand trap-mask exception-mask) - #+ppc - ,invalid-mask - #+mips ,(dpb 0 float-exceptions-byte #xffffffff)))))))) + their testing within, and restored on exit.") + + ;; WITH-FLOAT-TRAPS-ENABLED -- Public + (with-float-traps enabled logorc2 + _N"Execute BODY with the floating point exceptions listed in TRAPS + enabled. TRAPS should be a list of possible exceptions which + includes :UNDERFLOW, :OVERFLOW, :INEXACT, :INVALID and + :DIVIDE-BY-ZERO and on the X86 :DENORMALIZED-OPERAND. The respective + accrued exceptions are cleared at the start of the body to support + their testing within, and restored on exit.")) +
===================================== src/i18n/locale/cmucl.pot ===================================== --- a/src/i18n/locale/cmucl.pot +++ b/src/i18n/locale/cmucl.pot @@ -4797,6 +4797,16 @@ msgid "" " their testing within, and restored on exit." msgstr ""
+#: src/code/float-trap.lisp +msgid "" +"Execute BODY with the floating point exceptions listed in TRAPS\n" +" enabled. TRAPS should be a list of possible exceptions which\n" +" includes :UNDERFLOW, :OVERFLOW, :INEXACT, :INVALID and\n" +" :DIVIDE-BY-ZERO and on the X86 :DENORMALIZED-OPERAND. The respective\n" +" accrued exceptions are cleared at the start of the body to support\n" +" their testing within, and restored on exit." +msgstr "" + #: src/code/float.lisp msgid "Return true if the float X is denormalized." msgstr ""
===================================== src/lisp/e_asin.c ===================================== --- a/src/lisp/e_asin.c +++ b/src/lisp/e_asin.c @@ -50,7 +50,6 @@ static const double static double #endif one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */ -huge = 1.000e+300, pio2_hi = 1.57079632679489655800e+00, /* 0x3FF921FB, 0x54442D18 */ pio2_lo = 6.12323399573676603587e-17, /* 0x3C91A626, 0x33145C07 */ pio4_hi = 7.85398163397448278999e-01, /* 0x3FE921FB, 0x54442D18 */ @@ -89,7 +88,12 @@ qS4 = 7.70381505559019352791e-02; /* 0x3FB3B8C5, 0xB12E9282 */ return fdlibm_setexception(x, FDLIBM_INVALID); } else if (ix<0x3fe00000) { /* |x|<0.5 */ if(ix<0x3e400000) { /* if |x| < 2**-27 */ - if(huge+x>one) return x;/* return x with inexact if x!=0*/ + /* return x inexact except 0 */ + if (x != 0) { + fdlibm_setexception(x, FDLIBM_INEXACT); + } + + return x; } else t = x*x; p = t*(pS0+t*(pS1+t*(pS2+t*(pS3+t*(pS4+t*pS5)))));
===================================== src/lisp/e_cosh.c ===================================== --- a/src/lisp/e_cosh.c +++ b/src/lisp/e_cosh.c @@ -35,7 +35,7 @@ #include "fdlibm.h"
#ifdef __STDC__ -static const double one = 1.0, half=0.5, huge = 1.0e307; +static const double one = 1.0, half=0.5; #else static double one = 1.0, half=0.5, huge = 1.0e307; #endif
===================================== src/lisp/e_exp.c ===================================== --- a/src/lisp/e_exp.c +++ b/src/lisp/e_exp.c @@ -82,7 +82,6 @@ static double #endif one = 1.0, halF[2] = {0.5,-0.5,}, -huge = 1.0e+300, twom1000= 9.33263618503218878990e-302, /* 2**-1000=0x01700000,0*/ o_threshold= 7.09782712893383973096e+02, /* 0x40862E42, 0xFEFA39EF */ u_threshold= -7.45133219101941108420e+02, /* 0xc0874910, 0xD52D3051 */ @@ -161,7 +160,12 @@ P5 = 4.13813679705723846039e-08; /* 0x3E663769, 0x72BEA4D0 */ x = hi - lo; } else if(hx < 0x3e300000) { /* when |x|<2**-28 */ - if(huge+x>one) return one+x;/* trigger inexact */ + /* return x inexact except 0 */ + if (x != 0) { + fdlibm_setexception(x, FDLIBM_INEXACT); + } + + return one + x; } else k = 0;
===================================== src/lisp/e_sinh.c ===================================== --- a/src/lisp/e_sinh.c +++ b/src/lisp/e_sinh.c @@ -32,7 +32,7 @@ #include "fdlibm.h"
#ifdef __STDC__ -static const double one = 1.0, shuge = 1.0e307; +static const double one = 1.0; #else static double one = 1.0, shuge = 1.0e307; #endif @@ -67,8 +67,14 @@ static double one = 1.0, shuge = 1.0e307; if (jx<0) h = -h; /* |x| in [0,22], return sign(x)*0.5*(E+E/(E+1))) */ if (ix < 0x40360000) { /* |x|<22 */ - if (ix<0x3e300000) /* |x|<2**-28 */ - if(shuge+x>one) return x;/* sinh(tiny) = tiny with inexact */ + if (ix<0x3e300000) { /* |x|<2**-28 */ + /* sinh(tiny) = tiny with inexact */ + if (x != 0) { + fdlibm_setexception(x, FDLIBM_INEXACT); + } + + return x; + } t = fdlibm_expm1(fabs(x)); if(ix<0x3ff00000) return h*(2.0*t-t*t/(t+one)); return h*(t+t/(t+one));
===================================== src/lisp/fdlibm.h ===================================== --- a/src/lisp/fdlibm.h +++ b/src/lisp/fdlibm.h @@ -61,7 +61,8 @@ enum FDLIBM_EXCEPTION { FDLIBM_DIVIDE_BY_ZERO, FDLIBM_UNDERFLOW, FDLIBM_OVERFLOW, - FDLIBM_INVALID + FDLIBM_INVALID, + FDLIBM_INEXACT };
extern double fdlibm_setexception(double x, enum FDLIBM_EXCEPTION);
===================================== src/lisp/k_cos.c ===================================== --- a/src/lisp/k_cos.c +++ b/src/lisp/k_cos.c @@ -75,7 +75,12 @@ C6 = -1.13596475577881948265e-11; /* 0xBDA8FAE9, 0xBE8838D4 */ ux.d = x; ix = ux.i[HIWORD]&0x7fffffff; /* ix = |x|'s high word*/ if(ix<0x3e400000) { /* if x < 2**27 */ - if(((int)x)==0) return one; /* generate inexact */ + /* return 1 with inexact unless x == 0 */ + if (x != 0) { + fdlibm_setexception(x, FDLIBM_INEXACT); + } + + return one; } z = x*x; r = z*(C1+z*(C2+z*(C3+z*(C4+z*(C5+z*C6)))));
===================================== src/lisp/k_sin.c ===================================== --- a/src/lisp/k_sin.c +++ b/src/lisp/k_sin.c @@ -67,8 +67,14 @@ S6 = 1.58969099521155010221e-10; /* 0x3DE5D93A, 0x5ACFD57C */
ux.d = x; ix = ux.i[HIWORD]&0x7fffffff; /* high word of x */ - if(ix<0x3e400000) /* |x| < 2**-27 */ - {if((int)x==0) return x;} /* generate inexact */ + if(ix<0x3e400000) { /* |x| < 2**-27 */ + /* return x inexact except 0 */ + if (x != 0) { + fdlibm_setexception(x, FDLIBM_INEXACT); + } + + return x; + } z = x*x; v = z*x; r = S2+z*(S3+z*(S4+z*(S5+z*S6)));
===================================== src/lisp/k_tan.c ===================================== --- a/src/lisp/k_tan.c +++ b/src/lisp/k_tan.c @@ -78,31 +78,34 @@ __kernel_tan(double x, double y, int iy) { hx = ux.i[HIWORD]; /* high word of x */ ix = hx & 0x7fffffff; /* high word of |x| */ if (ix < 0x3e300000) { /* x < 2**-28 */ - if ((int) x == 0) { /* generate inexact */ - if (((ix | ux.i[LOWORD]) | (iy + 1)) == 0) - return one / fabs(x); - else { - if (iy == 1) - return x; - else { /* compute -1 / (x+y) carefully */ - double a, t; + /* return x inexact except 0 */ + if (x != 0) { + fdlibm_setexception(x, FDLIBM_INEXACT); + }
- z = w = x + y; - uz.d = z; - uz.i[LOWORD] = 0; - z = ux.d; + if (((ix | ux.i[LOWORD]) | (iy + 1)) == 0) + return one / fabs(x); + else { + if (iy == 1) + return x; + else { /* compute -1 / (x+y) carefully */ + double a, t; + + z = w = x + y; + uz.d = z; + uz.i[LOWORD] = 0; + z = ux.d;
- v = y - (z - x); - t = a = -one / w; - uz.d = t; - uz.i[LOWORD] = 0; - t = uz.d; + v = y - (z - x); + t = a = -one / w; + uz.d = t; + uz.i[LOWORD] = 0; + t = uz.d;
- s = one + t * z; - return t + a * (s + t * v); - } - } + s = one + t * z; + return t + a * (s + t * v); } + } } if (ix >= 0x3FE59428) { /* |x| >= 0.6744 */ if (hx < 0) {
===================================== src/lisp/s_asinh.c ===================================== --- a/src/lisp/s_asinh.c +++ b/src/lisp/s_asinh.c @@ -32,8 +32,7 @@ static const double static double #endif one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */ -ln2 = 6.93147180559945286227e-01, /* 0x3FE62E42, 0xFEFA39EF */ -huge= 1.00000000000000000000e+300; +ln2 = 6.93147180559945286227e-01; /* 0x3FE62E42, 0xFEFA39EF */
#ifdef __STDC__ double fdlibm_asinh(double x) @@ -59,7 +58,12 @@ huge= 1.00000000000000000000e+300; }
if(ix< 0x3e300000) { /* |x|<2**-28 */ - if(huge+x>one) return x; /* return x inexact except 0 */ + /* return x inexact except 0 */ + if (x != 0) { + fdlibm_setexception(x, FDLIBM_INEXACT); + } + + return x; } if(ix>0x41b00000) { /* |x| > 2**28 */ w = __ieee754_log(fabs(x))+ln2;
===================================== src/lisp/s_atan.c ===================================== --- a/src/lisp/s_atan.c +++ b/src/lisp/s_atan.c @@ -79,8 +79,7 @@ static double aT[] = { #else static double #endif -one = 1.0, -huge = 1.0e300; +one = 1.0;
#ifdef __STDC__ double fdlibm_atan(double x) @@ -104,7 +103,12 @@ huge = 1.0e300; else return -atanhi[3]-atanlo[3]; } if (ix < 0x3fdc0000) { /* |x| < 0.4375 */ if (ix < 0x3e200000) { /* |x| < 2^-29 */ - if(huge+x>one) return x; /* raise inexact */ + /* return x inexact except 0 */ + if (x != 0) { + fdlibm_setexception(x, FDLIBM_INEXACT); + } + + return x; } id = -1; } else {
===================================== src/lisp/s_expm1.c ===================================== --- a/src/lisp/s_expm1.c +++ b/src/lisp/s_expm1.c @@ -160,8 +160,8 @@ Q5 = -2.01099218183624371326e-07; /* BE8AFDB7 6E09C32D */
} if(xsb!=0) { /* x < -56*ln2, return -1.0 with inexact */ - if(x+tiny<0.0) /* raise inexact */ - return tiny-one; /* return -1 */ + fdlibm_setexception(x, FDLIBM_INEXACT); + return tiny - one; } }
===================================== src/lisp/s_log1p.c ===================================== --- a/src/lisp/s_log1p.c +++ b/src/lisp/s_log1p.c @@ -85,7 +85,6 @@ static double #endif ln2_hi = 6.93147180369123816490e-01, /* 3fe62e42 fee00000 */ ln2_lo = 1.90821492927058770002e-10, /* 3dea39ef 35793c76 */ -two54 = 1.80143985094819840000e+16, /* 43500000 00000000 */ Lp1 = 6.666666666666735130e-01, /* 3FE55555 55555593 */ Lp2 = 3.999999999940941908e-01, /* 3FD99999 9997FA04 */ Lp3 = 2.857142874366239149e-01, /* 3FD24924 94229359 */ @@ -123,9 +122,14 @@ static double zero = 0.0; } } if(ax<0x3e200000) { /* |x| < 2**-29 */ - if(two54+x>zero /* raise inexact */ - &&ax<0x3c900000) /* |x| < 2**-54 */ + if (ax < 0x3c900000) { /* |x| < 2**-54 */ + /* return x inexact except 0 */ + if (x != 0) { + fdlibm_setexception(x, FDLIBM_INEXACT); + } + return x; + } else return x - x*x*0.5; }
===================================== src/lisp/s_scalbn.c ===================================== --- a/src/lisp/s_scalbn.c +++ b/src/lisp/s_scalbn.c @@ -26,9 +26,7 @@ static const double static double #endif two54 = 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */ -twom54 = 5.55111512312578270212e-17, /* 0x3C900000, 0x00000000 */ -huge = 1.0e+300, -tiny = 1.0e-300; +twom54 = 5.55111512312578270212e-17; /* 0x3C900000, 0x00000000 */
#ifdef __STDC__ double fdlibm_scalbn (double x, int n)
===================================== src/lisp/s_tanh.c ===================================== --- a/src/lisp/s_tanh.c +++ b/src/lisp/s_tanh.c @@ -78,7 +78,9 @@ static double one=1.0, two=2.0, tiny = 1.0e-300; } /* |x| > 22, return +-1 */ } else { - z = one - tiny; /* raised inexact flag */ + /* Always raise inexact flag */ + fdlibm_setexception(x, FDLIBM_INEXACT); + z = one - tiny; } return (jx>=0)? z: -z; }
===================================== src/lisp/setexception.c ===================================== --- a/src/lisp/setexception.c +++ b/src/lisp/setexception.c @@ -88,6 +88,10 @@ fdlibm_setexception(double x, enum FDLIBM_EXCEPTION type)
break; } + case FDLIBM_INEXACT: + feraiseexcept(FE_INEXACT); + ret = x; + break; default: /* Shouldn't happen! */ ret = 0.0;
===================================== tests/fdlibm.lisp ===================================== --- a/tests/fdlibm.lisp +++ b/tests/fdlibm.lisp @@ -6,7 +6,7 @@ (in-package "FDLIBM-TESTS")
(defparameter *qnan* - (kernel::with-float-traps-masked (:invalid) + (ext:with-float-traps-masked (:invalid) (* 0 ext:double-float-positive-infinity)) "Some randon quiet MaN value")
@@ -14,15 +14,6 @@ (kernel:make-double-float #x7ff00000 1) "A randon signaling MaN value")
-(defmacro with-inexact-exception-enabled (&body body) - (let ((old-modes (gensym "OLD-MODES-"))) - `(let ((,old-modes (ext:get-floating-point-modes))) - (unwind-protect - (progn - (ext:set-floating-point-modes :traps '(:inexact)) - ,@body) - (apply 'ext:set-floating-point-modes ,old-modes))))) - (define-test %cosh.exceptions (:tag :fdlibm) (assert-error 'floating-point-overflow @@ -34,7 +25,7 @@ (assert-true (ext:float-nan-p (kernel:%cosh *qnan*)))
;; Same, but with overflow's masked - (kernel::with-float-traps-masked (:overflow) + (ext:with-float-traps-masked (:overflow) (assert-equal ext:double-float-positive-infinity (kernel:%cosh 1000d0)) (assert-equal ext:double-float-positive-infinity @@ -44,7 +35,7 @@ (assert-equal ext:double-float-positive-infinity (kernel:%cosh ext:double-float-negative-infinity))) ;; Test NaN - (kernel::with-float-traps-masked (:invalid) + (ext:with-float-traps-masked (:invalid) (assert-true (ext:float-nan-p (kernel:%cosh *snan*)))))
(define-test %sinh.exceptions @@ -57,7 +48,7 @@ (kernel:%sinh *snan*)) (assert-true (ext:float-nan-p (kernel:%sinh *qnan*))) ;; Same, but with overflow's masked - (kernel::with-float-traps-masked (:overflow) + (ext:with-float-traps-masked (:overflow) (assert-equal ext:double-float-positive-infinity (kernel:%sinh 1000d0)) (assert-equal ext:double-float-negative-infinity @@ -67,17 +58,35 @@ (assert-equal ext:double-float-negative-infinity (kernel:%sinh ext:double-float-negative-infinity))) ;; Test NaN - (kernel::with-float-traps-masked (:invalid) - (assert-true (ext:float-nan-p (kernel:%sinh *qnan*))))) - + (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)) + (ext:with-float-traps-enabled (:inexact) + ;; This must not throw an inexact exception because the result + ;; is exact when the arg is 0. + (assert-eql 0d0 (kernel:%sinh x0))) + (ext:with-float-traps-enabled (:inexact) + ;; 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*)) - (kernel::with-float-traps-masked (:invalid) - (assert-true (ext:float-nan-p (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)) + (ext:with-float-traps-enabled (:inexact) + ;; 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) @@ -85,10 +94,10 @@ (kernel:%acosh ext:double-float-positive-infinity)) (assert-error 'floating-point-invalid-operation (kernel:%acosh 0d0)) - (kernel::with-float-traps-masked (:overflow) + (ext:with-float-traps-masked (:overflow) (assert-equal ext:double-float-positive-infinity (kernel:%acosh ext:double-float-positive-infinity))) - (kernel::with-float-traps-masked (:invalid) + (ext:with-float-traps-masked (:invalid) (assert-true (ext:float-nan-p (kernel:%acosh 0d0)))))
(define-test %asinh.exceptions @@ -100,13 +109,24 @@ (assert-error 'floating-point-overflow (kernel:%asinh ext:double-float-negative-infinity)) (assert-true (ext:float-nan-p (kernel:%asinh *qnan*))) - (kernel::with-float-traps-masked (:overflow) + (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))) - (kernel::with-float-traps-masked (:invalid) - (assert-true (ext:float-nan-p (kernel:%asinh *snan*))))) + (ext:with-float-traps-masked (:invalid) + (assert-true (ext:float-nan-p (kernel:%asinh *snan*)))) + (let ((x (scale-float 1d0 -29)) + (x0 0d0)) + (ext:with-float-traps-enabled (:inexact) + ;; This must not throw an inexact exception because the result + ;; is exact when the arg is 0. + (assert-eql 0d0 (asinh x0))) + (ext:with-float-traps-enabled (:inexact) + ;; 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) @@ -118,10 +138,10 @@ (kernel:%atanh 1d0)) (assert-error 'division-by-zero (kernel:%atanh -1d0)) - (kernel::with-float-traps-masked (:invalid) + (ext:with-float-traps-masked (:invalid) (assert-true (ext:float-nan-p (kernel:%atanh 2d0))) (assert-true (ext:float-nan-p (kernel:%atanh -2d0)))) - (kernel::with-float-traps-masked (:divide-by-zero) + (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 @@ -136,12 +156,17 @@ (assert-error 'floating-point-invalid-operation (kernel:%expm1 *snan*)) (assert-true (ext:float-nan-p (kernel:%expm1 *qnan*))) - (kernel::with-float-traps-masked (:overflow) + (ext:with-float-traps-masked (:overflow) (assert-equal ext:double-float-positive-infinity (kernel:%expm1 709.8d0)) ) - (kernel::with-float-traps-masked (:invalid) - (assert-true (ext::float-nan-p (kernel:%expm1 *snan*))))) + (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)))) + (ext:with-float-traps-enabled (:inexact) + (assert-error 'floating-point-inexact + (kernel:%expm1 x)))))
(define-test %log1p.exceptions (:tag :fdlibm) @@ -150,11 +175,23 @@ (assert-error 'floating-point-overflow (kernel:%log1p -1d0)) (assert-true (ext:float-nan-p (kernel:%log1p *qnan*))) - (kernel::with-float-traps-masked (:overflow) + (ext:with-float-traps-masked (:overflow) (assert-equal ext:double-float-negative-infinity (kernel:%log1p -1d0))) - (kernel::with-float-traps-masked (:invalid) - (assert-true (ext:float-nan-p (kernel:%log1p *snan*))))) + (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)) + (ext:with-float-traps-enabled (:inexact) + ;; This must not throw an inexact exception because the result + ;; is exact when the arg is 0. + (assert-eql 0d0 (kernel:%log1p x0))) + (ext:with-float-traps-enabled (:inexact) + ;; 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) @@ -167,7 +204,7 @@ (kernel:%exp ext:double-float-positive-infinity)) (assert-equal 0d0 (kernel:%exp -1000d0)) - (kernel::with-float-traps-masked (:overflow) + (ext:with-float-traps-masked (:overflow) (assert-equal ext:double-float-positive-infinity (kernel:%exp 710d0))) (let ((modes (ext:get-floating-point-modes))) @@ -176,7 +213,19 @@ (ext:set-floating-point-modes :traps '(:underflow)) (assert-error 'floating-point-underflow (kernel:%exp -1000d0))) - (apply #'ext:set-floating-point-modes modes)))) + (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 + (ext:with-float-traps-enabled (:inexact) + ;; This must not throw an inexact exception because the result + ;; is exact when the arg is 0. + (assert-eql 1d0 (kernel:%exp x0))) + (ext:with-float-traps-enabled (:inexact) + ;; 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) @@ -189,12 +238,12 @@ (assert-error 'floating-point-invalid-operation (kernel:%log *snan*)) (assert-true (ext:float-nan-p (kernel:%log *qnan*))) - (kernel::with-float-traps-masked (:divide-by-zero) + (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))) - (kernel::with-float-traps-masked (:invalid) + (ext:with-float-traps-masked (:invalid) (assert-true (ext:float-nan-p (kernel:%log -1d0))) (assert-true (ext:float-nan-p (kernel:%log *snan*)))))
@@ -204,7 +253,7 @@ (kernel:%acos 2d0)) (assert-error 'floating-point-invalid-operation (kernel:%acos -2d0)) - (kernel::with-float-traps-masked (:invalid) + (ext:with-float-traps-masked (:invalid) (assert-true (ext:float-nan-p (kernel:%acos 2d0))) (assert-true (ext:float-nan-p (kernel:%acos -2d0)))))
@@ -214,7 +263,7 @@ (kernel:%asin 2d0)) (assert-error 'floating-point-invalid-operation (kernel:%asin -2d0)) - (kernel::with-float-traps-masked (:invalid) + (ext:with-float-traps-masked (:invalid) (assert-true (ext:float-nan-p (kernel:%asin 2d0))) (assert-true (ext:float-nan-p (kernel:%asin -2d0)))))
@@ -223,8 +272,20 @@ (assert-error 'floating-point-invalid-operation (kernel:%atan *snan*)) (assert-true (ext:float-nan-p (kernel:%atan *qnan*))) - (kernel::with-float-traps-masked (:invalid) - (assert-true (ext:float-nan-p (kernel:%atan *snan*))))) + (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)) + (ext:with-float-traps-enabled (:inexact) + ;; This must not throw an inexact exception because the result + ;; is exact when the arg is 0. + (assert-eql 0d0 (kernel:%atan x0))) + (ext:with-float-traps-enabled (:inexact) + ;; 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) @@ -239,12 +300,12 @@ (assert-true (ext:float-nan-p (kernel:%log10 *qnan*))) (assert-equal ext:double-float-positive-infinity (kernel:%log10 ext:double-float-positive-infinity)) - (kernel::with-float-traps-masked (:divide-by-zero) + (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))) - (kernel::with-float-traps-masked (:invalid) + (ext:with-float-traps-masked (:invalid) (assert-true (ext:float-nan-p (kernel:%log10 -1d0)))))
(define-test %scalbn.exceptions @@ -268,7 +329,7 @@ (kernel:%scalbn most-positive-double-float 2)) (assert-error 'floating-point-overflow (kernel:%scalbn most-negative-double-float 2)) - (kernel::with-float-traps-masked (:overflow) + (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 @@ -298,16 +359,7 @@ (x0 0d0)) ;; asinh(x) = x for x < 2^-28 (assert-eql x (asinh x)) - (assert-eql (- x) (asinh (- 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 (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)))) + (assert-eql (- x) (asinh (- x)))) (let ((x (scale-float 1d0 -28))) ;; Case 2 > |x| >= 2^-28 (assert-eql 3.725290298461914d-9 (asinh x)) @@ -556,4 +608,74 @@ (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)))) \ No newline at end of file + (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))) + (ext:with-float-traps-enabled (:inexact) + ;; This must not throw an inexact exception because the result + ;; is exact when the arg is 0. + (assert-eql 0d0 (kernel:%asin x0))) + (ext:with-float-traps-enabled (:inexact) + ;; 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)) + (ext:with-float-traps-enabled (:inexact) + ;; This must not throw an inexact exception because the result + ;; is exact when the arg is 0. + (assert-eql 1d0 (kernel:%cos x0))) + (ext:with-float-traps-enabled (:inexact) + ;; 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)) + (ext:with-float-traps-enabled (:inexact) + ;; This must not throw an inexact exception because the result + ;; is exact when the arg is 0. + (assert-eql 0d0 (kernel:%sin x0))) + (ext:with-float-traps-enabled (:inexact) + ;; 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)) + (ext:with-float-traps-enabled (:inexact) + ;; This must not throw an inexact exception because the result + ;; is exact when the arg is 0. + (assert-eql 0d0 (kernel:%tan x0))) + (ext:with-float-traps-enabled (:inexact) + ;; 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)))))
View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/compare/0fc0061b5302f516f051a5996...
-
Raymond Toy