Merge branch 'infinity-NaN-numeric-functions' into 'develop'

Make all numeric functions return sensible values for infinity/NaN

Closes #477

See merge request embeddable-common-lisp/ecl!160
This commit is contained in:
Daniel Kochmański 2019-08-16 18:52:23 +00:00
commit df339485eb
17 changed files with 767 additions and 69 deletions

View file

@ -64,17 +64,16 @@ si_float_infinity_p(cl_object x)
bool
ecl_float_nan_p(cl_object x)
{
return !ecl_number_equalp(x,x);
/* switch (ecl_t_of(x)) { */
/* case t_singlefloat: */
/* return !isnan(ecl_single_float(x)); */
/* case t_doublefloat: */
/* return !isnan(ecl_double_float(x)); */
/* case t_longfloat: */
/* return !isnan(ecl_long_float(x)); */
/* default: */
/* return 0; */
/* } */
switch (ecl_t_of(x)) {
case t_singlefloat:
return isnan(ecl_single_float(x));
case t_doublefloat:
return isnan(ecl_double_float(x));
case t_longfloat:
return isnan(ecl_long_float(x));
default:
return 0;
}
}
bool
@ -82,11 +81,11 @@ ecl_float_infinity_p(cl_object x)
{
switch (ecl_t_of(x)) {
case t_singlefloat:
return !isfinite(ecl_single_float(x));
return isinf(ecl_single_float(x));
case t_doublefloat:
return !isfinite(ecl_double_float(x));
return isinf(ecl_double_float(x));
case t_longfloat:
return !isfinite(ecl_long_float(x));
return isinf(ecl_long_float(x));
default:
return 0;
}

View file

@ -539,9 +539,9 @@ ecl_make_complex(cl_object r, cl_object i)
if (!ECL_REAL_TYPE_P(ti)) { ecl_type_error(@'complex', "imaginary part", i, @'real'); }
switch((tr > ti) ? tr : ti) {
#ifdef ECL_COMPLEX_FLOAT
case t_longfloat: return ecl_make_clfloat(ecl_to_long_double(r) + I * ecl_to_long_double(i));
case t_doublefloat: return ecl_make_cdfloat(ecl_to_double(r) + I * ecl_to_double(i));
case t_singlefloat: return ecl_make_csfloat(ecl_to_float(r) + I * ecl_to_float(i));
case t_longfloat: return ecl_make_clfloat(CMPLXL(ecl_to_long_double(r), ecl_to_long_double(i)));
case t_doublefloat: return ecl_make_cdfloat(CMPLX(ecl_to_double(r), ecl_to_double(i)));
case t_singlefloat: return ecl_make_csfloat(CMPLXF(ecl_to_float(r), ecl_to_float(i)));
#else
case t_singlefloat:
c = ecl_alloc_object(t_complex);
@ -599,17 +599,17 @@ ecl_make_complex_float(cl_object r, cl_object i)
case t_singlefloat:
if (ti != tr) { ecl_type_error(@'si::complex-float',"imag part", i, @'single-float'); }
result = ecl_alloc_object(t_csfloat);
ecl_csfloat(result) = ecl_single_float(r) + ecl_single_float(i) * I;
ecl_csfloat(result) = CMPLXF(ecl_single_float(r), ecl_single_float(i));
break;
case t_doublefloat:
if (ti != tr) { ecl_type_error(@'si::complex-float',"imag part", i, @'double-float'); }
result = ecl_alloc_object(t_cdfloat);
ecl_cdfloat(result) = ecl_double_float(r) + ecl_double_float(i) * I;
ecl_cdfloat(result) = CMPLX(ecl_double_float(r), ecl_double_float(i));
break;
case t_longfloat:
if (ti != tr) { ecl_type_error(@'si::complex-float',"imag part", i, @'long-float'); }
result = ecl_alloc_object(t_clfloat);
ecl_clfloat(result) = ecl_long_float(r) + ecl_long_float(i) * I;
ecl_clfloat(result) = CMPLXL(ecl_long_float(r), ecl_long_float(i));
break;
default:
ecl_type_error(@'si::complex-float',"real part", r, @'float');

View file

@ -21,7 +21,7 @@
@(defun ceiling (x &optional (y OBJNULL))
@
if (narg == 1)
return ecl_ceiling1(x);
return ecl_ceiling1(x);
else
return ecl_ceiling2(x, y);
@)

View file

@ -24,7 +24,7 @@
@(defun floor (x &optional (y OBJNULL))
@
if (narg == 1)
return ecl_floor1(x);
return ecl_floor1(x);
else
return ecl_floor2(x, y);
@)

View file

@ -24,7 +24,11 @@
}
} else do {
cl_object numi = ecl_va_arg(nums);
if (ecl_number_compare(max, numi) < 0)
if (ecl_lower(max, numi)
#ifdef ECL_IEEE_FP
|| ecl_float_nan_p(max)
#endif
)
max = numi;
} while (--narg);
@(return max);
@ -40,7 +44,11 @@
}
} else do {
cl_object numi = ecl_va_arg(nums);
if (ecl_number_compare(min, numi) > 0)
if (ecl_greater(min, numi)
#ifdef ECL_IEEE_FP
|| ecl_float_nan_p(min)
#endif
)
min = numi;
} while (--narg);
@(return min);

View file

@ -159,9 +159,18 @@ monotonic(int s, int t, int narg, ecl_va_list nums)
FEwrong_type_nth_arg(@[<], 1, c, @[real]);
}
/* INV: type check occurs in ecl_number_compare() */
for (c = ecl_va_arg(nums); --narg; c = d) {
c = ecl_va_arg(nums);
#ifdef ECL_IEEE_FP
if (ecl_float_nan_p(c))
return1(ECL_NIL);
#endif
for (; --narg; c = d) {
d = ecl_va_arg(nums);
if (s*ecl_number_compare(d, c) < t)
if (
#ifdef ECL_IEEE_FP
ecl_float_nan_p(d) ||
#endif
s*ecl_number_compare(d, c) < t)
return1(ECL_NIL);
}
return1(ECL_T);
@ -173,7 +182,7 @@ monotonic(int s, int t, int narg, ecl_va_list nums)
ecl_va_end(nums); \
return result; }
cl_object @<= MONOTONIC( 1, 0);
cl_object @>= MONOTONIC(-1, 0);
cl_object @< MONOTONIC( 1, 1);
cl_object @> MONOTONIC(-1, 1);
cl_object @<= MONOTONIC( 1, 0)
cl_object @>= MONOTONIC(-1, 0)
cl_object @< MONOTONIC( 1, 1)
cl_object @> MONOTONIC(-1, 1)

View file

@ -26,7 +26,7 @@
@(defun round (x &optional (y OBJNULL))
@
if (narg == 1)
return ecl_round1(x);
return ecl_round1(x);
else
return ecl_round2(x, y);
@)
@ -53,6 +53,8 @@ round_double(double d)
}
}
return q;
} else if (isnan(d)) {
return d;
} else {
return -round_double(-d);
}
@ -70,6 +72,8 @@ round_long_double(long double d)
}
}
return q;
} else if (isnan(d)) {
return d;
} else {
return -round_long_double(-d);
}

View file

@ -79,7 +79,7 @@ ecl_truncate2(cl_object x, cl_object y)
@(defun truncate (x &optional (y OBJNULL))
@
if (narg == 1)
return ecl_truncate1(x);
return ecl_truncate1(x);
else
return ecl_truncate2(x, y);
@)

View file

@ -258,7 +258,7 @@ cl_eq(cl_object x, cl_object y)
#define FLOAT_EQL(name, type) \
static bool name(type a, type b) { \
if (a == b) return signbit(a) == signbit(b); \
if (isnan(a) || isnan(b)) return !memcmp(&a, &b, sizeof(type)); \
if (isnan(a) || isnan(b)) return isnan(a) && isnan(b); \
return 0; \
}
#endif

View file

@ -385,30 +385,32 @@
(def-inline /= :always (t t) :bool "!ecl_number_equalp(#0,#1)")
(def-inline /= :always (fixnum-float fixnum-float) :bool "(#0)!=(#1)")
(def-inline < :always (t t) :bool "ecl_number_compare(#0,#1)<0")
(def-inline < :always (t t) :bool "ecl_lower(#0,#1)")
(def-inline < :always (fixnum-float fixnum-float) :bool "(#0)<(#1)")
(def-inline < :always (fixnum-float fixnum-float fixnum-float) :bool
"@012;((#0)<(#1) && (#1)<(#2))")
(def-inline > :always (t t) :bool "ecl_number_compare(#0,#1)>0")
(def-inline > :always (t t) :bool "ecl_greater(#0,#1)")
(def-inline > :always (fixnum-float fixnum-float) :bool "(#0)>(#1)")
(def-inline > :always (fixnum-float fixnum-float fixnum-float) :bool
"@012;((#0)>(#1) && (#1)>(#2))")
(def-inline <= :always (t t) :bool "ecl_number_compare(#0,#1)<=0")
(def-inline <= :always (t t) :bool "ecl_lowereq(#0,#1)")
(def-inline <= :always (fixnum-float fixnum-float) :bool "(#0)<=(#1)")
(def-inline <= :always (fixnum-float fixnum-float fixnum-float) :bool
"@012;((#0)<=(#1) && (#1)<=(#2))")
(def-inline >= :always (t t) :bool "ecl_number_compare(#0,#1)>=0")
(def-inline >= :always (t t) :bool "ecl_greatereq(#0,#1)")
(def-inline >= :always (fixnum-float fixnum-float) :bool "(#0)>=(#1)")
(def-inline >= :always (fixnum-float fixnum-float fixnum-float) :bool
"@012;((#0)>=(#1) && (#1)>=(#2))")
(def-inline max :always (t t) t "@01;(ecl_number_compare(#0,#1)>=0?#0:#1)")
#+ieee-floating-point (def-inline max :always (t t) t "@01;((ecl_float_nan_p(#1) || ecl_greatereq(#0,#1))?#0:#1)")
#-ieee-floating-point (def-inline max :always (t t) t "@01;(ecl_greatereq(#0,#1)?#0:#1)")
(def-inline max :always (fixnum fixnum) :fixnum "@01;(#0)>=(#1)?#0:#1")
(def-inline min :always (t t) t "@01;(ecl_number_compare(#0,#1)<=0?#0:#1)")
#+ieee-floating-point (def-inline min :always (t t) t "@01;((ecl_float_nan_p(#1) || ecl_lowereq(#0,#1))?#0:#1)")
#-ieee-floating-point (def-inline min :always (t t) t "@01;(ecl_lowereq(#0,#1)?#0:#1)")
(def-inline min :always (fixnum fixnum) :fixnum "@01;(#0)<=(#1)?#0:#1")
;; file num_log.d

View file

@ -5,6 +5,7 @@
* Numbers - Numeric types::
* Numbers - Floating point exceptions::
* Numbers - Random-States::
* Numbers - Infinity and Not a Number::
* Numbers - Dictionary::
* Numbers - C Reference::
@end menu
@ -85,8 +86,8 @@ a generalized boolean
If @var{condition} is @code{last}, @var{flag} is ignored and the
currently enabled floating point exceptions are returned in an
implementation depended format (currently an integer). Otherwise,
@var{flag} determines whether current thread will signal a floating
point exception for the conditions passed in @var{condition}.
@var{flag} determines whether the current thread will signal a
floating point exception for the conditions passed in @var{condition}.
@var{condition} can be either a symbol denoting a single condition,
@code{t} for all conditions that are enabled by default or a value
obtained from an earlier call to @code{ext:trap-fpe} with @code{last}.
@ -106,6 +107,44 @@ random seed (an integer), an array of random bytes (mainly used for
reading back printed random-state) and another random-state (syntactic
sugar for copying the random-state).
@node Numbers - Infinity and Not a Number
@subsection Infinity and Not a Number
The @ansi{} standard does not specify the behaviour of numeric
functions for infinite or not number valued floating point numbers. If
ECL is configured to support these special values (see the
@code{--with-ieee-fp} configure option) and floating point exceptions
are disabled, numeric functions generally return the same value as the
corresponding C function. This means, that the output will be a NaN
for a NaN input, and the ``mathematically correct'' value (which may
be NaN, e.g. for ∞-∞) for an infinite real input. For complex floats,
however, the return value of a numeric function called with a complex
number for which the real or imaginary part is infinite, is
undefined@footnote{The main reason for this is that some numeric
functions for C complex numbers return mathematically incorrect
values, for example sinh(i*∞) returns i*NaN instead of the
mathematically correct i*∞. Keeping this consistent with our own
implementation of complex arithmetic that is used when C complex
numbers are not available would require to much work. Furthermore,
complex arithmetic with infinities is unreliable anyway, since it
quickly leads to NaN values (consider i*∞ = (0+i*1)*(∞+i*0) = NaN+i*∞;
even this simple example is already mathematically incorrect).}.
For other functions dealing with numbers, we adopt the following
behaviour:
@table @asis
@item Comparison functions
All numeric comparisons with @code{=,<,<=,>,>=} involving NaN return
false. Comparing two NaNs of the same type with @code{eql} returns
true.
@item @code{min/max}
NaN values are ignored, i.e. the maximum/minimum is taken only over
the number valued parameters.
@item Rounding functions
All rounding functions signal an arithmetic-error if any of the given
parameters are not number valued or infinite.
@end table
@node Numbers - Dictionary
@subsection Dictionary

View file

@ -1194,10 +1194,17 @@ extern ECL_API cl_object cl_min _ECL_ARGS((cl_narg narg, cl_object min, ...));
extern ECL_API int ecl_number_equalp(cl_object x, cl_object y);
extern ECL_API int ecl_number_compare(cl_object x, cl_object y);
#ifdef ECL_IEEE_FP
#define ecl_lowereq(x,y) (!ecl_float_nan_p(x) && !ecl_float_nan_p(y) && ecl_number_compare((x),(y)) <= 0)
#define ecl_greatereq(x,y) (!ecl_float_nan_p(x) && !ecl_float_nan_p(y) && ecl_number_compare((x),(y)) >= 0)
#define ecl_lower(x,y) (!ecl_float_nan_p(x) && !ecl_float_nan_p(y) && ecl_number_compare((x),(y)) < 0)
#define ecl_greater(x,y) (!ecl_float_nan_p(x) && !ecl_float_nan_p(y) && ecl_number_compare((x),(y)) > 0)
#else
#define ecl_lowereq(x,y) (ecl_number_compare((x),(y)) <= 0)
#define ecl_greatereq(x,y) (ecl_number_compare((x),(y)) >= 0)
#define ecl_lower(x,y) (ecl_number_compare((x),(y)) < 0)
#define ecl_greater(x,y) (ecl_number_compare((x),(y)) > 0)
#endif
/* num_log.c */

View file

@ -646,6 +646,15 @@ static union {
#ifdef ECL_COMPLEX_FLOAT
#include <complex.h>
#ifndef CMPLXF
# define CMPLXF(x, y) ((float complex)((float)(x) + I * (float)(y)))
#endif
#ifndef CMPLX
# define CMPLX(x, y) ((double complex)((double)(x) + I * (double)(y)))
#endif
#ifndef CMPLXL
# define CMPLXL(x, y) ((long double complex)((long double)(x) + I * (long double)(y)))
#endif
#endif
#ifdef __cplusplus

View file

@ -126,7 +126,15 @@ Returns zero for non-complex numbers."
"Args: (number)
Returns a number that represents the sign of NUMBER. Returns NUMBER If it is
zero. Otherwise, returns the value of (/ NUMBER (ABS NUMBER))"
(if (zerop x) x (/ x (abs x))))
(if (complexp x)
(cis (atan (imagpart x) (realpart x)))
(let ((result (cond ((> x 0) 1)
((< x 0) -1)
(t ; x is 0 or NaN
x))))
(if (floatp x)
(float result x)
result))))
(defun cis (x)
"Args: (radians)
@ -143,7 +151,7 @@ RADIANS) and (SIN RADIANS) respectively."
(setf ,arg (float ,arg)))
(typecase ,arg
(single-float
(if ,restriction
(if (or ,restriction #+ieee-floating-point (ext:float-nan-p ,arg))
(ffi::c-inline (,arg) (:float) :float
,(format nil "~af(#0)" name)
:one-liner t)
@ -154,7 +162,7 @@ RADIANS) and (SIN RADIANS) respectively."
#-complex-float
(progn ,@gencomplex)))
(double-float
(if ,restriction
(if (or ,restriction #+ieee-floating-point (ext:float-nan-p ,arg))
(ffi::c-inline (,arg) (:double) :double
,(format nil "~a(#0)" name)
:one-liner t)
@ -165,7 +173,7 @@ RADIANS) and (SIN RADIANS) respectively."
#-complex-float
(progn ,@gencomplex)))
(long-float
(if ,restriction
(if (or ,restriction #+ieee-floating-point (ext:float-nan-p ,arg))
(ffi::c-inline (,arg) (:long-double) :long-double
,(format nil "~al(#0)" name)
:one-liner t)
@ -311,7 +319,7 @@ Returns the hyperbolic arc tangent of NUMBER."
(declare (number z) (si::c-local))
(/ (- (log (1+ z)) (log (- 1 z))) 2))
(defun ffloor (x &optional (y 1))
(defun ffloor (x &optional (y 1.0f0))
"Args: (number &optional (divisor 1))
Same as FLOOR, but returns a float as the first value."
(multiple-value-bind (i r) (floor x y)

View file

@ -121,7 +121,8 @@ as a second value."
;;; Approximate equality function
(defun approx= (x y &optional (eps (epsilon x)))
(<= (abs (/ (- x y) (max (abs x) 1))) eps))
(or (= x y)
(<= (abs (/ (- x y) (max (abs x) 1))) eps)))
(defun epsilon (number)
(etypecase number

View file

@ -18,6 +18,113 @@
,@body)
(si:trap-fpe bits t))))
(defconstant +number-subtypes-and-values+
'((fixnum 3)
(bignum #.(- most-negative-fixnum 13))
(ratio 3/5)
(single-float 1.23f0)
(double-float 4.56d0)
(long-float 7.89l0)
((complex fixnum) #C(2 3))
((complex integer) #.(complex 2 (+ most-positive-fixnum 3)))
((complex rational) #C(5/7 -4))
((complex single-float) #C(0.98f0 77.7f0))
((complex double-float) #C(0.87d0 88.8d0))
((complex long-float) #C(0.76l0 99.9l0))))
(defmacro for-all-number-subtypes ((variable &optional (type number) default-value) &body body)
(when (not (subtypep type 'number))
(error "~a must be a subtype of number" type))
`(progn ,@(loop :for (subtype value) :in +number-subtypes-and-values+
:if (subtypep subtype type)
:collect `(let ((,variable
,(if default-value
`(coerce ,default-value ',subtype)
value)))
,@body))))
(defmacro for-all-NaNs (variable &body body)
`(loop :for ,variable :in *ieee-fp.nan* :do
,@body))
(defmacro for-all-infinities (minus-infinity-var plus-infinity-var &body body)
`(loop :for ,minus-infinity-var :in (remove-if-not #'minusp *ieee-fp.inf*) :do
(loop :for ,plus-infinity-var :in (remove-if-not #'plusp *ieee-fp.inf*) :do
,@body)))
(defun type-contagion (&rest numbers)
(labels ((type-of-proper (number)
(etypecase number
(single-float 'single-float)
(double-float 'double-float)
(long-float 'long-float)
(ratio 'ratio)
(integer 'integer)
(complex `(complex ,(type-of-proper (realpart number))))))
(type-contagion1 (t1 &optional t2)
(cond ((not t2) t1)
((subtypep t1 'float)
(cond ((subtypep t2 'float)
(cond ((or (eql t1 'long-float)
(eql t2 'long-float))
'long-float)
((or (eql t1 'double-float)
(eql t2 'double-float))
'double-float)
((or (eql t1 'single-float)
(eql t2 'single-float))
'single-float)))
((subtypep t2 'rational)
t1)
((subtypep t2 'complex)
`(complex ,(type-contagion1 t1 (second t2))))))
((subtypep t1 'ratio)
(cond ((subtypep t2 'float) t2)
((subtypep t2 'rational) 'ratio)
((subtypep t2 'complex)
`(complex ,(type-contagion1 t1 (second t2))))))
((subtypep t1 'integer)
(cond ((or (subtypep t2 'float) (subtypep t2 'ratio)) t2)
((subtypep t2 'integer) 'integer)
((subtypep t2 'complex)
`(complex ,(type-contagion1 t1 (second t2))))))
((subtypep t1 'complex)
(cond ((or (subtypep t2 'float) (subtypep t2 'rational))
`(complex ,(type-contagion (second t1) t2)))
((subtypep t2 'complex)
`(complex ,(type-contagion1 (second t1) (second t2)))))))))
(reduce #'type-contagion1 (mapcar #'type-of-proper numbers))))
(defmacro type-and-value-check ((function &rest args) &rest values)
"Arguments ((function &rest args) &rest values)
Check that function applied to args returns numbers that are
mathematically equal to values and have the correct type according to
Common Lisp type contagion rules."
(if (eql function 'funcall)
(setf function (first args)
args (rest args))
(setf function `(quote ,function)))
`(is (every #'eql
(multiple-value-list (funcall ,function ,@args))
(mapcar #'(lambda (v)
(coerce v (type-contagion ,@args)))
(list ,@values)))))
;;; Constructing complex floats with infinities is difficult, because
;;; any finite real number x plus infinite imaginary number leads to
;;; an complex number with NaN real part due to a very literal
;;; application of calculation rules:
;;; x + I*∞ = x + (0 + I*1)*∞ = x + NaN + I*∞ = NaN + I*∞
;;; Therefore, when comparing the output of numeric functions
;;; returning complex infinities to expected results, we use the
;;; following function
(defmacro complex-equality (equality-comparison form realpart imagpart)
`(progn
(is (funcall ,equality-comparison (realpart ,form) ,realpart))
(is (funcall ,equality-comparison (imagpart ,form) ,imagpart))))
(test ieee-fp.0001.infinity-eql
(without-fpe-traps
(let ((sfni ext:single-float-negative-infinity)
@ -55,17 +162,280 @@
(is (ext:float-infinity-p ext:single-float-positive-infinity))
(is (ext:float-nan-p (si:nan))))
;;; Reported by: Robert Dodier
;;; URL: https://gitlab.com/embeddable-common-lisp/ecl/issues/299
(test ieee-fp.0003.b299
(finishes (< ext:double-float-negative-infinity 1/3))
(finishes (> ext:double-float-negative-infinity 1/3))
(finishes (< 1/3 ext:double-float-negative-infinity))
(finishes (> 1/3 ext:double-float-negative-infinity))
(finishes (< ext:double-float-negative-infinity (1+ most-positive-fixnum)))
(finishes (> ext:double-float-negative-infinity (1+ most-positive-fixnum)))
(finishes (< (1+ most-positive-fixnum) ext:double-float-negative-infinity))
(finishes (> (1+ most-positive-fixnum) ext:double-float-negative-infinity)))
(test ieee-fp.0004.nan-comparison
(for-all-number-subtypes (x real)
(for-all-NaNs NaN
(is (not (< NaN x)))
(is (not (< x NaN)))
(is (not (<= NaN x)))
(is (not (<= x NaN)))
(is (not (> NaN x)))
(is (not (> x NaN)))
(is (not (>= NaN x)))
(is (not (>= x NaN))))))
(test ieee-fp.0005.infinity-comparison
(for-all-number-subtypes (x real)
(for-all-infinities -infinity +infinity
(progn (is (< -infinity x))
(is (not (< x -infinity)))
(is (< x +infinity))
(is (not (< +infinity x))))
(progn (is (<= -infinity x))
(is (not (<= x -infinity)))
(is (<= x +infinity))
(is (not (<= +infinity x))))
(progn (is (not (> -infinity x)))
(is (> x -infinity))
(is (not (> x +infinity)))
(is (> +infinity x)))
(progn (is (not (>= -infinity x)))
(is (>= x -infinity))
(is (not (>= x +infinity)))
(is (>= +infinity x))))))
;;; We treat NaNs as missing data (like fmin/fmax C functions)
(test ieee-fp.0006.NaN-min/max
(for-all-NaNs NaN
(is (eql (min NaN) NaN))
(is (eql (max NaN) NaN))
(is (eql (min NaN NaN) NaN))
(is (eql (max NaN NaN) NaN))
(for-all-number-subtypes (x real)
(let* ((y (* 2 x))
(min (min x y))
(max (max x y)))
(is (= (min x NaN) x))
(is (= (max x NaN) x))
(is (= (min x y NaN) min))
(is (= (max x y NaN) max))))))
(test ieee-fp.0007.infinity-min/max
(for-all-number-subtypes (x real)
(let* ((y (* 2 x))
(max (max x y))
(min (min x y)))
(for-all-infinities -infinity +infinity
(is (= (min -infinity x y) -infinity))
(is (= (min +infinity x y) min))
(is (= (max -infinity x y) max))
(is (= (max +infinity x y) +infinity))))))
(test ieee-fp.0008.NaN-minusp/zerop/plusp
(without-fpe-traps
(for-all-NaNs NaN
(is (not (minusp NaN)))
(is (not (zerop NaN)))
(is (not (plusp NaN))))))
(test ieee-fp.0009.infinity-minusp/zerop/plusp
(for-all-infinities -infinity +infinity
(is (minusp -infinity))
(is (not (minusp +infinity)))
(is (not (zerop -infinity)))
(is (not (zerop +infinity)))
(is (not (plusp -infinity)))
(is (plusp +infinity))))
;;; If we supported NaN/infinity for rounding functions operating
;;; on floats, sensible behaviour would be as follows ...
#|
(test ieee-fp.0010.NaN-ffloor/fceiling/ftruncate/fround
(without-fpe-traps
(loop :for function :in '(ffloor fceiling ftruncate fround)
:do
(for-all-number-subtypes (x real)
(for-all-NaNs NaN
(type-and-value-check (funcall function NaN) NaN NaN)
(type-and-value-check (funcall function x NaN) NaN NaN)
(type-and-value-check (funcall function NaN x) NaN NaN)
(type-and-value-check (funcall function NaN NaN) NaN NaN)))
(for-all-number-subtypes (x float 0)
(type-and-value-check (funcall function x x) (ext:nan) (ext:nan))))))
(test ieee-fp.0011.infinity-ffloor/fceiling/ftruncate/fround
(without-fpe-traps
(loop :for function :in '(ffloor fceiling ftruncate fround)
:do
(for-all-infinities -infinity +infinity
;; +-infinity and 1
(type-and-value-check (funcall function -infinity)
-infinity (ext:nan))
(type-and-value-check (funcall function +infinity)
+infinity (ext:nan))
(for-all-number-subtypes (x real)
;; +-infinity and finite real
(type-and-value-check (funcall function -infinity (abs x))
-infinity (ext:nan))
(type-and-value-check (funcall function -infinity (- (abs x)))
+infinity (ext:nan))
(type-and-value-check (funcall function +infinity (abs x))
+infinity (ext:nan))
(type-and-value-check (funcall function +infinity (- (abs x)))
-infinity (ext:nan))
;; finite real and +-infinity
(type-and-value-check (funcall function (abs x) -infinity)
(case function
(ffloor -1)
(t 0))
(case function
(ffloor -infinity)
(t x)))
(type-and-value-check (funcall function (- (abs x)) -infinity)
(case function
(fceiling 1)
(t 0))
(case function
(fceiling +infinity)
(t x)))
(type-and-value-check (funcall function (abs x) +infinity)
(case function
(fceiling 1)
(t 0))
(case function
(fceiling -infinity)
(t x)))
(type-and-value-check (funcall function (- (abs x)) +infinity)
(case function
(ffloor -1)
(t 0))
(case function
(ffloor +infinity)
(t x))))
;; +-infinity and +-infinity
(type-and-value-check (funcall function -infinity -infinity)
(ext:nan) (ext:nan))
(type-and-value-check (funcall function -infinity +infinity)
(ext:nan) (ext:nan))
(type-and-value-check (funcall function +infinity -infinity)
(ext:nan) (ext:nan))
(type-and-value-check (funcall function +infinity +infinity)
(ext:nan) (ext:nan))
(for-all-number-subtypes (x float 0)
;; finite real and 0
(for-all-number-subtypes (y real)
(type-and-value-check (funcall function (abs y) x)
+infinity (ext:nan))
(type-and-value-check (funcall function (- (abs y)) x)
-infinity (ext:nan)))
;; +-infinity and 0
(type-and-value-check (funcall function -infinity 0)
-infinity (ext:nan))
(type-and-value-check (funcall function +infinity 0)
+infinity (ext:nan)))))))
|#
;;; ... but we don't, therefore everything throws arithmetic errors.
(test ieee-fp.0010.NaN-floor/ceiling/truncate/round/mod/rem
(loop :for function :in '(floor ceiling truncate round
ffloor fceiling ftruncate fround
mod rem)
:do
(for-all-NaNs NaN
(unless (member function '(mod rem))
(signals arithmetic-error (funcall function NaN)))
(signals arithmetic-error (funcall function NaN NaN))
(for-all-number-subtypes (x real)
(signals arithmetic-error (funcall function x NaN))
(signals arithmetic-error (funcall function NaN x))))
(for-all-number-subtypes (x float 0)
(signals arithmetic-error (funcall function x x)))))
(test ieee-fp.0011.infinity-floor/ceiling/truncate/round
(loop :for function :in '(floor ceiling truncate round
ffloor fceiling ftruncate fround
mod rem)
:do
(for-all-infinities -infinity +infinity
;; +-infinity and 1
(unless (member function '(mod rem))
(signals arithmetic-error (funcall function -infinity))
(signals arithmetic-error (funcall function +infinity)))
(for-all-number-subtypes (x real)
;; +-infinity and finite real
(signals arithmetic-error (funcall function -infinity x))
(signals arithmetic-error (funcall function +infinity x))
;; finite real and +-infinity
(signals arithmetic-error (funcall function x -infinity))
(signals arithmetic-error (funcall function x +infinity)))
;; +-infinity and +-infinity
(signals arithmetic-error (funcall function -infinity -infinity))
(signals arithmetic-error (funcall function -infinity +infinity))
(signals arithmetic-error (funcall function +infinity -infinity))
(signals arithmetic-error (funcall function +infinity +infinity))
(for-all-number-subtypes (x float 0)
;; finite real and 0
(for-all-number-subtypes (y real)
(signals arithmetic-error (funcall function y x)))
;; +-infinity and 0
(signals arithmetic-error (funcall function -infinity x))
(signals arithmetic-error (funcall function +infinity x))))))
;;; All numeric functions should return a NaN if given one
(test ieee-fp.0012.sticky-NaN
(without-fpe-traps
(for-all-NaNs NaN
(macrolet ((f-NaN (&rest functions)
`(progn ,@(loop :for f :in functions :collect
`(is (eql (,f NaN) NaN))))))
(f-NaN sin cos tan
asin acos atan
sinh cosh tanh
asinh acosh atanh
* + - / 1+ 1-
abs signum exp log sqrt
conjugate phase))
(for-all-number-subtypes (x real)
(macrolet ((f-NaN (&rest functions)
`(progn ,@(loop :for f :in functions :append
`((type-and-value-check (,f x NaN) NaN)
(type-and-value-check (,f NaN x) NaN)
(type-and-value-check (,f NaN NaN) NaN))))))
(f-NaN atan * + - /))
(macrolet ((f-NaN (&rest functions)
`(progn ,@(loop :for f :in functions :append
`((type-and-value-check (,f (abs x) NaN) NaN)
(type-and-value-check (,f NaN (abs x)) NaN)
(type-and-value-check (,f NaN NaN) NaN))))))
(f-NaN expt log)))
(is (eql (cis NaN)
(complex (coerce (ext:nan) (type-of NaN))
(coerce (ext:nan) (type-of NaN))))))))
(test ieee-fp.0013.infinity-sin/cos/tan
(without-fpe-traps
(for-all-infinities -infinity +infinity
(type-and-value-check (sin -infinity) (ext:nan))
(type-and-value-check (sin +infinity) (ext:nan))
#|(type-and-value-check (sin (complex 0 -infinity)) (complex 0 -infinity))
(type-and-value-check (sin (complex 0 +infinity)) (complex 0 +infinity))|#
(type-and-value-check (cos -infinity) (ext:nan))
(type-and-value-check (cos +infinity) (ext:nan))
#|(type-and-value-check (cos (complex 0 -infinity)) +infinity)
(type-and-value-check (cos (complex 0 +infinity)) +infinity)|#
(type-and-value-check (tan -infinity) (ext:nan))
(type-and-value-check (tan +infinity) (ext:nan))
#|(type-and-value-check (tan (complex 0 -infinity)) (complex 0 -1))
(type-and-value-check (tan (complex 0 +infinity)) (complex 0 +1))|#)))
#|(test ieee-fp.0014.infinity-asin/acos/atan
(for-all-infinities -infinity +infinity
(complex-equality #'approx= (asin -infinity) (- (/ pi 2)) +infinity)
(complex-equality #'approx= (asin +infinity) (/ pi 2) +infinity)
(complex-equality #'approx= (acos -infinity) pi -infinity)
(complex-equality #'approx= (acos +infinity) 0 -infinity)
(is (approx= (atan -infinity) (- (/ pi 2))))
(is (approx= (atan +infinity) (/ pi 2)))))|#
;;; Reported by: Raymond Toy
@ -79,16 +449,7 @@
(+zerop (elt) (and (zerop elt) (plusp (float-sign elt))))
(-zerop (elt) (and (zerop elt) (minusp (float-sign elt)))))
(test ieee-fp.0004.atan2-special-cases.nan-arg
(without-fpe-traps
;; (atan (anything) nan) -> nan
;; (atan nan (anything)) -> nan
(map () (lambda (n)
(is (si:float-nan-p (atan n (si:nan))))
(is (si:float-nan-p (atan (si:nan) n))))
(append *floats* *ieee-fp.inf* (list (si:nan))))))
(test ieee-fp.0005.atan2-special-case.zero-arg
(test ieee-fp.0015.atan2-special-case.zero-arg
(+zerop (atan +0.0 +0.0))
(-zerop (atan -0.0 +0.0))
(+pi-p (atan +0.0 -0.0))
@ -108,7 +469,7 @@
(is (-pi/2-p (atan (- n) -0.0))))
(remove-if-not #'plusp (append *floats* *ieee-fp.inf*)))))
(test ieee-fp.0006.atan2-special-case.inf-arg
(test ieee-fp.0016.atan2-special-case.inf-arg
;; (atan +-inf +inf) -> +-pi/4
(let ((+inf ext:single-float-positive-infinity)
(-inf ext:single-float-negative-infinity))
@ -125,3 +486,248 @@
(+pi/2-p (atan +inf n))
(-pi/2-p (atan -inf n))))
(remove-if-not #'plusp *floats*)))))
(test ieee-fp.0017.infinity-sinh/cosh/tanh
(for-all-infinities -infinity +infinity
(type-and-value-check (sinh -infinity) -infinity)
(type-and-value-check (sinh +infinity) +infinity)
(type-and-value-check (cosh -infinity) +infinity)
(type-and-value-check (cosh +infinity) +infinity)
(type-and-value-check (tanh -infinity) -1)
(type-and-value-check (tanh +infinity) +1)
#|(without-fpe-traps
(type-and-value-check (sinh (complex 0 -infinity)) (complex (ext:nan) (ext:nan)))
(type-and-value-check (sinh (complex 0 +infinity)) (complex (ext:nan) (ext:nan)))
(type-and-value-check (cosh (complex 0 -infinity)) (complex (ext:nan) (ext:nan)))
(type-and-value-check (cosh (complex 0 +infinity)) (complex (ext:nan) (ext:nan)))
(type-and-value-check (tanh (complex 0 -infinity)) (complex (ext:nan) (ext:nan)))
(type-and-value-check (tanh (complex 0 +infinity)) (complex (ext:nan) (ext:nan))))|#))
(test ieee-fp.0018.infinity-asinh/acosh/atanh
(for-all-infinities -infinity +infinity
(type-and-value-check (asinh -infinity) -infinity)
(type-and-value-check (asinh +infinity) +infinity)
(without-fpe-traps
(complex-equality #'approx= (acosh -infinity) +infinity pi))
(type-and-value-check (acosh +infinity) +infinity)
#| (complex-equality #'approx= (atanh -infinity) 0 (/ pi 2))
(complex-equality #'approx= (atanh +infinity) 0 (/ pi 2))|#))
(test ieee-fp.0019.infinity-+/-
(for-all-infinities -infinity +infinity
(type-and-value-check (+ -infinity) -infinity)
(type-and-value-check (+ +infinity) +infinity)
(type-and-value-check (+ -infinity -infinity) -infinity)
(without-fpe-traps
(type-and-value-check (+ -infinity +infinity) (ext:nan))
(type-and-value-check (+ +infinity -infinity) (ext:nan)))
(type-and-value-check (+ +infinity +infinity) +infinity)
(type-and-value-check (- -infinity) +infinity)
(type-and-value-check (- +infinity) -infinity)
(without-fpe-traps
(type-and-value-check (- -infinity -infinity) (ext:nan)))
(type-and-value-check (- +infinity -infinity) +infinity)
(type-and-value-check (- -infinity +infinity) -infinity)
(without-fpe-traps
(type-and-value-check (- +infinity +infinity) (ext:nan)))
(for-all-number-subtypes (x real)
(type-and-value-check (+ x -infinity) -infinity)
(type-and-value-check (+ x +infinity) +infinity)
(type-and-value-check (+ -infinity x) -infinity)
(type-and-value-check (+ +infinity x) +infinity)
(type-and-value-check (- x -infinity) +infinity)
(type-and-value-check (- x +infinity) -infinity)
(type-and-value-check (- -infinity x) -infinity)
(type-and-value-check (- +infinity x) +infinity))))
(test ieee-fp.0020.infinity-*//
(for-all-infinities -infinity +infinity
(type-and-value-check (* -infinity) -infinity)
(type-and-value-check (* +infinity) +infinity)
(type-and-value-check (* -infinity -infinity) +infinity)
(type-and-value-check (* -infinity +infinity) -infinity)
(type-and-value-check (* +infinity -infinity) -infinity)
(type-and-value-check (* +infinity +infinity) +infinity)
(is (eql (/ -infinity) (coerce -0.0 (type-of -infinity))))
(is (eql (/ +infinity) (coerce +0.0 (type-of +infinity))))
(without-fpe-traps
(type-and-value-check (/ -infinity -infinity) (ext:nan))
(type-and-value-check (/ -infinity +infinity) (ext:nan))
(type-and-value-check (/ +infinity -infinity) (ext:nan))
(type-and-value-check (/ +infinity +infinity) (ext:nan)))
(for-all-number-subtypes (x real)
(type-and-value-check (* (abs x) -infinity) -infinity)
(type-and-value-check (* (- (abs x)) -infinity) +infinity)
(type-and-value-check (* (abs x) +infinity) +infinity)
(type-and-value-check (* (- (abs x)) +infinity) -infinity)
(type-and-value-check (* -infinity (abs x)) -infinity)
(type-and-value-check (* -infinity (- (abs x))) +infinity)
(type-and-value-check (* +infinity (abs x)) +infinity)
(type-and-value-check (* +infinity (- (abs x))) -infinity)
(is (eql (/ (abs x) -infinity) (coerce -0.0 (type-contagion (abs x) -infinity))))
(is (eql (/ (- (abs x)) -infinity) (coerce +0.0 (type-contagion (- (abs x)) -infinity))))
(is (eql (/ (abs x) +infinity) (coerce +0.0 (type-contagion (abs x) +infinity))))
(is (eql (/ (- (abs x)) +infinity) (coerce -0.0 (type-contagion (- (abs x)) +infinity))))
(type-and-value-check (/ -infinity (abs x)) -infinity)
(type-and-value-check (/ -infinity (- (abs x))) +infinity)
(type-and-value-check (/ +infinity (abs x)) +infinity)
(type-and-value-check (/ +infinity (- (abs x))) -infinity))))
(test ieee-fp.0021.infinity-1+/1-
(for-all-infinities -infinity +infinity
(type-and-value-check (1- -infinity) -infinity)
(type-and-value-check (1- +infinity) +infinity)
(type-and-value-check (1+ -infinity) -infinity)
(type-and-value-check (1+ +infinity) +infinity)))
(test ieee-fp.0022.infinity-abs
(for-all-infinities -infinity +infinity
(type-and-value-check (abs -infinity) +infinity)
(type-and-value-check (abs +infinity) +infinity)))
(test ieee-fp.0023.infinity-exp
(for-all-infinities -infinity +infinity
(type-and-value-check (exp -infinity) 0)
(type-and-value-check (exp +infinity) +infinity)))
(test ieee-fp.0024.infinity-expt
(without-fpe-traps
(flet ((relaxed-eql (x y)
;; (relaxed-eql +0.0 -0.0) => true
(and (eql (type-of x) (type-of y))
(= x y))))
(for-all-infinities -infinity +infinity
;; x^∞: |x| < 1
(for-all-number-subtypes (x (and real (not integer)) 1/2)
(type-and-value-check (expt x -infinity) +infinity)
#| (complex-equality #'eql (expt (- x) -infinity)
(coerce +infinity (type-contagion x -infinity)) ;
(coerce (ext:nan) (type-contagion x -infinity))) |#
(type-and-value-check (expt x +infinity) 0)
#| (complex-equality #'relaxed-eql (expt (- x) +infinity)
(coerce 0.0 (type-contagion x +infinity)) ;
(coerce 0.0 (type-contagion x +infinity))) |#)
;; x^∞: |x| = 1
(for-all-number-subtypes (x (or fixnum float) 1)
(type-and-value-check (expt x -infinity) 1)
(without-fpe-traps
(complex-equality #'eql (expt (- x) -infinity)
(coerce (ext:nan) (type-contagion x -infinity))
(coerce (ext:nan) (type-contagion x -infinity))))
(type-and-value-check (expt x +infinity) 1)
(without-fpe-traps
(complex-equality #'eql (expt (- x) +infinity)
(coerce (ext:nan) (type-contagion x +infinity))
(coerce (ext:nan) (type-contagion x +infinity)))))
;; x^∞: |x| > 1
(for-all-number-subtypes (x (or fixnum float) 2)
(type-and-value-check (expt x -infinity) 0)
#| (complex-equality #'relaxed-eql (expt (- x) -infinity)
(coerce 0.0 (type-contagion x -infinity)) ;
(coerce 0.0 (type-contagion x -infinity))) |#
(type-and-value-check (expt x +infinity) +infinity)
#| (complex-equality #'eql (expt (- x) +infinity)
(coerce +infinity (type-contagion x +infinity)) ;
(coerce (ext:nan) (type-contagion x +infinity))) |#)
;; ∞^y: even y
(for-all-number-subtypes (y float 2)
(is (eql (realpart (expt -infinity y))
(coerce +infinity (type-contagion -infinity y))))
(complex-equality #'relaxed-eql (expt -infinity (- y))
(coerce 0.0 (type-contagion -infinity y))
(coerce 0.0 (type-contagion -infinity y)))
(type-and-value-check (expt +infinity y) +infinity)
(type-and-value-check (expt +infinity (- y)) 0.0))
(for-all-number-subtypes (y fixnum 2)
(type-and-value-check (expt -infinity y) +infinity)
(type-and-value-check (expt -infinity (- y)) 0.0)
(type-and-value-check (expt +infinity y) +infinity)
(type-and-value-check (expt +infinity (- y)) 0.0))
;; ∞^y: odd y
(for-all-number-subtypes (y float 3)
;; Should be -infinity, but C99 complex functions return +infinity
#| (is (eql (realpart (expt -infinity y))
(coerce -infinity (type-contagion -infinity y)))) |#
(complex-equality #'relaxed-eql (expt -infinity (- y))
(coerce 0.0 (type-contagion -infinity y))
(coerce 0.0 (type-contagion -infinity y)))
(type-and-value-check (expt +infinity y) +infinity)
(type-and-value-check (expt +infinity (- y)) 0.0))
(for-all-number-subtypes (y fixnum 3)
(type-and-value-check (expt -infinity y) -infinity)
(complex-equality #'relaxed-eql (expt -infinity (- y))
(coerce 0.0 (type-contagion -infinity y))
(coerce 0.0 (type-contagion -infinity y)))
(type-and-value-check (expt +infinity y) +infinity)
(type-and-value-check (expt +infinity (- y)) 0.0))
;; ∞^y: non-integer y
(for-all-number-subtypes (y (and real (not integer)) 5/2)
(complex-equality #'relaxed-eql (expt -infinity (- y))
(coerce 0.0 (type-contagion -infinity y))
(coerce 0.0 (type-contagion -infinity y)))
(type-and-value-check (expt +infinity y) +infinity)
(type-and-value-check (expt +infinity (- y)) 0.0))
;; ∞^∞
#| (complex-equality #'eql (expt -infinity -infinity)
(coerce 0.0 (type-contagion -infinity -infinity)) ;
(coerce -0.0 (type-contagion -infinity -infinity))) |#
#| (complex-equality #'eql (expt -infinity +infinity)
(coerce +infinity (type-contagion -infinity +infinity)) ;
(coerce (ext:nan) (type-contagion -infinity +infinity))) |#
(type-and-value-check (expt +infinity -infinity) 0.0)
(type-and-value-check (expt +infinity +infinity) +infinity)))))
(test ieee-fp.0025.infinity-log
(for-all-infinities -infinity +infinity
(complex-equality #'approx= (log -infinity) +infinity pi)
(type-and-value-check (log +infinity) +infinity)))
(test ieee-fp.0026.infinity-signum
(for-all-infinities -infinity +infinity
(type-and-value-check (signum -infinity) -1)
(type-and-value-check (signum +infinity) +1)))
(test ieee-fp.0027.infinity-sqrt
(for-all-infinities -infinity +infinity
(without-fpe-traps
(is (= (imagpart (sqrt -infinity)) +infinity)))
(type-and-value-check (sqrt +infinity) +infinity)))
(test ieee-fp.0028.infinity-cis
(without-fpe-traps
(for-all-infinities -infinity +infinity
(is (eql (cis -infinity)
(complex (coerce (ext:nan) (type-of -infinity))
(coerce (ext:nan) (type-of -infinity)))))
(is (eql (cis +infinity)
(complex (coerce (ext:nan) (type-of +infinity))
(coerce (ext:nan) (type-of +infinity))))))))
(test ieee-fp.0029.infinity-conjugate
(for-all-infinities -infinity +infinity
(type-and-value-check (conjugate -infinity) -infinity)
(type-and-value-check (conjugate +infinity) +infinity)
#| (for-all-number-subtypes (x real)
(type-and-value-check (conjugate (complex x -infinity)) (complex x +infinity))
(type-and-value-check (conjugate (complex x +infinity)) (complex x -infinity))) |#))
(test ieee-fp.0030.infinity-phase
(for-all-infinities -infinity +infinity
(is (approx= (phase -infinity) (- pi)))
(is (approx= (phase +infinity) 0))
#| (is (approx= (phase (complex 0 -infinity)) (- (/ pi 2))))
(is (approx= (phase (complex 0 +infinity)) (+ (/ pi 2))))
(is (approx= (phase (complex -infinity -infinity)) (- (* pi 3/4))))
(is (approx= (phase (complex -infinity +infinity)) (+ (* pi 3/4))))
(is (approx= (phase (complex +infinity -infinity)) (- (* pi 1/4))))
(is (approx= (phase (complex +infinity +infinity)) (+ (* pi 1/4)))) |#))

View file

@ -129,6 +129,12 @@
ext:double-float-positive-infinity
ext:long-float-positive-infinity))
#+ieee-floating-point
(defparameter *ieee-fp.nan*
(list (coerce (ext:nan) 'single-float)
(coerce (ext:nan) 'double-float)
(coerce (ext:nan) 'long-float)))
(defparameter *ratios*
'(1/3 1/1000 1/1000000000000000 -10/3 -1000/7 -987129387912381/13612986912361
189729874978126783786123/1234678123487612347896123467851234671234))