diff --git a/pkgs/racket-test-core/tests/racket/number.rktl b/pkgs/racket-test-core/tests/racket/number.rktl index da1ecb72aa..78c6a53791 100644 --- a/pkgs/racket-test-core/tests/racket/number.rktl +++ b/pkgs/racket-test-core/tests/racket/number.rktl @@ -2144,8 +2144,54 @@ (err/rt-test (atan 0 0) exn:fail:contract:divide-by-zero?) (err/rt-test (atan 0+i) exn:fail:contract:divide-by-zero?) (err/rt-test (atan 0-i) exn:fail:contract:divide-by-zero?) -(test -inf.0 atan 0+1.0i) -(test -inf.0 atan 0-1.0i) + +(when has-exact-zero-inexact-complex? + (test -inf.0 atan 0+1.0i) + (test -inf.0 atan 0-1.0i)) + +(test 79.0+17710.0i z-round (* 100 (atan 0.0+1.0i))) +(test 79.0-17710.0i z-round (* 100 (atan 0.0-1.0i))) +(test 157.0+55.0i z-round (* 100 (atan 0.0+2.0i))) +(test 157.0-55.0i z-round (* 100 (atan 0.0-2.0i))) +(test -157.0+55.0i z-round (* 100 (atan -0.0+2.0i))) +(test 134.0-40.0i z-round (* 100 (atan 1.0-2.0i))) +(test 157.0+0.0i z-round (* 100 (atan 1.0-4e153i))) +(test 157.0+0.0i z-round (* 100 (atan 1.0+4e153i))) +(test -2.5e-154 imag-part (atan 1.0-4e153i)) +(test +2.5e-154 imag-part (atan 1.0+4e153i)) +(test 157.0+0.0i z-round (* 100 (atan 5e153+4e153i))) +(test 125.0 round (* 1e156 (imag-part (atan 5e153+4e153i)))) +(test 157.0+0.0i z-round (* 100 (atan 4e153+4e153i))) +(test 125.0 round (* 1e156 (imag-part (atan 4e153+4e153i)))) + +#reader "maybe-single.rkt" +(when has-single-flonum? + (test 79.0f0+17710.0f0i z-round (* 100 (atan 0.0f0+1.0f0i)))) + +(test 157.0-inf.0i z-round (* 100 (asin +inf.0))) +(test -157.0+inf.0i z-round (* 100 (asin -inf.0))) +(test 0.0+inf.0i z-round (* 100 (acos +inf.0))) +(test 314.0-inf.0i z-round (* 100 (acos -inf.0))) +(test 0.0-inf.0i asin -inf.0i) +(test 157.0+inf.0i z-round (* 100 (acos -inf.0i))) + +(test 63.0+231.0i z-round (* 100 (asin 3+4i))) +(test 63.0+231.0i z-round (* 100 (asin 3.0+4.0i))) +(test 94.0-231.0i z-round (* 100 (acos 3+4i))) +(test 94.0-231.0i z-round (* 100 (acos 3.0+4.0i))) + +(test 157.0-46300.0i z-round (* 100 (asin 6e200))) +(test 0.0+46300.0i z-round (* 100 (acos 6e200))) + +#reader "maybe-single.rkt" +(when has-single-flonum? + (test 157.0f0-inf.fi z-round (* 100 (asin +inf.f))) + (test -157.0f0+inf.fi z-round (* 100 (asin -inf.f))) + (test 0.0f0+inf.fi z-round (* 100 (acos +inf.f))) + (test 314.0f0-inf.fi z-round (* 100 (acos -inf.f))) + (test 63.0f0+231.0f0i z-round (* 100 (asin 3.0f0+4.0f0i))) + (test 94.0f0-231.0f0i z-round (* 100 (acos 3.0f0+4.0f0i)))) + (test 1024.0 round (expt 2.0 10.0)) (test 1024.0 round (expt -2.0 10.0)) (test -512.0 round (expt -2.0 9.0)) @@ -2250,8 +2296,6 @@ (test-inf-bad tan) (test-inf-bad sin) (test-inf-bad cos) -(test-inf-bad asin) -(test-inf-bad acos) (test 11/7 rationalize (inexact->exact (atan +inf.0 1)) 1/100) (test -11/7 rationalize (inexact->exact (atan -inf.0 1)) 1/100) @@ -2303,7 +2347,7 @@ (test 693.+3142.i z-round (* 1000 (log -2))) (test 1571.-1317.i z-round (* 1000 (asin 2))) (test -1571.+1317.i z-round (* 1000 (asin -2))) -(test 0+3688.i z-round (* 1000 (acos 20))) +(test 0.0+3688.i z-round (* 1000 (acos 20))) (test 3142.-3688.i z-round (* 1000 (acos -20))) (define (cs2 c) (+ (* (cos c) (cos c)) (* (sin c) (sin c)))) diff --git a/racket/src/racket/src/complex.c b/racket/src/racket/src/complex.c index aefe409d45..b132008a4f 100644 --- a/racket/src/racket/src/complex.c +++ b/racket/src/racket/src/complex.c @@ -402,3 +402,138 @@ Scheme_Object *scheme_complex_sqrt(const Scheme_Object *o) return scheme_make_complex(ni, nr); } + +Scheme_Object *scheme_complex_atan(const Scheme_Object *c) +{ + /* From Chez Scheme, which implements the principal expression from Kahan: + (log(1+z) - log(1-z))/2 + */ +# define OMEGA 1.7976931348623157e308 +# define THETA 3.351951982485649e+153 /* = sqrt(OMEGA)/4 */ +# define RHO 2.9833362924800834e-154 /* = 1/THETA */ +# define HALF_PI 1.5707963267948966 +# define IS_NEG(x) (((x) < 0.0) || ((x == 0.0) && scheme_minus_zero_p(x))) + double x, y, ay, r, i; + int negate; + + /* Get components after multiplication by i */ + x = -scheme_real_to_double(_scheme_complex_imaginary_part(c)); + y = scheme_real_to_double(_scheme_complex_real_part(c)); + + /* Compute atanh */ + + if (x < 0.0) { + x = -x; + negate = 1; + } else { + y = -y; + negate = 0; + } + + ay = fabs(y); + + if ((x > THETA) || (y > THETA)) { + /* RP(1/z) +/- (pi/2)i */ + if (x > ay) + r = 1/(x + ((y/x) * y)); + else if (x < ay) { + r = y/x; + r = r/((x * r)+ y); + } else + r = 1/(x+ay); + + i = (IS_NEG(y) ? HALF_PI : (-HALF_PI)); + } else if (x == 1.0) { + double k = ay + RHO; + r = log(sqrt(sqrt((y * y) + 4.0)) / sqrt(k)); + i = (HALF_PI + atan(k/2.0)) / (IS_NEG(y) ? 2.0 : -2.0); + } else { + double mx = 1.0 - x; + double k = ay + RHO; + k = k * k; + + r = log(((4.0 * x) / ((mx * mx) + k)) + 1.0) / 4.0; + i = atan2(2.0 * y, (mx * (1.0 + x)) - k) / -2.0; + } + + if (negate) { + i = -i; + r = -r; + } + + /* Multiply by -i to get atan */ + x = i; + y = -r; + +#ifdef MZ_USE_SINGLE_FLOATS + if (SCHEME_FLTP(_scheme_complex_real_part(c)) + || SCHEME_FLTP(_scheme_complex_imaginary_part(c))) { + return scheme_make_complex(scheme_make_float(x), scheme_make_float(y)); + } +#endif + + return scheme_make_complex(scheme_make_double(x), scheme_make_double(y)); +} + +Scheme_Object *scheme_complex_asin_or_acos(const Scheme_Object *z, int get_asin) +{ + /* From Chez Scheme, which implements the principal expression from Kahan */ + Scheme_Object *zp, *zm, *aa[1]; + double a, b, c, d, r, i; + + aa[0] = scheme_bin_minus(scheme_make_integer(1), z); + zm = scheme_sqrt(1, aa); + + aa[0] = scheme_bin_plus(scheme_make_integer(1), z); + zp = scheme_sqrt(1, aa); + + if (SCHEME_COMPLEXP(zm)) { + a = scheme_real_to_double(_scheme_complex_real_part(zm)); + b = scheme_real_to_double(_scheme_complex_imaginary_part(zm)); + } else { + a = scheme_real_to_double(zm); + b = 0.0; + } + + if (SCHEME_COMPLEXP(zp)) { + c = scheme_real_to_double(_scheme_complex_real_part(zp)); + d = scheme_real_to_double(_scheme_complex_imaginary_part(zp)); + } else { + c = scheme_real_to_double(zp); + d = 0.0; + } + + if (get_asin) { + if (SCHEME_COMPLEXP(z)) { + r = scheme_real_to_double(_scheme_complex_real_part(z)); + r = atan2(r, (a*c)-(b*d)); + } else { + r = scheme_real_to_double((Scheme_Object *)z); + r = atan2(r, 0.0); /* void +nan.0 from (a*c)-(b*d) */ + } + + i = asinh((a*d)-(b*c)); + } else { + r = 2.0 * atan2(a, c); + i = asinh((b*c) - (a*d)); + } + +#ifdef MZ_USE_SINGLE_FLOATS + if (SCHEME_FLTP(_scheme_complex_real_part(z)) + || SCHEME_FLTP(_scheme_complex_imaginary_part(z))) { + return scheme_make_complex(scheme_make_float(r), scheme_make_float(i)); + } +#endif + + return scheme_make_complex(scheme_make_double(r), scheme_make_double(i)); +} + +Scheme_Object *scheme_complex_asin(const Scheme_Object *c) +{ + return scheme_complex_asin_or_acos(c, 1); +} + +Scheme_Object *scheme_complex_acos(const Scheme_Object *c) +{ + return scheme_complex_asin_or_acos(c, 0); +} diff --git a/racket/src/racket/src/number.c b/racket/src/racket/src/number.c index a54d58a5c5..82880bb000 100644 --- a/racket/src/racket/src/number.c +++ b/racket/src/racket/src/number.c @@ -224,9 +224,10 @@ READ_ONLY Scheme_Object *scheme_inf_object, *scheme_minus_inf_object, *scheme_na #define zeroi scheme_exact_zero -READ_ONLY Scheme_Object *scheme_zerod, *scheme_nzerod, *scheme_pi, *scheme_half_pi, *scheme_plus_i, *scheme_minus_i; +READ_ONLY Scheme_Object *scheme_zerod, *scheme_nzerod, *scheme_pi, *scheme_half_pi, *scheme_minus_half_pi; +READ_ONLY Scheme_Object *scheme_plus_i, *scheme_minus_i; #ifdef MZ_USE_SINGLE_FLOATS -READ_ONLY Scheme_Object *scheme_zerof, *scheme_nzerof, *scheme_single_pi, *scheme_single_half_pi; +READ_ONLY Scheme_Object *scheme_zerof, *scheme_nzerof, *scheme_single_pi, *scheme_single_half_pi, *scheme_single_minus_half_pi; READ_ONLY Scheme_Object *scheme_single_inf_object, *scheme_single_minus_inf_object, *scheme_single_nan_object; #endif @@ -330,11 +331,13 @@ scheme_init_number (Scheme_Startup_Env *env) REGISTER_SO(scheme_pi); REGISTER_SO(scheme_half_pi); + REGISTER_SO(scheme_minus_half_pi); REGISTER_SO(scheme_zerod); REGISTER_SO(scheme_nzerod); #ifdef MZ_USE_SINGLE_FLOATS REGISTER_SO(scheme_single_pi); REGISTER_SO(scheme_single_half_pi); + REGISTER_SO(scheme_single_minus_half_pi); REGISTER_SO(scheme_zerof); REGISTER_SO(scheme_nzerof); #endif @@ -390,15 +393,16 @@ scheme_init_number (Scheme_Startup_Env *env) scheme_pi = scheme_make_double(atan2(0.0, -1.0)); scheme_half_pi = scheme_make_double(atan2(0.0, -1.0)/2); + scheme_minus_half_pi = scheme_make_double(-SCHEME_DBL_VAL(scheme_half_pi)); #ifdef MZ_USE_SINGLE_FLOATS scheme_zerof = scheme_make_float(0.0f); scheme_nzerof = scheme_make_float(-0.0f); scheme_single_pi = scheme_make_float(SCHEME_DBL_VAL(scheme_pi)); scheme_single_half_pi = scheme_make_float(SCHEME_DBL_VAL(scheme_half_pi)); + scheme_single_minus_half_pi = scheme_make_float(SCHEME_DBL_VAL(scheme_minus_half_pi)); #endif scheme_plus_i = scheme_make_complex(scheme_make_integer(0), scheme_make_integer(1)); scheme_minus_i = scheme_make_complex(scheme_make_integer(0), scheme_make_integer(-1)); - scheme_inf_object = scheme_make_double(scheme_infinity_val); scheme_minus_inf_object = scheme_make_double(scheme_minus_infinity_val); #ifdef NAN_EQUALS_ANYTHING @@ -2762,18 +2766,6 @@ static Scheme_Object *get_frac(char *name, int low_p, return n; } -#ifdef USE_SINGLE_FLOATS -XFORM_NONGCING static int complex_single_flonum_p(Scheme_Object *c) -{ - Scheme_Complex *cb = (Scheme_Complex *)c; - - if (SCHEME_FLTP(cb->r) || SCHEME_FLTP(cb->i)) - return 1; - - return 0; -} -#endif - static Scheme_Object *un_exp(Scheme_Object *o); static Scheme_Object *un_log(Scheme_Object *o); @@ -2853,8 +2845,6 @@ static Scheme_Object *bignum_log(Scheme_Object *b) return scheme_make_double(d); } -static Scheme_Object *complex_sin(Scheme_Object *c); - static Scheme_Object *complex_sin(Scheme_Object *c) { Scheme_Object *i_c; @@ -2866,8 +2856,6 @@ static Scheme_Object *complex_sin(Scheme_Object *c) scheme_bin_mult(scheme_make_integer(2), scheme_plus_i)); } -static Scheme_Object *complex_cos(Scheme_Object *c); - static Scheme_Object *complex_cos(Scheme_Object *c) { Scheme_Object *i_c; @@ -2879,65 +2867,23 @@ static Scheme_Object *complex_cos(Scheme_Object *c) scheme_make_integer(2)); } -static Scheme_Object *complex_tan(Scheme_Object *c); - static Scheme_Object *complex_tan(Scheme_Object *c) { return scheme_bin_div(complex_sin(c), complex_cos(c)); } -static Scheme_Object *complex_asin(Scheme_Object *c); -static Scheme_Object *complex_atan(Scheme_Object *c); - static Scheme_Object *complex_asin(Scheme_Object *c) { - Scheme_Object *one_minus_c_sq, *sqrt_1_minus_c_sq; - - one_minus_c_sq = scheme_bin_minus(scheme_make_integer(1), - scheme_bin_mult(c, c)); - sqrt_1_minus_c_sq = scheme_sqrt(1, &one_minus_c_sq); - return scheme_bin_mult(scheme_make_integer(2), - complex_atan(scheme_bin_div(c, - scheme_bin_plus(scheme_make_integer(1), - sqrt_1_minus_c_sq)))); + return scheme_complex_asin(c); } -static Scheme_Object *complex_acos(Scheme_Object *c); - static Scheme_Object *complex_acos(Scheme_Object *c) { - Scheme_Object *a, *r; - a = complex_asin(c); - if (scheme_is_zero(_scheme_complex_imaginary_part(c)) - && (scheme_bin_gt(_scheme_complex_real_part(c), scheme_make_integer(1)) - || scheme_bin_lt(_scheme_complex_real_part(c), scheme_make_integer(-1)))) { - /* Make sure real part is 0 or pi */ - if (scheme_is_negative(_scheme_complex_real_part(c))) { -#ifdef MZ_USE_SINGLE_FLOATS - if (complex_single_flonum_p(c)) - r = scheme_single_pi; - else -#endif - r = scheme_pi; - } else - r = scheme_make_integer(0); - return scheme_make_complex(r, scheme_bin_minus(scheme_make_integer(0), - _scheme_complex_imaginary_part(a))); - } else { -#ifdef MZ_USE_SINGLE_FLOATS - if (complex_single_flonum_p(c)) - r = scheme_single_half_pi; - else -#endif - r = scheme_half_pi; - return scheme_bin_minus(r, a); - } + return scheme_complex_acos(c); } static Scheme_Object *complex_atan(Scheme_Object *c) { - Scheme_Object *one_half = NULL; - if (SAME_OBJ(_scheme_complex_real_part(c), scheme_make_integer(0))) { Scheme_Object *i = _scheme_complex_imaginary_part(c); if (SAME_OBJ(i, scheme_make_integer(1)) || SAME_OBJ(i, scheme_make_integer(-1))) { @@ -2956,19 +2902,7 @@ static Scheme_Object *complex_atan(Scheme_Object *c) } } - /* select single versus complex: */ -#ifdef MZ_USE_SINGLE_FLOATS - if (complex_single_flonum_p(c)) - one_half = scheme_make_float(0.5); - else -#endif - one_half = scheme_make_double(0.5); - - return scheme_bin_mult(scheme_plus_i, - scheme_bin_mult(one_half, - un_log(scheme_bin_div(scheme_bin_plus(scheme_plus_i, c), - scheme_bin_plus(scheme_plus_i, - scheme_bin_minus(zeroi, c)))))); + return scheme_complex_atan(c); } #define GEN_ZERO_IS_ZERO() if (o == zeroi) return zeroi; @@ -3060,25 +2994,65 @@ long_double scheme_long_double_exp(long_double x) { return long_double_exp(x); } #endif -static Scheme_Object *scheme_inf_plus_pi() +static Scheme_Object *inf_plus_pi_i() { return scheme_make_complex(scheme_inf_object, scheme_pi); } +static Scheme_Object *half_pi_minus_inf_i() +{ + return scheme_make_complex(scheme_half_pi, scheme_minus_inf_object); +} + +static Scheme_Object *minus_half_pi_plus_inf_i() +{ + return scheme_make_complex(scheme_minus_half_pi, scheme_inf_object); +} + +static Scheme_Object *zero_plus_inf_i() +{ + return scheme_make_complex(scheme_zerod, scheme_inf_object); +} + +static Scheme_Object *pi_minus_inf_i() +{ + return scheme_make_complex(scheme_pi, scheme_minus_inf_object); +} + #ifdef MZ_USE_SINGLE_FLOATS -static Scheme_Object *scheme_single_inf_plus_pi() +static Scheme_Object *single_inf_plus_pi_i() { return scheme_make_complex(scheme_single_inf_object, scheme_single_pi); } + +static Scheme_Object *single_half_pi_minus_inf_i() +{ + return scheme_make_complex(scheme_single_half_pi, scheme_single_minus_inf_object); +} + +static Scheme_Object *single_minus_half_pi_plus_inf_i() +{ + return scheme_make_complex(scheme_single_minus_half_pi, scheme_single_inf_object); +} + +static Scheme_Object *single_zero_plus_inf_i() +{ + return scheme_make_complex(scheme_zerof, scheme_single_inf_object); +} + +static Scheme_Object *single_pi_minus_inf_i() +{ + return scheme_make_complex(scheme_single_pi, scheme_single_minus_inf_object); +} #endif GEN_UNARY_OP(exp_prim, exp, exp, scheme_inf_object, scheme_single_inf_object, scheme_zerod, scheme_zerof, scheme_nan_object, scheme_single_nan_object, complex_exp, GEN_ZERO_IS_ONE, NEVER_RESORT_TO_COMPLEX, BIGNUMS_AS_DOUBLES) -GEN_UNARY_OP(log_e_prim, log, SCH_LOG, scheme_inf_object, scheme_single_inf_object, scheme_inf_plus_pi(), scheme_single_inf_plus_pi(), scheme_nan_object, scheme_single_nan_object, complex_log, GEN_ONE_IS_ZERO_AND_ZERO_IS_ERR, NEGATIVE_USES_COMPLEX, BIGNUM_LOG) +GEN_UNARY_OP(log_e_prim, log, SCH_LOG, scheme_inf_object, scheme_single_inf_object, inf_plus_pi_i(), single_inf_plus_pi_i(), scheme_nan_object, scheme_single_nan_object, complex_log, GEN_ONE_IS_ZERO_AND_ZERO_IS_ERR, NEGATIVE_USES_COMPLEX, BIGNUM_LOG) GEN_UNARY_OP(sin_prim, sin, SCH_SIN, scheme_nan_object, scheme_single_nan_object, scheme_nan_object, scheme_single_nan_object, scheme_nan_object, scheme_single_nan_object, complex_sin, GEN_ZERO_IS_ZERO, NEVER_RESORT_TO_COMPLEX, BIGNUMS_AS_DOUBLES) GEN_UNARY_OP(cos_prim, cos, SCH_COS, scheme_nan_object, scheme_single_nan_object, scheme_nan_object, scheme_single_nan_object, scheme_nan_object, scheme_single_nan_object, complex_cos, GEN_ZERO_IS_ONE, NEVER_RESORT_TO_COMPLEX, BIGNUMS_AS_DOUBLES) GEN_UNARY_OP(tan_prim, tan, SCH_TAN, scheme_nan_object, scheme_single_nan_object, scheme_nan_object, scheme_single_nan_object, scheme_nan_object, scheme_single_nan_object, complex_tan, GEN_ZERO_IS_ZERO, NEVER_RESORT_TO_COMPLEX, BIGNUMS_AS_DOUBLES) -GEN_UNARY_OP(asin_prim, asin, SCH_ASIN, scheme_nan_object, scheme_single_nan_object, scheme_nan_object, scheme_single_nan_object, scheme_nan_object, scheme_single_nan_object, complex_asin, GEN_ZERO_IS_ZERO, OVER_ONE_MAG_USES_COMPLEX, BIGNUMS_AS_DOUBLES) -GEN_UNARY_OP(acos_prim, acos, acos, scheme_nan_object, scheme_single_nan_object, scheme_nan_object, scheme_single_nan_object, scheme_nan_object, scheme_single_nan_object, complex_acos, GEN_ONE_IS_ZERO, OVER_ONE_MAG_USES_COMPLEX, BIGNUMS_AS_DOUBLES) +GEN_UNARY_OP(asin_prim, asin, SCH_ASIN, half_pi_minus_inf_i(), single_half_pi_minus_inf_i(), minus_half_pi_plus_inf_i(), single_minus_half_pi_plus_inf_i(), scheme_nan_object, scheme_single_nan_object, complex_asin, GEN_ZERO_IS_ZERO, OVER_ONE_MAG_USES_COMPLEX, BIGNUMS_AS_DOUBLES) +GEN_UNARY_OP(acos_prim, acos, acos, zero_plus_inf_i(), single_zero_plus_inf_i(), pi_minus_inf_i(), single_pi_minus_inf_i(), scheme_nan_object, scheme_single_nan_object, complex_acos, GEN_ONE_IS_ZERO, OVER_ONE_MAG_USES_COMPLEX, BIGNUMS_AS_DOUBLES) static Scheme_Object * log_prim (int argc, Scheme_Object *argv[]) diff --git a/racket/src/racket/src/schpriv.h b/racket/src/racket/src/schpriv.h index 1ad348e48a..759bfcb47b 100644 --- a/racket/src/racket/src/schpriv.h +++ b/racket/src/racket/src/schpriv.h @@ -2300,6 +2300,9 @@ Scheme_Object *scheme_complex_multiply(const Scheme_Object *a, const Scheme_Obje Scheme_Object *scheme_complex_divide(const Scheme_Object *n, const Scheme_Object *d); Scheme_Object *scheme_complex_power(const Scheme_Object *a, const Scheme_Object *b); Scheme_Object *scheme_complex_sqrt(const Scheme_Object *a); +Scheme_Object *scheme_complex_atan(const Scheme_Object *c); +Scheme_Object *scheme_complex_asin(const Scheme_Object *c); +Scheme_Object *scheme_complex_acos(const Scheme_Object *c); XFORM_NONGCING int scheme_is_complex_exact(const Scheme_Object *o); /****** Inexacts ******/ @@ -2415,7 +2418,8 @@ extern int scheme_is_nan(double); extern double scheme_infinity_val, scheme_minus_infinity_val; extern double scheme_floating_point_zero; extern double scheme_floating_point_nzero; -extern Scheme_Object *scheme_zerod, *scheme_nzerod, *scheme_pi, *scheme_half_pi, *scheme_plus_i, *scheme_minus_i; +extern Scheme_Object *scheme_zerod, *scheme_nzerod, *scheme_pi, *scheme_half_pi, *scheme_minus_half_pi; +extern Scheme_Object *scheme_plus_i, *scheme_minus_i; extern Scheme_Object *scheme_inf_object, *scheme_minus_inf_object, *scheme_nan_object; #ifdef MZ_LONG_DOUBLE extern long_double scheme_long_infinity_val, scheme_long_minus_infinity_val;