diff options
author | Ludovic Courtès <ludo@gnu.org> | 2008-09-10 22:27:30 +0200 |
---|---|---|
committer | Ludovic Courtès <ludo@gnu.org> | 2008-09-10 22:27:30 +0200 |
commit | 35a9197ccc91a3663313e1bf7d369101754a1075 (patch) | |
tree | 0389f04f4267d26b70c8e8cd54329be28ccf438b /libguile/numbers.c | |
parent | e7bca22779c68b800c75fdad8841440dfeb32f8d (diff) | |
parent | f30e1bdf97ae8b2b2918da585f887a4d3a23a347 (diff) | |
download | guile-35a9197ccc91a3663313e1bf7d369101754a1075.tar.gz |
Merge commit 'f30e1bdf97ae8b2b2918da585f887a4d3a23a347' into boehm-demers-weiser-gc
Conflicts:
libguile/Makefile.am
libguile/coop-pthreads.c
libguile/gc-freelist.c
libguile/gc-segment.c
libguile/gc.c
libguile/private-gc.h
test-suite/tests/environments.nottest
Diffstat (limited to 'libguile/numbers.c')
-rw-r--r-- | libguile/numbers.c | 165 |
1 files changed, 164 insertions, 1 deletions
diff --git a/libguile/numbers.c b/libguile/numbers.c index 1a662b015..3be5314e0 100644 --- a/libguile/numbers.c +++ b/libguile/numbers.c @@ -40,7 +40,7 @@ */ -/* tell glibc (2.3) to give prototype for C99 trunc() */ +/* tell glibc (2.3) to give prototype for C99 trunc(), csqrt(), etc */ #define _GNU_SOURCE #if HAVE_CONFIG_H @@ -51,6 +51,10 @@ #include <ctype.h> #include <string.h> +#if HAVE_COMPLEX_H +#include <complex.h> +#endif + #include "libguile/_scm.h" #include "libguile/feature.h" #include "libguile/ports.h" @@ -66,6 +70,14 @@ #include "libguile/discouraged.h" +/* values per glibc, if not already defined */ +#ifndef M_LOG10E +#define M_LOG10E 0.43429448190325182765 +#endif +#ifndef M_PI +#define M_PI 3.14159265358979323846 +#endif + /* @@ -150,6 +162,21 @@ xisnan (double x) #endif } + +/* For an SCM object Z which is a complex number (ie. satisfies + SCM_COMPLEXP), return its value as a C level "complex double". */ +#define SCM_COMPLEX_VALUE(z) \ + (SCM_COMPLEX_REAL (z) + _Complex_I * SCM_COMPLEX_IMAG (z)) + +/* Convert a C "complex double" to an SCM value. */ +#if HAVE_COMPLEX_DOUBLE +static SCM +scm_from_complex_double (complex double z) +{ + return scm_c_make_rectangular (creal (z), cimag (z)); +} +#endif /* HAVE_COMPLEX_DOUBLE */ + static mpz_t z_negative_one; @@ -5978,6 +6005,142 @@ scm_is_number (SCM z) return scm_is_true (scm_number_p (z)); } + +/* In the following functions we dispatch to the real-arg funcs like log() + when we know the arg is real, instead of just handing everything to + clog() for instance. This is in case clog() doesn't optimize for a + real-only case, and because we have to test SCM_COMPLEXP anyway so may as + well use it to go straight to the applicable C func. */ + +SCM_DEFINE (scm_log, "log", 1, 0, 0, + (SCM z), + "Return the natural logarithm of @var{z}.") +#define FUNC_NAME s_scm_log +{ + if (SCM_COMPLEXP (z)) + { +#if HAVE_COMPLEX_DOUBLE + return scm_from_complex_double (clog (SCM_COMPLEX_VALUE (z))); +#else + double re = SCM_COMPLEX_REAL (z); + double im = SCM_COMPLEX_IMAG (z); + return scm_c_make_rectangular (log (hypot (re, im)), + atan2 (im, re)); +#endif + } + else + { + /* ENHANCE-ME: When z is a bignum the logarithm will fit a double + although the value itself overflows. */ + double re = scm_to_double (z); + double l = log (fabs (re)); + if (re >= 0.0) + return scm_from_double (l); + else + return scm_c_make_rectangular (l, M_PI); + } +} +#undef FUNC_NAME + + +SCM_DEFINE (scm_log10, "log10", 1, 0, 0, + (SCM z), + "Return the base 10 logarithm of @var{z}.") +#define FUNC_NAME s_scm_log10 +{ + if (SCM_COMPLEXP (z)) + { + /* Mingw has clog() but not clog10(). (Maybe it'd be worth using + clog() and a multiply by M_LOG10E, rather than the fallback + log10+hypot+atan2.) */ +#if HAVE_COMPLEX_DOUBLE && HAVE_CLOG10 + return scm_from_complex_double (clog10 (SCM_COMPLEX_VALUE (z))); +#else + double re = SCM_COMPLEX_REAL (z); + double im = SCM_COMPLEX_IMAG (z); + return scm_c_make_rectangular (log10 (hypot (re, im)), + M_LOG10E * atan2 (im, re)); +#endif + } + else + { + /* ENHANCE-ME: When z is a bignum the logarithm will fit a double + although the value itself overflows. */ + double re = scm_to_double (z); + double l = log10 (fabs (re)); + if (re >= 0.0) + return scm_from_double (l); + else + return scm_c_make_rectangular (l, M_LOG10E * M_PI); + } +} +#undef FUNC_NAME + + +SCM_DEFINE (scm_exp, "exp", 1, 0, 0, + (SCM z), + "Return @math{e} to the power of @var{z}, where @math{e} is the\n" + "base of natural logarithms (2.71828@dots{}).") +#define FUNC_NAME s_scm_exp +{ + if (SCM_COMPLEXP (z)) + { +#if HAVE_COMPLEX_DOUBLE + return scm_from_complex_double (cexp (SCM_COMPLEX_VALUE (z))); +#else + return scm_c_make_polar (exp (SCM_COMPLEX_REAL (z)), + SCM_COMPLEX_IMAG (z)); +#endif + } + else + { + /* When z is a negative bignum the conversion to double overflows, + giving -infinity, but that's ok, the exp is still 0.0. */ + return scm_from_double (exp (scm_to_double (z))); + } +} +#undef FUNC_NAME + + +SCM_DEFINE (scm_sqrt, "sqrt", 1, 0, 0, + (SCM x), + "Return the square root of @var{z}. Of the two possible roots\n" + "(positive and negative), the one with the a positive real part\n" + "is returned, or if that's zero then a positive imaginary part.\n" + "Thus,\n" + "\n" + "@example\n" + "(sqrt 9.0) @result{} 3.0\n" + "(sqrt -9.0) @result{} 0.0+3.0i\n" + "(sqrt 1.0+1.0i) @result{} 1.09868411346781+0.455089860562227i\n" + "(sqrt -1.0-1.0i) @result{} 0.455089860562227-1.09868411346781i\n" + "@end example") +#define FUNC_NAME s_scm_sqrt +{ + if (SCM_COMPLEXP (x)) + { +#if HAVE_COMPLEX_DOUBLE && HAVE_USABLE_CSQRT + return scm_from_complex_double (csqrt (SCM_COMPLEX_VALUE (x))); +#else + double re = SCM_COMPLEX_REAL (x); + double im = SCM_COMPLEX_IMAG (x); + return scm_c_make_polar (sqrt (hypot (re, im)), + 0.5 * atan2 (im, re)); +#endif + } + else + { + double xx = scm_to_double (x); + if (xx < 0) + return scm_c_make_rectangular (0.0, sqrt (-xx)); + else + return scm_from_double (sqrt (xx)); + } +} +#undef FUNC_NAME + + + void scm_init_numbers () { |