flushing out trig functions for complex numbers

eval.c

@@ -1149,7 +1149,14 @@ sexp sexp_register_optimization (sexp ctx sexp_api_params(self, n), sexp f, sexp
 #define maybe_convert_ratio(z)
 #endif
 
-#define define_math_op(name, cname)                                     \
+#if SEXP_USE_COMPLEX
+#define maybe_convert_complex(z, f)             \
+  else if (sexp_complexp(z)) return sexp_complex_normalize(f(ctx, z));
+#else
+#define maybe_convert_complex(z, f)
+#endif
+
+#define define_math_op(name, cname, complexf)                           \
   static sexp name (sexp ctx sexp_api_params(self, n), sexp z) {        \
     double d;                                                           \
     if (sexp_flonump(z))                                                \
@@ -1158,23 +1165,24 @@ sexp sexp_register_optimization (sexp ctx sexp_api_params(self, n), sexp f, sexp
       d = (double)sexp_unbox_fixnum(z);                                 \
     maybe_convert_ratio(z)                                              \
     maybe_convert_bignum(z)                                             \
+    maybe_convert_complex(z, complexf)                                  \
     else                                                                \
       return sexp_type_exception(ctx, self, SEXP_NUMBER, z);            \
     return sexp_make_flonum(ctx, cname(d));                             \
   }
 
-define_math_op(sexp_exp, exp)
+define_math_op(sexp_exp, exp, sexp_complex_exp)
-define_math_op(sexp_log, log)
+define_math_op(sexp_log, log, sexp_complex_log)
-define_math_op(sexp_sin, sin)
+define_math_op(sexp_sin, sin, sexp_complex_sin)
-define_math_op(sexp_cos, cos)
+define_math_op(sexp_cos, cos, sexp_complex_cos)
-define_math_op(sexp_tan, tan)
+define_math_op(sexp_tan, tan, sexp_complex_math_error)
-define_math_op(sexp_asin, asin)
+define_math_op(sexp_asin, asin, sexp_complex_asin)
-define_math_op(sexp_acos, acos)
+define_math_op(sexp_acos, acos, sexp_complex_acos)
-define_math_op(sexp_atan, atan)
+define_math_op(sexp_atan, atan, sexp_complex_math_error)
-define_math_op(sexp_round, round)
+define_math_op(sexp_round, round, sexp_complex_math_error)
-define_math_op(sexp_trunc, trunc)
+define_math_op(sexp_trunc, trunc, sexp_complex_math_error)
-define_math_op(sexp_floor, floor)
+define_math_op(sexp_floor, floor, sexp_complex_math_error)
-define_math_op(sexp_ceiling, ceil)
+define_math_op(sexp_ceiling, ceil, sexp_complex_math_error)
 
 static sexp sexp_sqrt (sexp ctx sexp_api_params(self, n), sexp z) {
   int negativep = 0;
@@ -1186,6 +1194,7 @@ static sexp sexp_sqrt (sexp ctx sexp_api_params(self, n), sexp z) {
     d = (double)sexp_unbox_fixnum(z);
   maybe_convert_bignum(z)       /* XXXX add bignum sqrt */
   maybe_convert_ratio(z)        /* XXXX add ratio sqrt */
+  maybe_convert_complex(z, sexp_complex_sqrt)
   else
     return sexp_type_exception(ctx, self, SEXP_NUMBER, z);
   sexp_gc_preserve1(ctx, res);

include/chibi/bignum.h

@@ -46,6 +46,14 @@ sexp sexp_ratio_normalize (sexp ctx, sexp rat, sexp in);
 #if SEXP_USE_COMPLEX
 sexp sexp_make_complex (sexp ctx, sexp real, sexp image);
 sexp sexp_complex_normalize (sexp real);
+sexp sexp_complex_math_error (sexp ctx, sexp z);
+sexp sexp_complex_sqrt (sexp ctx, sexp z);
+sexp sexp_complex_exp (sexp ctx, sexp z);
+sexp sexp_complex_log (sexp ctx, sexp z);
+sexp sexp_complex_sin (sexp ctx, sexp z);
+sexp sexp_complex_cos (sexp ctx, sexp z);
+sexp sexp_complex_asin (sexp ctx, sexp z);
+sexp sexp_complex_acos (sexp ctx, sexp z);
 #endif
 
 #endif  /* ! SEXP_BIGNUM_H */

include/chibi/sexp.h

@@ -662,12 +662,18 @@ SEXP_API sexp sexp_make_unsigned_integer(sexp ctx, sexp_luint_t x);
 #define sexp_negativep(x) (sexp_exact_negativep(x) ||                   \
                            (sexp_flonump(x) && sexp_flonum_value(x) < 0))
 
-#define sexp_negate(x)                                  \
+#define sexp_negate_exact(x)                            \
   if (sexp_bignump(x))                                  \
     sexp_bignum_sign(x) = -sexp_bignum_sign(x);         \
   else if (sexp_fixnump(x))                             \
     x = sexp_fx_neg(x);
 
+#define sexp_negate(x)                                  \
+  if (sexp_flonump(x))                                  \
+    sexp_flonum_value(x) = -sexp_flonum_value(x);       \
+  else                                                  \
+    sexp_negate_exact(x)
+
 #if SEXP_USE_FLONUMS || SEXP_USE_BIGNUMS
 #define sexp_uint_value(x) ((sexp_uint_t)(sexp_fixnump(x) ? sexp_unbox_fixnum(x) : sexp_bignum_data(x)[0]))
 #define sexp_sint_value(x) ((sexp_sint_t)(sexp_fixnump(x) ? sexp_unbox_fixnum(x) : sexp_bignum_sign(x)*sexp_bignum_data(x)[0]))

lib/init.scm

@@ -894,7 +894,7 @@
 (cond-expand
  (complex
   (define (exact-complex? x)
-    (and (complex? x) (exact? (complex-real x)) (exact? (complex-imag x))))))
+    (and (%complex? x) (exact? (complex-real x)) (exact? (complex-imag x))))))
 
 (cond-expand
  (ratios
@@ -931,7 +931,7 @@
 (cond-expand
  (complex
   (define (inexact? x)
-    (if (flonum? x) #t (and (complex? x) (not (exact-complex? x))))))
+    (if (flonum? x) #t (and (%complex? x) (not (exact-complex? x))))))
  (else (define inexact? flonum?)))
 (define (exact-integer? x) (if (fixnum? x) #t (bignum? x)))
 (define (integer? x)
@@ -939,7 +939,7 @@
 (define (number? x) (if (inexact? x) #t (exact? x)))
 (define complex? number?)
 (cond-expand
- (complex (define (rational? x) (and (number? x) (not (complex? x)))))
+ (complex (define (rational? x) (and (number? x) (not (%complex? x)))))
  (else (define rational? number?)))
 (define real? rational?)
 
@@ -986,12 +986,16 @@
 
 (cond-expand
  (complex
-  (define (real-part z) (if (complex? z) (complex-real z) z))
+  (define (real-part z) (if (%complex? z) (complex-real z) z))
-  (define (imag-part z) (if (complex? z) (complex-imag z) 0.0))
+  (define (imag-part z) (if (%complex? z) (complex-imag z) 0.0))
   (define (magnitude z)
     (sqrt (+ (* (real-part z) (real-part z))
              (* (imag-part z) (imag-part z)))))
-  (define (angle z) (atan (imag-part z) (real-part z))))
+  (define (angle z) (atan (imag-part z) (real-part z)))
+  (define (make-rectangular x y)
+    (+ x (* z 0+1i)))
+  (define (make-polar r phi)
+    (make-rectangular (* r (cos phi)) (* r (sin phi)))))
  (else
   (define (real-part z) z)
   (define (imag-part z) 0.0)

opcodes.c

@@ -88,7 +88,7 @@ _FN1(_I(SEXP_FIXNUM), _I(SEXP_RATIO), "ratio-numerator", 0, sexp_ratio_numerator
 _FN1(_I(SEXP_FIXNUM), _I(SEXP_RATIO), "ratio-denominator", 0, sexp_ratio_denominator_op),
 #endif
 #if SEXP_USE_COMPLEX
-_OP(SEXP_OPC_TYPE_PREDICATE, SEXP_OP_TYPEP,  1, 0, _I(SEXP_BOOLEAN), _I(SEXP_OBJECT), SEXP_FALSE, SEXP_FALSE, 0, "complex?", _I(SEXP_COMPLEX), 0),
+_OP(SEXP_OPC_TYPE_PREDICATE, SEXP_OP_TYPEP,  1, 0, _I(SEXP_BOOLEAN), _I(SEXP_OBJECT), SEXP_FALSE, SEXP_FALSE, 0, "%complex?", _I(SEXP_COMPLEX), 0),
 _FN1(_I(SEXP_FLONUM), _I(SEXP_RATIO), "complex-real", 0, sexp_complex_real_op),
 _FN1(_I(SEXP_FLONUM), _I(SEXP_RATIO), "complex-imag", 0, sexp_complex_imag_op),
 #endif

opt/bignum.c

@@ -423,7 +423,7 @@ sexp sexp_bignum_quot_rem (sexp ctx, sexp *rem, sexp a, sexp b) {
   res = quot_step(ctx, rem, a1, b1, k, i);
   sexp_bignum_sign(res) = sexp_bignum_sign(a) * sexp_bignum_sign(b);
   if (sexp_bignum_sign(a) < 0) {
-    sexp_negate(*rem);
+    sexp_negate_exact(*rem);
   }
   sexp_gc_release4(ctx);
   return res;
@@ -529,6 +529,17 @@ sexp sexp_complex_add (sexp ctx, sexp a, sexp b) {
   return sexp_complex_normalize(res);
 }
 
+sexp sexp_complex_sub (sexp ctx, sexp a, sexp b) {
+  sexp_gc_var2(res, tmp);
+  sexp_gc_preserve2(ctx, res, tmp);
+  tmp = sexp_make_complex(ctx, sexp_complex_real(b), sexp_complex_imag(b));
+  sexp_negate(sexp_complex_real(tmp));
+  sexp_negate(sexp_complex_imag(tmp));
+  res = sexp_complex_add(ctx, a, tmp);
+  sexp_gc_release2(ctx);
+  return res;
+}
+
 sexp sexp_complex_mul (sexp ctx, sexp a, sexp b) {
   sexp_gc_var3(res, real, imag);
   sexp_gc_preserve3(ctx, res, real, imag);
@@ -569,6 +580,121 @@ sexp sexp_complex_div (sexp ctx, sexp a, sexp b) {
   return sexp_complex_normalize(res);
 }
 
+#if SEXP_USE_MATH
+
+static double sexp_to_double (sexp x) {
+  if (sexp_flonump(x))
+    return sexp_flonum_value(x);
+  else if (sexp_fixnump(x))
+    return sexp_fixnum_to_double(x);
+  else if (sexp_bignump(x))
+    return sexp_bignum_to_double(x);
+#if SEXP_USE_RATIOS
+  else if (sexp_ratiop(x))
+    return sexp_ratio_to_double(x);
+#endif
+  else
+    return 0.0;
+}
+
+sexp sexp_complex_math_error (sexp ctx, sexp z) {
+  return sexp_type_exception(ctx, NULL, SEXP_NUMBER, z);
+}
+
+sexp sexp_complex_sqrt (sexp ctx, sexp z) {
+  double x = sexp_to_double(sexp_complex_real(z)),
+    y = sexp_to_double(sexp_complex_imag(z)), r;
+  sexp_gc_var1(res);
+  sexp_gc_preserve1(ctx, res);
+  r = sqrt(x*x + y*y);
+  res = sexp_make_complex(ctx, SEXP_ZERO, SEXP_ZERO);
+  sexp_complex_real(res) = sexp_make_flonum(ctx, sqrt((x+r)/2));
+  sexp_complex_imag(res) = sexp_make_flonum(ctx, (y<0?-1:1)*sqrt((-x+r)/2));
+  sexp_gc_release1(ctx);
+  return res;
+}
+
+sexp sexp_complex_exp (sexp ctx, sexp z) {
+  double x = sexp_to_double(sexp_complex_real(z)),
+    y = sexp_to_double(sexp_complex_imag(z));
+  sexp_gc_var1(res);
+  sexp_gc_preserve1(ctx, res);
+  res = sexp_make_complex(ctx, SEXP_ZERO, SEXP_ZERO);
+  sexp_complex_real(res) = sexp_make_flonum(ctx, exp(x)*cos(y));
+  sexp_complex_imag(res) = sexp_make_flonum(ctx, sin(y));
+  sexp_gc_release1(ctx);
+  return res;
+}
+
+sexp sexp_complex_log (sexp ctx, sexp z) {
+  double x = sexp_to_double(sexp_complex_real(z)),
+    y = sexp_to_double(sexp_complex_imag(z));
+  sexp_gc_var1(res);
+  sexp_gc_preserve1(ctx, res);
+  res = sexp_make_complex(ctx, SEXP_ZERO, SEXP_ZERO);
+  sexp_complex_real(res) = sexp_make_flonum(ctx, log(sqrt(x*x + y*y)));
+  sexp_complex_imag(res) = sexp_make_flonum(ctx, atan2(y, x));
+  sexp_gc_release1(ctx);
+  return res;
+}
+
+sexp sexp_complex_sin (sexp ctx, sexp z) {
+  double x = sexp_to_double(sexp_complex_real(z)),
+    y = sexp_to_double(sexp_complex_imag(z));
+  sexp_gc_var1(res);
+  sexp_gc_preserve1(ctx, res);
+  res = sexp_make_complex(ctx, SEXP_ZERO, SEXP_ZERO);
+  sexp_complex_real(res) = sexp_make_flonum(ctx, sin(x)*cosh(y));
+  sexp_complex_imag(res) = sexp_make_flonum(ctx, cos(x)*sinh(y));
+  sexp_gc_release1(ctx);
+  return res;
+}
+
+sexp sexp_complex_cos (sexp ctx, sexp z) {
+  double x = sexp_to_double(sexp_complex_real(z)),
+    y = sexp_to_double(sexp_complex_imag(z));
+  sexp_gc_var1(res);
+  sexp_gc_preserve1(ctx, res);
+  res = sexp_make_complex(ctx, SEXP_ZERO, SEXP_ZERO);
+  sexp_complex_real(res) = sexp_make_flonum(ctx, cos(x)*cosh(y));
+  sexp_complex_imag(res) = sexp_make_flonum(ctx, -sin(x)*sinh(y));
+  sexp_gc_release1(ctx);
+  return res;
+}
+
+sexp sexp_complex_asin (sexp ctx, sexp z) {
+  sexp_gc_var2(res, tmp);
+  sexp_gc_preserve2(ctx, res, tmp);
+  res = sexp_complex_mul(ctx, z, z);
+  tmp = sexp_make_complex(ctx, SEXP_ONE, SEXP_ZERO);
+  res = sexp_complex_sub(ctx, tmp, res);
+  res = sexp_complex_sqrt(ctx, res);
+  /* tmp = iz */
+  sexp_complex_real(tmp) = sexp_complex_imag(z);
+  sexp_negate(sexp_complex_real(tmp));
+  sexp_complex_imag(tmp) = sexp_complex_real(z);
+  res = sexp_complex_add(ctx, tmp, res);
+  tmp = sexp_complex_log(ctx, res);
+  /* res = -i*tmp */
+  sexp_complex_real(res) = sexp_complex_imag(tmp);
+  sexp_complex_imag(res) = sexp_complex_real(tmp);
+  sexp_negate(sexp_complex_imag(res));
+  sexp_gc_release2(ctx);
+  return res;
+}
+
+sexp sexp_complex_acos (sexp ctx, sexp z) {
+  sexp_gc_var2(res, tmp);
+  sexp_gc_preserve2(ctx, res, tmp);
+  res = sexp_complex_asin(ctx, z);
+  tmp = sexp_make_complex(ctx, SEXP_ZERO, SEXP_ZERO);
+  sexp_complex_real(tmp) = sexp_make_flonum(ctx, acos(-1)/2);
+  res = sexp_sub(ctx, tmp, res);
+  sexp_gc_release2(ctx);
+  return res;
+}
+
+#endif
 #endif
 
 /****************** generic arithmetic ************************/
@@ -757,7 +883,7 @@ sexp sexp_sub (sexp ctx, sexp a, sexp b) {
   case SEXP_NUM_FIX_BIG:
     tmp1 = sexp_fixnum_to_bignum(ctx, a);
     r = sexp_bignum_sub(ctx, NULL, b, tmp1);
-    sexp_negate(r);
+    sexp_negate_exact(r);
     r = sexp_bignum_normalize(r);
     break;
   case SEXP_NUM_FLO_FIX:
@@ -797,13 +923,13 @@ sexp sexp_sub (sexp ctx, sexp a, sexp b) {
     /* ... FALLTHROUGH ... */
   case SEXP_NUM_RAT_RAT:
     tmp2 = sexp_make_ratio(ctx, sexp_ratio_numerator(b), sexp_ratio_denominator(b));
-    sexp_negate(sexp_ratio_numerator(tmp2));
+    sexp_negate_exact(sexp_ratio_numerator(tmp2));
     r = sexp_ratio_add(ctx, a, tmp2);
     if (negatep) {
       if (sexp_ratiop(r)) {
-        sexp_negate(sexp_ratio_numerator(r));
+        sexp_negate_exact(sexp_ratio_numerator(r));
       } else {
-        sexp_negate(r);
+        sexp_negate_exact(r);
       }
     }
     break;
@@ -828,10 +954,7 @@ sexp sexp_sub (sexp ctx, sexp a, sexp b) {
     /* ... FALLTHROUGH ... */
   case SEXP_NUM_CPX_CPX:
   complex_sub:
-    tmp2 = sexp_make_complex(ctx, sexp_complex_real(b), sexp_complex_imag(b));
+    r = sexp_complex_sub(ctx, a, b);
-    sexp_negate(sexp_complex_real(tmp2));
-    sexp_negate(sexp_complex_imag(tmp2));
-    r = sexp_complex_add(ctx, a, tmp2);
     if (negatep) {
       if (sexp_complexp(r)) {
         sexp_negate(sexp_complex_real(r));