|
From: | Pat Lasswell |
Subject: | Re: log10 etc in C |
Date: | Thu, 7 Sep 2006 16:12:28 -0700 |
I'm looking at moving some of log, log10, etc from boot-9.scm into
numbers.c. For log10 it means using the log10() func directly, and
for the others gives the chance to use the libc clog() etc.
--- numbers.c.~1.281.2.4.~ 2006-07-13 11:35:24.000000000 +1000
+++ numbers.c 2006-09-08 09:02:11.000000000 +1000
@@ -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,11 @@
#include "libguile/discouraged.h"
+/* value per glibc, if not already defined */
+#ifndef M_LOG10E
+#define M_LOG10E 0.43429448190325182765
+#endif
+
/*
@@ -150,6 +159,21 @@
#endif
}
+#if HAVE_COMPLEX_DOUBLE
+
+/* 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))
+
+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;
@@ -5977,6 +6001,122 @@
return scm_is_true (scm_number_p (z));
}
+
+/* In the following it might save a bit of code to use scm_to_complex_double
+ and call the "clog" or whatever function for both real and complex
+ values. But that's not done in case the complex funcs aren't optimized
+ to look for real-only arguments. We have to test SCM_COMPLEXP anyway, so
+ may as well use that test to dispatch direct to the real or complex libc
+ function. */
+
+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. */
+ return scm_from_double (log (scm_to_double (z)));
+ }
+}
+#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))
+ {
+#if HAVE_COMPLEX_DOUBLE
+ 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. */
+ return scm_from_double (log10 (scm_to_double (z)));
+ }
+}
+#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 ( address@hidden).")
+#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
+ double re = SCM_COMPLEX_REAL (z);
+ double im = SCM_COMPLEX_IMAG (z);
+ double ee = exp (re);
+ return scm_c_make_rectangular (ee * cos(im), ee * sin(im));
+#endif
+ }
+ else
+ {
+ /* When z is a negative bignum, conversion to double overflows, giving
+ -inf, 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),
+ "")
+#define FUNC_NAME s_scm_sqrt
+{
+ if (SCM_COMPLEXP (x))
+ {
+#if HAVE_COMPLEX_DOUBLE
+ 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 ()
{
_______________________________________________
Guile-devel mailing list
address@hidden
http://lists.gnu.org/mailman/listinfo/guile-devel
[Prev in Thread] | Current Thread | [Next in Thread] |