From b0117b1d5ae73916c6f0d289be1f693bb5f46824 Mon Sep 17 00:00:00 2001 From: mbkma Date: Thu, 5 Mar 2020 13:06:45 +0100 Subject: Port to GNU MPFR/MPC Library For further information please visit: https://www.mpfr.org/ http://www.multiprecision.org/mpc --- src/mp-trigonometric.c | 578 ++++++++++--------------------------------------- 1 file changed, 111 insertions(+), 467 deletions(-) (limited to 'src/mp-trigonometric.c') diff --git a/src/mp-trigonometric.c b/src/mp-trigonometric.c index 72aba0c..ca19ca4 100644 --- a/src/mp-trigonometric.c +++ b/src/mp-trigonometric.c @@ -15,614 +15,258 @@ #include #include "mp.h" -#include "mp-private.h" - -static MPNumber pi; -static gboolean have_pi = FALSE; - -static int -mp_compare_mp_to_int(const MPNumber *x, int i) -{ - MPNumber t; - mp_set_from_integer(i, &t); - return mp_compare_mp_to_mp(x, &t); -} - /* Convert x to radians */ void convert_to_radians(const MPNumber *x, MPAngleUnit unit, MPNumber *z) { - MPNumber t1, t2; + int i; switch(unit) { default: case MP_RADIANS: mp_set_from_mp(x, z); - break; + return; case MP_DEGREES: - mp_get_pi(&t1); - mp_multiply(x, &t1, &t2); - mp_divide_integer(&t2, 180, z); + i = 180; break; case MP_GRADIANS: - mp_get_pi(&t1); - mp_multiply(x, &t1, &t2); - mp_divide_integer(&t2, 200, z); + i = 200; break; } + mpfr_t scale; + mpfr_init2(scale, PRECISION); + mpfr_const_pi(scale, MPFR_RNDN); + mpfr_div_si(scale, scale, i, MPFR_RNDN); + mpc_mul_fr(z->num, x->num, scale, MPC_RNDNN); + mpfr_clear(scale); } - -void -mp_get_pi(MPNumber *z) -{ - if (mp_is_zero(&pi)) { - mp_set_from_string("3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679", 10, &pi); - have_pi = TRUE; - } - mp_set_from_mp(&pi, z); -} - - void convert_from_radians(const MPNumber *x, MPAngleUnit unit, MPNumber *z) { - MPNumber t1, t2; + int i; - switch (unit) { + switch(unit) { default: case MP_RADIANS: mp_set_from_mp(x, z); - break; + return; case MP_DEGREES: - mp_multiply_integer(x, 180, &t2); - mp_get_pi(&t1); - mp_divide(&t2, &t1, z); + i = 180; break; case MP_GRADIANS: - mp_multiply_integer(x, 200, &t2); - mp_get_pi(&t1); - mp_divide(&t2, &t1, z); + i = 200; break; } + mpfr_t scale; + mpfr_init2(scale, PRECISION); + mpfr_const_pi(scale, MPFR_RNDN); + mpfr_si_div(scale, i, scale, MPFR_RNDN); + mpc_mul_fr(z->num, x->num, scale, MPC_RNDNN); + mpfr_clear(scale); } -/* z = sin(x) -1 >= x >= 1, do_sin = 1 - * z = cos(x) -1 >= x >= 1, do_sin = 0 - */ -static void -mpsin1(const MPNumber *x, MPNumber *z, int do_sin) -{ - int i, b2; - MPNumber t1, t2; - - /* sin(0) = 0, cos(0) = 1 */ - if (mp_is_zero(x)) { - if (do_sin == 0) - mp_set_from_integer(1, z); - else - mp_set_from_integer(0, z); - return; - } - - mp_multiply(x, x, &t2); - if (mp_compare_mp_to_int(&t2, 1) > 0) { - mperr("*** ABS(X) > 1 IN CALL TO MPSIN1 ***"); - } - - if (do_sin == 0) { - mp_set_from_integer(1, &t1); - mp_set_from_integer(0, z); - i = 1; - } else { - mp_set_from_mp(x, &t1); - mp_set_from_mp(&t1, z); - i = 2; - } - - /* Taylor series */ - /* POWER SERIES LOOP. REDUCE T IF POSSIBLE */ - b2 = 2 * max(MP_BASE, 64); - do { - if (MP_T + t1.exponent <= 0) - break; - - /* IF I*(I+1) IS NOT REPRESENTABLE AS AN INTEGER, THE FOLLOWING - * DIVISION BY I*(I+1) HAS TO BE SPLIT UP. - */ - mp_multiply(&t2, &t1, &t1); - if (i > b2) { - mp_divide_integer(&t1, -i, &t1); - mp_divide_integer(&t1, i + 1, &t1); - } else { - mp_divide_integer(&t1, -i * (i + 1), &t1); - } - mp_add(&t1, z, z); - - i += 2; - } while (t1.sign != 0); - - if (do_sin == 0) - mp_add_integer(z, 1, z); -} - - -static void -mp_sin_real(const MPNumber *x, MPAngleUnit unit, MPNumber *z) -{ - int xs; - MPNumber x_radians; - - /* sin(0) = 0 */ - if (mp_is_zero(x)) { - mp_set_from_integer(0, z); - return; - } - - convert_to_radians(x, unit, &x_radians); - - xs = x_radians.sign; - mp_abs(&x_radians, &x_radians); - - /* USE MPSIN1 IF ABS(X) <= 1 */ - if (mp_compare_mp_to_int(&x_radians, 1) <= 0) { - mpsin1(&x_radians, z, 1); - } - /* FIND ABS(X) MODULO 2PI */ - else { - mp_get_pi(z); - mp_divide_integer(z, 4, z); - mp_divide(&x_radians, z, &x_radians); - mp_divide_integer(&x_radians, 8, &x_radians); - mp_fractional_component(&x_radians, &x_radians); - - /* SUBTRACT 1/2, SAVE SIGN AND TAKE ABS */ - mp_add_fraction(&x_radians, -1, 2, &x_radians); - xs = -xs * x_radians.sign; - if (xs == 0) { - mp_set_from_integer(0, z); - return; - } - - x_radians.sign = 1; - mp_multiply_integer(&x_radians, 4, &x_radians); - - /* IF NOT LESS THAN 1, SUBTRACT FROM 2 */ - if (x_radians.exponent > 0) - mp_add_integer(&x_radians, -2, &x_radians); - - if (mp_is_zero(&x_radians)) { - mp_set_from_integer(0, z); - return; - } - - x_radians.sign = 1; - mp_multiply_integer(&x_radians, 2, &x_radians); - - /* NOW REDUCED TO FIRST QUADRANT, IF LESS THAN PI/4 USE - * POWER SERIES, ELSE COMPUTE COS OF COMPLEMENT - */ - if (x_radians.exponent > 0) { - mp_add_integer(&x_radians, -2, &x_radians); - mp_multiply(&x_radians, z, &x_radians); - mpsin1(&x_radians, z, 0); - } else { - mp_multiply(&x_radians, z, &x_radians); - mpsin1(&x_radians, z, 1); - } - } - - z->sign = xs; -} - - -static void -mp_cos_real(const MPNumber *x, MPAngleUnit unit, MPNumber *z) +void +mp_get_pi (MPNumber *z) { - /* cos(0) = 1 */ - if (mp_is_zero(x)) { - mp_set_from_integer(1, z); - return; - } - - convert_to_radians(x, unit, z); - - /* Use power series if |x| <= 1 */ - mp_abs(z, z); - if (mp_compare_mp_to_int(z, 1) <= 0) { - mpsin1(z, z, 0); - } else { - MPNumber t; - - /* cos(x) = sin(π/2 - |x|) */ - mp_get_pi(&t); - mp_divide_integer(&t, 2, &t); - mp_subtract(&t, z, z); - mp_sin(z, MP_RADIANS, z); - } + mpfr_const_pi(mpc_realref(z->num), MPFR_RNDN); + mpfr_set_zero(mpc_imagref(z->num), 0); } void mp_sin(const MPNumber *x, MPAngleUnit unit, MPNumber *z) { - if (mp_is_complex(x)) { - MPNumber x_real, x_im, z_real, z_im, t; - - mp_real_component(x, &x_real); - mp_imaginary_component(x, &x_im); - - mp_sin_real(&x_real, unit, &z_real); - mp_cosh(&x_im, &t); - mp_multiply(&z_real, &t, &z_real); - - mp_cos_real(&x_real, unit, &z_im); - mp_sinh(&x_im, &t); - mp_multiply(&z_im, &t, &z_im); - - mp_set_from_complex(&z_real, &z_im, z); - } + if (mp_is_complex(x)) + mp_set_from_mp(x, z); else - mp_sin_real(x, unit, z); + convert_to_radians(x, unit, z); + mpc_sin(z->num, z->num, MPC_RNDNN); } void mp_cos(const MPNumber *x, MPAngleUnit unit, MPNumber *z) { - if (mp_is_complex(x)) { - MPNumber x_real, x_im, z_real, z_im, t; - - mp_real_component(x, &x_real); - mp_imaginary_component(x, &x_im); - - mp_cos_real(&x_real, unit, &z_real); - mp_cosh(&x_im, &t); - mp_multiply(&z_real, &t, &z_real); - - mp_sin_real(&x_real, unit, &z_im); - mp_sinh(&x_im, &t); - mp_multiply(&z_im, &t, &z_im); - mp_invert_sign(&z_im, &z_im); - - mp_set_from_complex(&z_real, &z_im, z); - } + if (mp_is_complex(x)) + mp_set_from_mp(x, z); else - mp_cos_real(x, unit, z); + convert_to_radians(x, unit, z); + mpc_cos(z->num, z->num, MPC_RNDNN); } void mp_tan(const MPNumber *x, MPAngleUnit unit, MPNumber *z) { - MPNumber cos_x, sin_x; + MPNumber x_radians = mp_new(); + MPNumber pi = mp_new(); + MPNumber t1 = mp_new(); + + convert_to_radians(x, unit, &x_radians); + mp_get_pi(&pi); + mp_divide_integer(&pi, 2, &t1); + mp_subtract(&x_radians, &t1, &t1); + mp_divide(&t1, &pi, &t1); - /* Check for undefined values */ - mp_cos(x, unit, &cos_x); - if (mp_is_zero(&cos_x)) { + if (mp_is_integer(&t1)) { /* Translators: Error displayed when tangent value is undefined */ mperr(_("Tangent is undefined for angles that are multiples of π (180°) from π∕2 (90°)")); mp_set_from_integer(0, z); return; } - /* tan(x) = sin(x) / cos(x) */ - mp_sin(x, unit, &sin_x); - mp_divide(&sin_x, &cos_x, z); + if (mp_is_complex(x)) + mp_set_from_mp(x, z); + else + mp_set_from_mp(&x_radians, z); + mpc_tan(z->num, z->num, MPC_RNDNN); + mp_clear(&x_radians); + mp_clear(&pi); + mp_clear(&t1); } void mp_asin(const MPNumber *x, MPAngleUnit unit, MPNumber *z) { - MPNumber t1, t2; + MPNumber x_max = mp_new(); + MPNumber x_min = mp_new(); + mp_set_from_integer(1, &x_max); + mp_set_from_integer(-1, &x_min); - /* asin⁻¹(0) = 0 */ - if (mp_is_zero(x)) { + if (mp_compare(x, &x_max) > 0 || mp_compare(x, &x_min) < 0) + { + /* Translators: Error displayed when inverse sine value is undefined */ + mperr(_("Inverse sine is undefined for values outside [-1, 1]")); mp_set_from_integer(0, z); return; } - - /* sin⁻¹(x) = tan⁻¹(x / √(1 - x²)), |x| < 1 */ - if (x->exponent <= 0) { - mp_set_from_integer(1, &t1); - mp_set_from_mp(&t1, &t2); - mp_subtract(&t1, x, &t1); - mp_add(&t2, x, &t2); - mp_multiply(&t1, &t2, &t2); - mp_root(&t2, -2, &t2); - mp_multiply(x, &t2, z); - mp_atan(z, unit, z); - return; - } - - /* sin⁻¹(1) = π/2, sin⁻¹(-1) = -π/2 */ - mp_set_from_integer(x->sign, &t2); - if (mp_is_equal(x, &t2)) { - mp_get_pi(z); - mp_divide_integer(z, 2 * t2.sign, z); + mpc_asin(z->num, x->num, MPC_RNDNN); + if (!mp_is_complex(z)) convert_from_radians(z, unit, z); - return; - } - - /* Translators: Error displayed when inverse sine value is undefined */ - mperr(_("Inverse sine is undefined for values outside [-1, 1]")); - mp_set_from_integer(0, z); + mp_clear(&x_max); + mp_clear(&x_min); } void mp_acos(const MPNumber *x, MPAngleUnit unit, MPNumber *z) { - MPNumber t1, t2; - MPNumber MPn1, pi, MPy; - - mp_get_pi(&pi); - mp_set_from_integer(1, &t1); - mp_set_from_integer(-1, &MPn1); + MPNumber x_max = mp_new(); + MPNumber x_min = mp_new(); + mp_set_from_integer(1, &x_max); + mp_set_from_integer(-1, &x_min); - if (mp_is_greater_than(x, &t1) || mp_is_less_than(x, &MPn1)) { - /* Translators: Error displayed when inverse cosine value is undefined */ + if (mp_compare(x, &x_max) > 0 || mp_compare(x, &x_min) < 0) + { + /* Translators: Error displayed when inverse sine value is undefined */ mperr(_("Inverse cosine is undefined for values outside [-1, 1]")); mp_set_from_integer(0, z); - } else if (mp_is_zero(x)) { - mp_divide_integer(&pi, 2, z); - } else if (mp_is_equal(x, &t1)) { - mp_set_from_integer(0, z); - } else if (mp_is_equal(x, &MPn1)) { - mp_set_from_mp(&pi, z); - } else { - /* cos⁻¹(x) = tan⁻¹(√(1 - x²) / x) */ - mp_multiply(x, x, &t2); - mp_subtract(&t1, &t2, &t2); - mp_sqrt(&t2, &t2); - mp_divide(&t2, x, &t2); - mp_atan(&t2, MP_RADIANS, &MPy); - if (x->sign > 0) { - mp_set_from_mp(&MPy, z); - } else { - mp_add(&MPy, &pi, z); - } + return; } - - convert_from_radians(z, unit, z); + mpc_acos(z->num, x->num, MPC_RNDNN); + if (!mp_is_complex(z)) + convert_from_radians(z, unit, z); + mp_clear(&x_max); + mp_clear(&x_min); } void mp_atan(const MPNumber *x, MPAngleUnit unit, MPNumber *z) { - int i, q; - float rx = 0.0; - MPNumber t1, t2; + MPNumber i = mp_new(); + MPNumber minus_i = mp_new(); + mpc_set_si_si(i.num, 0, 1, MPC_RNDNN); + mpc_set_si_si(minus_i.num, 0, -1, MPC_RNDNN); - if (mp_is_zero(x)) { + /* Check x != i and x != -i */ + if (mp_is_equal(x, &i) || mp_is_equal(x, &minus_i)) + { + /* Translators: Error displayed when inverse sine value is undefined */ + mperr(_("Arctangent function is undefined for values i and -i")); mp_set_from_integer(0, z); return; } - - mp_set_from_mp(x, &t2); - if (abs(x->exponent) <= 2) - rx = mp_cast_to_float(x); - - /* REDUCE ARGUMENT IF NECESSARY BEFORE USING SERIES */ - q = 1; - while (t2.exponent >= 0) - { - if (t2.exponent == 0 && 2 * (t2.fraction[0] + 1) <= MP_BASE) - break; - - q *= 2; - - /* t = t / (√(t² + 1) + 1) */ - mp_multiply(&t2, &t2, z); - mp_add_integer(z, 1, z); - mp_sqrt(z, z); - mp_add_integer(z, 1, z); - mp_divide(&t2, z, &t2); - } - - /* USE POWER SERIES NOW ARGUMENT IN (-0.5, 0.5) */ - mp_set_from_mp(&t2, z); - mp_multiply(&t2, &t2, &t1); - - /* SERIES LOOP. REDUCE T IF POSSIBLE. */ - for (i = 1; ; i += 2) { - if (MP_T + 2 + t2.exponent <= 1) - break; - - mp_multiply(&t2, &t1, &t2); - mp_multiply_fraction(&t2, -i, i + 2, &t2); - - mp_add(z, &t2, z); - if (mp_is_zero(&t2)) - break; - } - - /* CORRECT FOR ARGUMENT REDUCTION */ - mp_multiply_integer(z, q, z); - - /* CHECK THAT RELATIVE ERROR LESS THAN 0.01 UNLESS EXPONENT - * OF X IS LARGE (WHEN ATAN MIGHT NOT WORK) - */ - if (abs(x->exponent) <= 2) { - float ry = mp_cast_to_float(z); - /* THE FOLLOWING MESSAGE MAY INDICATE THAT B**(T-1) IS TOO SMALL. */ - if (fabs(ry - atan(rx)) >= fabs(ry) * 0.01) - mperr("*** ERROR OCCURRED IN MP_ATAN, RESULT INCORRECT ***"); - } - - convert_from_radians(z, unit, z); + mpc_atan(z->num, x->num, MPC_RNDNN); + if (!mp_is_complex(z)) + convert_from_radians(z, unit, z); + mp_clear(&i); + mp_clear(&minus_i); } void mp_sinh(const MPNumber *x, MPNumber *z) { - MPNumber abs_x; - - /* sinh(0) = 0 */ - if (mp_is_zero(x)) { - mp_set_from_integer(0, z); - return; - } - - /* WORK WITH ABS(X) */ - mp_abs(x, &abs_x); - - /* If |x| < 1 USE MPEXP TO AVOID CANCELLATION, otherwise IF TOO LARGE MP_EPOWY GIVES ERROR MESSAGE */ - if (abs_x.exponent <= 0) { - MPNumber exp_x, a, b; - - /* ((e^|x| + 1) * (e^|x| - 1)) / e^|x| */ - // FIXME: Solves to e^|x| - e^-|x|, why not lower branch always? */ - mp_epowy(&abs_x, &exp_x); - mp_add_integer(&exp_x, 1, &a); - mp_add_integer(&exp_x, -1, &b); - mp_multiply(&a, &b, z); - mp_divide(z, &exp_x, z); - } - else { - MPNumber exp_x; - - /* e^|x| - e^-|x| */ - mp_epowy(&abs_x, &exp_x); - mp_reciprocal(&exp_x, z); - mp_subtract(&exp_x, z, z); - } - - /* DIVIDE BY TWO AND RESTORE SIGN */ - mp_divide_integer(z, 2, z); - mp_multiply_integer(z, x->sign, z); + mpc_sinh(z->num, x->num, MPC_RNDNN); } void mp_cosh(const MPNumber *x, MPNumber *z) { - MPNumber t; - - /* cosh(0) = 1 */ - if (mp_is_zero(x)) { - mp_set_from_integer(1, z); - return; - } - - /* cosh(x) = (e^x + e^-x) / 2 */ - mp_abs(x, &t); - mp_epowy(&t, &t); - mp_reciprocal(&t, z); - mp_add(&t, z, z); - mp_divide_integer(z, 2, z); + mpc_cosh(z->num, x->num, MPC_RNDNN); } void mp_tanh(const MPNumber *x, MPNumber *z) { - float r__1; - MPNumber t; - - /* tanh(0) = 0 */ - if (mp_is_zero(x)) { - mp_set_from_integer(0, z); - return; - } - - mp_abs(x, &t); - - /* SEE IF ABS(X) SO LARGE THAT RESULT IS +-1 */ - r__1 = (float) MP_T * 0.5 * log((float) MP_BASE); - mp_set_from_float(r__1, z); - if (mp_compare_mp_to_mp(&t, z) > 0) { - mp_set_from_integer(x->sign, z); - return; - } - - /* If |x| >= 1/2 use ?, otherwise use ? to avoid cancellation */ - /* |tanh(x)| = (e^|2x| - 1) / (e^|2x| + 1) */ - mp_multiply_integer(&t, 2, &t); - if (t.exponent > 0) { - mp_epowy(&t, &t); - mp_add_integer(&t, -1, z); - mp_add_integer(&t, 1, &t); - mp_divide(z, &t, z); - } else { - mp_epowy(&t, &t); - mp_add_integer(&t, 1, z); - mp_add_integer(&t, -1, &t); - mp_divide(&t, z, z); - } - - /* Restore sign */ - z->sign = x->sign * z->sign; + mpc_tanh(z->num, x->num, MPC_RNDNN); } void mp_asinh(const MPNumber *x, MPNumber *z) { - MPNumber t; - - /* sinh⁻¹(x) = ln(x + √(x² + 1)) */ - mp_multiply(x, x, &t); - mp_add_integer(&t, 1, &t); - mp_sqrt(&t, &t); - mp_add(x, &t, &t); - mp_ln(&t, z); + mpc_asinh(z->num, x->num, MPC_RNDNN); } void mp_acosh(const MPNumber *x, MPNumber *z) { - MPNumber t; + MPNumber t = mp_new(); /* Check x >= 1 */ mp_set_from_integer(1, &t); - if (mp_is_less_than(x, &t)) { + if (mp_is_less_than(x, &t)) + { /* Translators: Error displayed when inverse hyperbolic cosine value is undefined */ mperr(_("Inverse hyperbolic cosine is undefined for values less than one")); mp_set_from_integer(0, z); return; } - /* cosh⁻¹(x) = ln(x + √(x² - 1)) */ - mp_multiply(x, x, &t); - mp_add_integer(&t, -1, &t); - mp_sqrt(&t, &t); - mp_add(x, &t, &t); - mp_ln(&t, z); + mpc_acosh(z->num, x->num, MPC_RNDNN); + mp_clear(&t); } void mp_atanh(const MPNumber *x, MPNumber *z) { - MPNumber one, minus_one, n, d; - - /* Check -1 <= x <= 1 */ - mp_set_from_integer(1, &one); - mp_set_from_integer(-1, &minus_one); - if (mp_is_greater_equal(x, &one) || mp_is_less_equal(x, &minus_one)) { - /* Translators: Error displayed when inverse hyperbolic tangent value is undefined */ - mperr(_("Inverse hyperbolic tangent is undefined for values outside [-1, 1]")); + MPNumber x_max = mp_new(); + MPNumber x_min = mp_new(); + mp_set_from_integer(1, &x_max); + mp_set_from_integer(-1, &x_min); + + if (mp_compare(x, &x_max) >= 0 || mp_compare(x, &x_min) <= 0) + { + /* Translators: Error displayed when inverse hyperbolic tangent value is undefined */ + mperr(_("Inverse hyperbolic tangent is undefined for values outside (-1, 1)")); mp_set_from_integer(0, z); return; } - - /* atanh(x) = 0.5 * ln((1 + x) / (1 - x)) */ - mp_add_integer(x, 1, &n); - mp_set_from_mp(x, &d); - mp_invert_sign(&d, &d); - mp_add_integer(&d, 1, &d); - mp_divide(&n, &d, z); - mp_ln(z, z); - mp_divide_integer(z, 2, z); + mpc_atanh(z->num, x->num, MPC_RNDNN); + mp_clear(&x_max); + mp_clear(&x_min); } -- cgit v1.2.1