improve atan, asin, and acos on complex numbers

Replace naive calculations with ones based on Kahan's "Branch Cuts for
Complex Elementary Functions" as implemented in Chez Scheme.
This commit is contained in:
Matthew Flatt 2019-06-28 14:02:30 -06:00
parent 8e1b27592f
commit 2847d1d22a
4 changed files with 244 additions and 87 deletions

View File

@ -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?)
(when has-exact-zero-inexact-complex?
(test -inf.0 atan 0+1.0i)
(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))))

View File

@ -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);
}

View File

@ -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[])

View File

@ -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;