Increased precision and fixed problems with long-float/long double and special functions

This commit is contained in:
jjgarcia 2008-08-02 18:31:10 +00:00
parent 2ef8b05d4b
commit f83ba9a73d
2 changed files with 49 additions and 56 deletions

View file

@ -542,16 +542,11 @@ cl_sin(cl_object x)
z = x + I y
sin(z) = sinh(I z) = sinh(-y + I x)
*/
double dx = ecl_to_double(x->complex.real);
double dy = ecl_to_double(x->complex.imag);
double a = sin(dx) * cosh(dy);
double b = cos(dx) * sinh(dy);
if (type_of(x->complex.real) != t_doublefloat)
output = ecl_make_complex(ecl_make_singlefloat(a),
ecl_make_singlefloat(b));
else
output = ecl_make_complex(ecl_make_doublefloat(a),
ecl_make_doublefloat(b));
cl_object dx = x->complex.real;
cl_object dy = x->complex.imag;
cl_object a = ecl_times(cl_sin(dx), cl_cosh(dy));
cl_object b = ecl_times(cl_cos(dx), cl_sinh(dy));
output = ecl_make_complex(a, b);
break;
}
default:
@ -588,16 +583,11 @@ cl_cos(cl_object x)
z = x + I y
cos(z) = cosh(I z) = cosh(-y + I x)
*/
double dx = ecl_to_double(x->complex.real);
double dy = ecl_to_double(x->complex.imag);
double a = cos(dx) * cosh(dy);
double b = -sin(dx) * sinh(dy);
if (type_of(x->complex.real) != t_doublefloat)
output = ecl_make_complex(ecl_make_singlefloat(a),
ecl_make_singlefloat(b));
else
output = ecl_make_complex(ecl_make_doublefloat(a),
ecl_make_doublefloat(b));
cl_object dx = x->complex.real;
cl_object dy = x->complex.imag;
cl_object a = ecl_times(cl_cos(dx), cl_cosh(dy));
cl_object b = ecl_times(ecl_negate(cl_sin(dx)), cl_sinh(dy));
output = ecl_make_complex(a, b);
break;
}
default:
@ -671,16 +661,11 @@ cl_sinh(cl_object x)
= (exp(x)*(cos(y)+Isin(y))-exp(-x)*(cos(y)-Isin(y)))/2
= sinh(x)*cos(y) + Icosh(x)*sin(y);
*/
double dx = ecl_to_double(x->complex.real);
double dy = ecl_to_double(x->complex.imag);
double a = sinh(dx) * cos(dy);
double b = cosh(dx) * sin(dy);
if (type_of(x->complex.real) != t_doublefloat)
output = ecl_make_complex(ecl_make_singlefloat(a),
ecl_make_singlefloat(b));
else
output = ecl_make_complex(ecl_make_doublefloat(a),
ecl_make_doublefloat(b));
cl_object dx = x->complex.real;
cl_object dy = x->complex.imag;
cl_object a = ecl_times(cl_sinh(dx), cl_cos(dy));
cl_object b = ecl_times(cl_cosh(dx), cl_sin(dy));
output = ecl_make_complex(a, b);
break;
}
default:
@ -719,16 +704,11 @@ cl_cosh(cl_object x)
= (exp(x)*(cos(y)+Isin(y))+exp(-x)*(cos(y)-Isin(y)))/2
= cosh(x)*cos(y) + Isinh(x)*sin(y);
*/
double dx = ecl_to_double(x->complex.real);
double dy = ecl_to_double(x->complex.imag);
double a = cosh(dx) * cos(dy);
double b = sinh(dx) * sin(dy);
if (type_of(x->complex.real) != t_doublefloat)
output = ecl_make_complex(ecl_make_singlefloat(a),
ecl_make_singlefloat(b));
else
output = ecl_make_complex(ecl_make_doublefloat(a),
ecl_make_doublefloat(b));
cl_object dx = x->complex.real;
cl_object dy = x->complex.imag;
cl_object a = ecl_times(cl_cosh(dx), cl_cos(dy));
cl_object b = ecl_times(cl_sinh(dx), cl_sin(dy));
output = ecl_make_complex(a, b);
break;
}
default:
@ -758,7 +738,7 @@ cl_tanh(cl_object x)
output = ecl_make_doublefloat(tanh(df(x))); break;
#ifdef ECL_LONG_FLOAT
case t_longfloat:
output = make_longfloat(coshl(ecl_long_float(x))); break;
output = make_longfloat(tanhl(ecl_long_float(x))); break;
#endif
case t_complex: {
cl_object a = cl_sinh(x);

View file

@ -89,6 +89,18 @@ Returns a complex number whose realpart and imagpart are the values of (COS
RADIANS) and (SIN RADIANS) respectively."
(exp (* imag-one x)))
#-ecl-min
(eval-when (:compile-toplevel)
(defmacro c-num-op (name arg)
#+long-float
`(ffi::c-inline (,arg) (:long-double) :long-double
,(format nil "~al(#0)" name)
:one-liner t)
#-long-float
`(ffi::c-inline (,arg) (:double) :double
,(format nil "~a(#0)" name)
:one-liner t)))
(defun asin (x)
"Args: (number)
Returns the arc sine of NUMBER."
@ -96,11 +108,10 @@ Returns the arc sine of NUMBER."
(complex-asin x)
#-ecl-min
(let* ((x (float x))
(xr (float x 1d0)))
(declare (double-float xr))
(xr (float x 1l0)))
(declare (long-float xr))
(if (and (<= -1.0 xr) (<= xr 1.0))
(float (ffi::c-inline (xr) (:double) :double "asin(#0)" :one-liner t)
x)
(float (c-num-op "asin" xr) x)
(complex-asin x)))))
;; Ported from CMUCL
@ -120,11 +131,10 @@ Returns the arc cosine of NUMBER."
(complex-acos x)
#-ecl-min
(let* ((x (float x))
(xr (float x 1d0)))
(declare (double-float xr))
(xr (float x 1l0)))
(declare (long-float xr))
(if (and (<= -1.0 xr) (<= xr 1.0))
(float (ffi::c-inline (xr) (:double) :double "acos(#0)" :one-liner t)
(float x))
(float (c-num-op "acos" xr) (float x))
(complex-acos x)))))
;; Ported from CMUCL
@ -143,6 +153,12 @@ Returns the arc cosine of NUMBER."
(ffi:clines "double acosh(double x) { return log(x+sqrt((x-1)*(x+1))); }")
(ffi:clines "double atanh(double x) { return log((1+x)/(1-x))/2; }"))
#+(and long-float (not ecl-min) win32)
(progn
(ffi:clines "double asinhl(long double x) { return logl(x+sqrtl(1.0+x*x)); }")
(ffi:clines "double acoshl(long double x) { return logl(x+sqrtl((x-1)*(x+1))); }")
(ffi:clines "double atanhl(long double x) { return logl((1+x)/(1-x))/2; }"))
;; Ported from CMUCL
(defun asinh (x)
"Args: (number)
@ -154,8 +170,7 @@ Returns the hyperbolic arc sine of NUMBER."
(complex (imagpart result)
(- (realpart result))))
#-(or ecl-min)
(float (ffi:c-inline (x) (:double) :double "asinh(#0)" :one-liner t)
(float x))))
(float (c-num-op "asinh" x) (float x))))
;; Ported from CMUCL
(defun acosh (x)
@ -169,8 +184,7 @@ Returns the hyperbolic arc cosine of NUMBER."
(xr (float x 1d0)))
(declare (double-float xr))
(if (<= 1.0 xr)
(float (ffi::c-inline (xr) (:double) :double "acosh(#0)" :one-liner t)
(float x))
(float (c-num-op "acosh" xr) (float x))
(complex-acosh x)))))
(defun complex-acosh (z)
@ -192,12 +206,11 @@ Returns the hyperbolic arc tangent of NUMBER."
(xr (float x 1d0)))
(declare (double-float xr))
(if (and (<= -1.0 xr) (<= xr 1.0))
(float (ffi::c-inline (xr) (:double) :double "atanh(#0)" :one-liner t)
(float x))
(float (c-num-op "atanh" xr) (float x))
(complex-atanh x)))))
(defun complex-atanh (z)
(declare (number x) (si::c-local))
(declare (number z) (si::c-local))
(/ (- (log (1+ z)) (log (- 1 z))) 2))
(defun ffloor (x &optional (y 1))