log10 etc in C

[Kevin Ryde]

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 ()
 {