From f83ba9a73da09573f666fde5e7c4d883b48da7e0 Mon Sep 17 00:00:00 2001 From: jjgarcia Date: Sat, 2 Aug 2008 18:31:10 +0000 Subject: [PATCH] Increased precision and fixed problems with long-float/long double and special functions --- src/c/num_sfun.d | 62 ++++++++++++++++------------------------------ src/lsp/numlib.lsp | 43 +++++++++++++++++++++----------- 2 files changed, 49 insertions(+), 56 deletions(-) diff --git a/src/c/num_sfun.d b/src/c/num_sfun.d index 1cfc93c04..1c9cdaaab 100644 --- a/src/c/num_sfun.d +++ b/src/c/num_sfun.d @@ -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); diff --git a/src/lsp/numlib.lsp b/src/lsp/numlib.lsp index bb6d8eef8..d29bf8947 100644 --- a/src/lsp/numlib.lsp +++ b/src/lsp/numlib.lsp @@ -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))