diff --git a/eval.c b/eval.c index 4d3dfe5c..09abed32 100644 --- a/eval.c +++ b/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_log, log) -define_math_op(sexp_sin, sin) -define_math_op(sexp_cos, cos) -define_math_op(sexp_tan, tan) -define_math_op(sexp_asin, asin) -define_math_op(sexp_acos, acos) -define_math_op(sexp_atan, atan) -define_math_op(sexp_round, round) -define_math_op(sexp_trunc, trunc) -define_math_op(sexp_floor, floor) -define_math_op(sexp_ceiling, ceil) +define_math_op(sexp_exp, exp, sexp_complex_exp) +define_math_op(sexp_log, log, sexp_complex_log) +define_math_op(sexp_sin, sin, sexp_complex_sin) +define_math_op(sexp_cos, cos, sexp_complex_cos) +define_math_op(sexp_tan, tan, sexp_complex_math_error) +define_math_op(sexp_asin, asin, sexp_complex_asin) +define_math_op(sexp_acos, acos, sexp_complex_acos) +define_math_op(sexp_atan, atan, sexp_complex_math_error) +define_math_op(sexp_round, round, sexp_complex_math_error) +define_math_op(sexp_trunc, trunc, sexp_complex_math_error) +define_math_op(sexp_floor, floor, sexp_complex_math_error) +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); diff --git a/include/chibi/bignum.h b/include/chibi/bignum.h index 8d75fe28..82950294 100644 --- a/include/chibi/bignum.h +++ b/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 */ diff --git a/include/chibi/sexp.h b/include/chibi/sexp.h index 480280b1..68b355be 100755 --- a/include/chibi/sexp.h +++ b/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])) diff --git a/lib/init.scm b/lib/init.scm index 74e7386d..0359747b 100644 --- a/lib/init.scm +++ b/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 (imag-part z) (if (complex? z) (complex-imag z) 0.0)) + (define (real-part z) (if (%complex? z) (complex-real z) z)) + (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) diff --git a/opcodes.c b/opcodes.c index 8a3ae052..c257ff8b 100644 --- a/opcodes.c +++ b/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 diff --git a/opt/bignum.c b/opt/bignum.c index 47f8aa21..5b38389f 100644 --- a/opt/bignum.c +++ b/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)); - sexp_negate(sexp_complex_real(tmp2)); - sexp_negate(sexp_complex_imag(tmp2)); - r = sexp_complex_add(ctx, a, tmp2); + r = sexp_complex_sub(ctx, a, b); if (negatep) { if (sexp_complexp(r)) { sexp_negate(sexp_complex_real(r));