- <range min="0" max="15"/>

- mp-private.h \

- -lm

- -lm

- -lm

- MPNumber value;

- MPNumber r, v;

- MPNumber r, v;

- MPNumber MP1, MP2, MP3, MP4;

- MPNumber MPbv, MP1, MP2;

- len = mp_cast_to_int(period);

- MPNumber MP1, MP2, MP3, MP4;

- MPNumber MP1, MP2;

- MPNumber MP1, MP2, MP3, MP4;

- MPNumber MP1, MP2, MP3, MP4;

- MPNumber MP1, MP2, MP3, MP4;

- MPNumber MP1;

- MPNumber MP1, MP2, MP3, MP4;

- MPNumber MP1, MP2, MP3, MP4;

- MPNumber result;

- MPNumber tmp;

- MPNumber z;

- MPNumber result;

- MPNumber x;

- MPNumber max, fraction;

- bits = mp_cast_to_unsigned_int(&x);

- MPNumber x;

-

-

- MPNumber x, z;

- MPNumber x, z;

- MPNumber t;

-

- MPNumber z;

- MPNumber x;

- MPNumber x;

-

- MPNumber z;

- MPNumber x;

- MPNumber max;

- bits = mp_cast_to_unsigned_int(&x);

- MPNumber z;

- MPNumber z;

-#include "mp-private.h"

- MPNumber tmp1, tmp2;

- return mp_is_greater_than (&tmp2, x);

- MPNumber temp;

- mp_set_from_integer(0, &temp);

- int i;

- MPNumber multiplier;

-

- mp_set_from_integer(1, &multiplier);

-

- for (i = 0; i < count; i++)

- for (i = 0; i < -count; i++)

- MPNumber t;

-#include "mp-private.h"

- if (z != x)

- memcpy(z, x, sizeof(MPNumber));

-}

-

-

-void

-mp_set_from_float(float rx, MPNumber *z)

-{

- int i, k, ib, ie, tp;

- float rj;

-

- mp_set_from_integer(0, z);

-

- /* CHECK SIGN */

- if (rx < 0.0f) {

- z->sign = -1;

- rj = -(double)(rx);

- } else if (rx > 0.0f) {

- z->sign = 1;

- rj = rx;

- } else {

- /* IF RX = 0E0 RETURN 0 */

- mp_set_from_integer(0, z);

- return;

- }

-

- /* INCREASE IE AND DIVIDE RJ BY 16. */

- ie = 0;

- while (rj >= 1.0f) {

- ++ie;

- rj *= 0.0625f;

- }

- while (rj < 0.0625f) {

- --ie;

- rj *= 16.0f;

- }

-

- /* NOW RJ IS DY DIVIDED BY SUITABLE POWER OF 16.

- * SET EXPONENT TO 0

- */

- z->exponent = 0;

-

- /* CONVERSION LOOP (ASSUME SINGLE-PRECISION OPS. EXACT) */

- for (i = 0; i < MP_T + 4; i++) {

- rj *= (float) MP_BASE;

- z->fraction[i] = (int) rj;

- rj -= (float) z->fraction[i];

- }

-

- /* NORMALIZE RESULT */

- mp_normalize(z);

-

- /* Computing MAX */

- ib = max(MP_BASE * 7 * MP_BASE, 32767) / 16;

- tp = 1;

-

- /* NOW MULTIPLY BY 16**IE */

- if (ie < 0) {

- k = -ie;

- for (i = 1; i <= k; ++i) {

- tp <<= 4;

- if (tp <= ib && tp != MP_BASE && i < k)

- continue;

- mp_divide_integer(z, tp, z);

- tp = 1;

- }

- } else if (ie > 0) {

- for (i = 1; i <= ie; ++i) {

- tp <<= 4;

- if (tp <= ib && tp != MP_BASE && i < ie)

- continue;

- mp_multiply_integer(z, tp, z);

- tp = 1;

- }

- }

- int i, k, ib, ie, tp;

- double dj;

-

- mp_set_from_integer(0, z);

-

- /* CHECK SIGN */

- if (dx < 0.0) {

- z->sign = -1;

- dj = -dx;

- } else if (dx > 0.0) {

- z->sign = 1;

- dj = dx;

- } else {

- mp_set_from_integer(0, z);

- return;

- }

-

- /* INCREASE IE AND DIVIDE DJ BY 16. */

- for (ie = 0; dj >= 1.0; ie++)

- dj *= 1.0/16.0;

-

- for ( ; dj < 1.0/16.0; ie--)

- dj *= 16.0;

-

- /* NOW DJ IS DY DIVIDED BY SUITABLE POWER OF 16

- * SET EXPONENT TO 0

- */

- z->exponent = 0;

-

- /* CONVERSION LOOP (ASSUME DOUBLE-PRECISION OPS. EXACT) */

- for (i = 0; i < MP_T + 4; i++) {

- dj *= (double) MP_BASE;

- z->fraction[i] = (int) dj;

- dj -= (double) z->fraction[i];

- }

-

- /* NORMALIZE RESULT */

- mp_normalize(z);

-

- /* Computing MAX */

- ib = max(MP_BASE * 7 * MP_BASE, 32767) / 16;

- tp = 1;

-

- /* NOW MULTIPLY BY 16**IE */

- if (ie < 0) {

- k = -ie;

- for (i = 1; i <= k; ++i) {

- tp <<= 4;

- if (tp <= ib && tp != MP_BASE && i < k)

- continue;

- mp_divide_integer(z, tp, z);

- tp = 1;

- }

- } else if (ie > 0) {

- for (i = 1; i <= ie; ++i) {

- tp <<= 4;

- if (tp <= ib && tp != MP_BASE && i < ie)

- continue;

- mp_multiply_integer(z, tp, z);

- tp = 1;

- }

- }

- int i;

-

- memset(z, 0, sizeof(MPNumber));

-

- if (x == 0) {

- z->sign = 0;

- return;

- }

-

- if (x < 0) {

- x = -x;

- z->sign = -1;

- }

- else if (x > 0)

- z->sign = 1;

-

- while (x != 0) {

- z->fraction[z->exponent] = x % MP_BASE;

- z->exponent++;

- x /= MP_BASE;

- }

- for (i = 0; i < z->exponent / 2; i++) {

- int t = z->fraction[i];

- z->fraction[i] = z->fraction[z->exponent - i - 1];

- z->fraction[z->exponent - i - 1] = t;

- }

- int i;

-

- mp_set_from_integer(0, z);

-

- if (x == 0) {

- z->sign = 0;

- return;

- }

- z->sign = 1;

-

- while (x != 0) {

- z->fraction[z->exponent] = x % MP_BASE;

- x = x / MP_BASE;

- z->exponent++;

- }

- for (i = 0; i < z->exponent / 2; i++) {

- int t = z->fraction[i];

- z->fraction[i] = z->fraction[z->exponent - i - 1];

- z->fraction[z->exponent - i - 1] = t;

- }

- mp_gcd(&numerator, &denominator);

-

- if (denominator == 0) {

- mperr("*** J == 0 IN CALL TO MP_SET_FROM_FRACTION ***\n");

- mp_set_from_integer(0, z);

- return;

- }

-

- MPNumber x, y;

-}

- /* NOTE: Do imaginary component first as z may be x or y */

- z->im_sign = y->sign;

- z->im_exponent = y->exponent;

- memcpy(z->im_fraction, y->fraction, sizeof(int) * MP_SIZE);

-

- z->sign = x->sign;

- z->exponent = x->exponent;

- if (z != x)

- memcpy(z->fraction, x->fraction, sizeof(int) * MP_SIZE);

-

-

-mp_cast_to_int(const MPNumber *x)

- int i;

- int64_t z = 0, v;

-

- /* |x| <= 1 */

- if (x->sign == 0 || x->exponent <= 0)

- return 0;

-

- /* Multiply digits together */

- for (i = 0; i < x->exponent; i++) {

- int64_t t;

-

- t = z;

- z = z * MP_BASE + x->fraction[i];

-

- /* Check for overflow */

- if (z <= t)

- return 0;

- }

-

- /* Validate result */

- v = z;

- for (i = x->exponent - 1; i >= 0; i--) {

- int64_t digit;

-

- /* Get last digit */

- digit = v - (v / MP_BASE) * MP_BASE;

- if (x->fraction[i] != digit)

- return 0;

-

- v /= MP_BASE;

- }

- if (v != 0)

- return 0;

-

- return x->sign * z;

-mp_cast_to_unsigned_int(const MPNumber *x)

-{

- int i;

- uint64_t z = 0, v;

-

- /* x <= 1 */

- if (x->sign <= 0 || x->exponent <= 0)

- return 0;

-

- /* Multiply digits together */

- for (i = 0; i < x->exponent; i++) {

- uint64_t t;

-

- t = z;

- z = z * MP_BASE + x->fraction[i];

-

- /* Check for overflow */

- if (z <= t)

- return 0;

- }

-

- /* Validate result */

- v = z;

- for (i = x->exponent - 1; i >= 0; i--) {

- uint64_t digit;

-

- /* Get last digit */

- digit = v - (v / MP_BASE) * MP_BASE;

- if (x->fraction[i] != digit)

- return 0;

-

- v /= MP_BASE;

- }

- if (v != 0)

- return 0;

-

- return z;

-}

-

-

-static double

-mppow_ri(float ap, int bp)

- double pow;

-

- if (bp == 0)

- return 1.0;

-

- if (bp < 0) {

- if (ap == 0)

- return 1.0;

- bp = -bp;

- ap = 1 / ap;

- }

-

- pow = 1.0;

- for (;;) {

- if (bp & 01)

- pow *= ap;

- if (bp >>= 1)

- ap *= ap;

- else

- break;

- }

-

- return pow;

-

-mp_cast_to_float(const MPNumber *x)

-{

- int i;

- float rz = 0.0;

-

- if (mp_is_zero(x))

- return 0.0;

-

- for (i = 0; i < MP_T; i++) {

- rz = (float) MP_BASE * rz + (float)x->fraction[i];

-

- /* CHECK IF FULL SINGLE-PRECISION ACCURACY ATTAINED */

- if (rz + 1.0f <= rz)

- break;

- }

-

- /* NOW ALLOW FOR EXPONENT */

- rz *= mppow_ri((float) MP_BASE, x->exponent - i - 1);

-

- /* CHECK REASONABLENESS OF RESULT */

- /* LHS SHOULD BE <= 0.5, BUT ALLOW FOR SOME ERROR IN ALOG */

- if (rz <= (float)0. ||

- fabs((float) x->exponent - (log(rz) / log((float) MP_BASE) + (float).5)) > (float).6) {

- /* FOLLOWING MESSAGE INDICATES THAT X IS TOO LARGE OR SMALL -

- * TRY USING MPCMRE INSTEAD.

- */

- mperr("*** FLOATING-POINT OVER/UNDER-FLOW IN MP_CAST_TO_FLOAT ***\n");

- return 0.0;

- }

-

- if (x->sign < 0)

- rz = -(double)(rz);

-

- return rz;

-}

-

-

-static double

-mppow_di(double ap, int bp)

- double pow = 1.0;

-

- if (bp != 0) {

- if (bp < 0) {

- if (ap == 0) return(pow);

- bp = -bp;

- ap = 1/ap;

- }

- for (;;) {

- if (bp & 01) pow *= ap;

- if (bp >>= 1) ap *= ap;

- else break;

- }

- }

-

- return(pow);

-

-mp_cast_to_double(const MPNumber *x)

- int i, tm = 0;

- double d__1, dz2, ret_val = 0.0;

-

- if (mp_is_zero(x))

- return 0.0;

-

- for (i = 0; i < MP_T; i++) {

- ret_val = (double) MP_BASE * ret_val + (double) x->fraction[i];

- tm = i;

-

- /* CHECK IF FULL DOUBLE-PRECISION ACCURACY ATTAINED */

- dz2 = ret_val + 1.0;

-

- /* TEST BELOW NOT ALWAYS EQUIVALENT TO - IF (DZ2.LE.DZ) GO TO 20,

- * FOR EXAMPLE ON CYBER 76.

- */

- if (dz2 - ret_val <= 0.0)

- break;

- }

-

- /* NOW ALLOW FOR EXPONENT */

- ret_val *= mppow_di((double) MP_BASE, x->exponent - tm - 1);

-

- /* CHECK REASONABLENESS OF RESULT. */

- /* LHS SHOULD BE .LE. 0.5 BUT ALLOW FOR SOME ERROR IN DLOG */

- if (ret_val <= 0. ||

- ((d__1 = (double) ((float) x->exponent) - (log(ret_val) / log((double)

- ((float) MP_BASE)) + .5), abs(d__1)) > .6)) {

- /* FOLLOWING MESSAGE INDICATES THAT X IS TOO LARGE OR SMALL -

- * TRY USING MPCMDE INSTEAD.

- */

- mperr("*** FLOATING-POINT OVER/UNDER-FLOW IN MP_CAST_TO_DOUBLE ***\n");

- return 0.0;

- }

- else

- {

- if (x->sign < 0)

- ret_val = -ret_val;

- return ret_val;

- }

- MPNumber t;

- MPNumber fraction;

- MPNumber numerator, denominator;

- MPNumber t;

-/*

- * Copyright (C) 1987-2008 Sun Microsystems, Inc. All Rights Reserved.

- * Copyright (C) 2008-2009 Robert Ancell

- *

- * This program is free software: you can redistribute it and/or modify it under

- * the terms of the GNU General Public License as published by the Free Software

- * Foundation, either version 2 of the License, or (at your option) any later

- * version. See http://www.gnu.org/copyleft/gpl.html the full text of the

- * license.

- */

-

-#ifndef MP_INTERNAL_H

-#define MP_INTERNAL_H

-

-#include <glib/gi18n.h>

-

-/* If we're not using GNU C, elide __attribute__ */

-#ifndef __GNUC__

-# define __attribute__(x) /*NOTHING*/

-#endif

-

-#define min(a, b) ((a) <= (b) ? (a) : (b))

-#define max(a, b) ((a) >= (b) ? (a) : (b))

-

-//2E0 BELOW ENSURES AT LEAST ONE GUARD DIGIT

-//MP.t = (int) ((float) (accuracy) * log((float)10.) / log((float) MP_BASE) + (float) 2.0);

-//if (MP.t > MP_SIZE) {

-// mperr("MP_SIZE TOO SMALL IN CALL TO MPSET, INCREASE MP_SIZE AND DIMENSIONS OF MP ARRAYS TO AT LEAST %d ***", MP.t);

-// MP.t = MP_SIZE;

-//}

-#define MP_T 100

-

-void mperr(const char *format, ...) __attribute__((format(printf, 1, 2)));

-void mp_gcd(int64_t *, int64_t *);

-void mp_normalize(MPNumber *);

-void convert_to_radians(const MPNumber *x, MPAngleUnit unit, MPNumber *z);

-void convert_from_radians(const MPNumber *x, MPAngleUnit unit, MPNumber *z);

-

-#endif /* MP_INTERNAL_H */

-mp_cast_to_string_real(MpSerializer *serializer, const MPNumber *x, int base, gboolean force_sign, int *n_digits, GString *string)

- MPNumber number, integer_component, fractional_component, temp;

- MPNumber t, t2, t3;

- d = mp_cast_to_int(&t3);

- g_string_prepend_c(string, d < 16 ? digits[d] : '?');

- MPNumber digit;

- d = mp_cast_to_int(&digit);

-mp_cast_to_string(MpSerializer *serializer, const MPNumber *x, int *n_digits)

- MPNumber x_real;

- mp_cast_to_string_real(serializer, &x_real, serializer->priv->base, FALSE, n_digits, string);

- MPNumber x_im;

- mp_cast_to_string_real(serializer, &x_im, 10, force_sign, &n_complex_digits, s);

-mp_cast_to_exponential_string_real(MpSerializer *serializer, const MPNumber *x, GString *string, gboolean eng_format, int *n_digits)

- MPNumber t, z, base, base3, base10, base10inv, mantissa;

- fixed = mp_cast_to_string(serializer, &mantissa, n_digits);

-mp_cast_to_exponential_string(MpSerializer *serializer, const MPNumber *x, gboolean eng_format, int *n_digits)

- MPNumber x_real;

- exponent = mp_cast_to_exponential_string_real(serializer, &x_real, string, eng_format, n_digits);

- MPNumber x_im;

- exponent = mp_cast_to_exponential_string_real(serializer, &x_im, s, eng_format, &n_complex_digits);

- MPNumber cmp, xcmp;

- mp_set_from_integer(10, &cmp);

- mp_xpowy_integer(&cmp, -(serializer->priv->trailing_digits), &cmp);

- mp_real_component(x, &xcmp);

- mp_abs(&xcmp, &xcmp);

- s0 = mp_cast_to_string(serializer, x, &n_digits);

- if (n_digits <= serializer->priv->leading_digits &&

- mp_is_greater_equal(&xcmp, &cmp))

- return mp_cast_to_exponential_string(serializer, x, FALSE, &n_digits);

- return mp_cast_to_string(serializer, x, &n_digits);

- return mp_cast_to_exponential_string(serializer, x, FALSE, &n_digits);

- return mp_cast_to_exponential_string(serializer, x, TRUE, &n_digits);

-#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);

-}

-

- MPNumber t1, t2;

- break;

- mp_get_pi(&t1);

- mp_multiply(x, &t1, &t2);

- mp_divide_integer(&t2, 180, z);

- mp_get_pi(&t1);

- mp_multiply(x, &t1, &t2);

- mp_divide_integer(&t2, 200, z);

-

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

-}

-

-

- MPNumber t1, t2;

- switch (unit) {

- break;

- mp_multiply_integer(x, 180, &t2);

- mp_get_pi(&t1);

- mp_divide(&t2, &t1, z);

- mp_multiply_integer(x, 200, &t2);

- mp_get_pi(&t1);

- mp_divide(&t2, &t1, z);

-/* 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)

- /* 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);

- }

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

- }

- mp_sin_real(x, unit, 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);

- }

- mp_cos_real(x, unit, z);

- MPNumber cos_x, sin_x;

- /* Check for undefined values */

- mp_cos(x, unit, &cos_x);

- if (mp_is_zero(&cos_x)) {

- /* tan(x) = sin(x) / cos(x) */

- mp_sin(x, unit, &sin_x);

- mp_divide(&sin_x, &cos_x, z);

- MPNumber t1, t2;

- /* asin⁻¹(0) = 0 */

- if (mp_is_zero(x)) {

-

- /* 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);

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

- MPNumber t1, t2;

- MPNumber MPn1, pi, MPy;

-

- mp_get_pi(&pi);

- mp_set_from_integer(1, &t1);

- mp_set_from_integer(-1, &MPn1);

- if (mp_is_greater_than(x, &t1) || mp_is_less_than(x, &MPn1)) {

- /* Translators: Error displayed when inverse cosine value is undefined */

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

- }

-

- convert_from_radians(z, unit, z);

- int i, q;

- float rx = 0.0;

- MPNumber t1, t2;

- if (mp_is_zero(x)) {

-

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

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

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

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

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

- MPNumber t;

- if (mp_is_less_than(x, &t)) {

- /* 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);

- 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]"));

-

- /* 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);

-#include "mp-private.h"

-

-static MPNumber eulers_number;

-static gboolean have_eulers_number = FALSE;

-

-// FIXME: Re-add overflow and underflow detection

-/* ROUTINE CALLED BY MP_DIVIDE AND MP_SQRT TO ENSURE THAT

- * RESULTS ARE REPRESENTED EXACTLY IN T-2 DIGITS IF THEY

- * CAN BE. X IS AN MP NUMBER, I AND J ARE INTEGERS.

- */

-static void

-mp_ext(int i, int j, MPNumber *x)

- int q, s;

- if (mp_is_zero(x) || MP_T <= 2 || i == 0)

- return;

- /* COMPUTE MAXIMUM POSSIBLE ERROR IN THE LAST PLACE */

- q = (j + 1) / i + 1;

- s = MP_BASE * x->fraction[MP_T - 2] + x->fraction[MP_T - 1];

- /* SET LAST TWO DIGITS TO ZERO */

- if (s <= q) {

- x->fraction[MP_T - 2] = 0;

- x->fraction[MP_T - 1] = 0;

- return;

- }

- if (s + q < MP_BASE * MP_BASE)

- return;

- /* ROUND UP HERE */

- x->fraction[MP_T - 2] = MP_BASE - 1;

- x->fraction[MP_T - 1] = MP_BASE;

- /* NORMALIZE X (LAST DIGIT B IS OK IN MP_MULTIPLY_INTEGER) */

- mp_multiply_integer(x, 1, x);

- if (!have_eulers_number) {

- MPNumber t;

- mp_set_from_integer(1, &t);

- mp_epowy(&t, &eulers_number);

- have_eulers_number = TRUE;

- }

- mp_set_from_mp(&eulers_number, z);

- mp_set_from_integer(0, z);

- z->im_sign = 1;

- z->im_exponent = 1;

- z->im_fraction[0] = 1;

- if (mp_is_complex(x)){

- MPNumber x_real, x_im;

-

- mp_real_component(x, &x_real);

- mp_imaginary_component(x, &x_im);

-

- mp_multiply(&x_real, &x_real, &x_real);

- mp_multiply(&x_im, &x_im, &x_im);

- mp_add(&x_real, &x_im, z);

- mp_root(z, 2, z);

- }

- else {

- mp_set_from_mp(x, z);

- if (z->sign < 0)

- z->sign = -z->sign;

- }

- MPNumber x_real, x_im, pi;

-

- if (mp_is_zero(x)) {

- mp_real_component(x, &x_real);

- mp_imaginary_component(x, &x_im);

- mp_get_pi(&pi);

-

- if (mp_is_zero(&x_im)) {

- if (mp_is_negative(&x_real))

- convert_from_radians(&pi, MP_RADIANS, z);

- else

- mp_set_from_integer(0, z);

- }

- else if (mp_is_zero(&x_real)) {

- mp_set_from_mp(&pi, z);

- if (mp_is_negative(&x_im))

- mp_divide_integer(z, -2, z);

- else

- mp_divide_integer(z, 2, z);

- }

- else if (mp_is_negative(&x_real)) {

- mp_divide(&x_im, &x_real, z);

- mp_atan(z, MP_RADIANS, z);

- if (mp_is_negative(&x_im))

- mp_subtract(z, &pi, z);

- else

- mp_add(z, &pi, z);

- }

- else {

- mp_divide(&x_im, &x_real, z);

- mp_atan(z, MP_RADIANS, z);

- }

-

- mp_set_from_mp(x, z);

- z->im_sign = -z->im_sign;

- mp_set_from_mp(x, z);

-

- /* Clear imaginary component */

- z->im_sign = 0;

- z->im_exponent = 0;

- memset(z->im_fraction, 0, sizeof(int) * MP_SIZE);

- /* Copy imaginary component to real component */

- z->sign = x->im_sign;

- z->exponent = x->im_exponent;

- memcpy(z->fraction, x->im_fraction, sizeof(int) * MP_SIZE);

-

- /* Clear (old) imaginary component */

- z->im_sign = 0;

- z->im_exponent = 0;

- memset(z->im_fraction, 0, sizeof(int) * MP_SIZE);

-

-static void

-mp_add_real(const MPNumber *x, int y_sign, const MPNumber *y, MPNumber *z)

-{

- int sign_prod, i, c;

- int exp_diff, med;

- bool x_largest = false;

- const int *big_fraction, *small_fraction;

- MPNumber x_copy, y_copy;

-

- /* 0 + y = y */

- if (mp_is_zero(x)) {

- mp_set_from_mp(y, z);

- z->sign = y_sign;

- return;

- }

- /* x + 0 = x */

- else if (mp_is_zero(y)) {

- mp_set_from_mp(x, z);

- return;

- }

-

- sign_prod = y_sign * x->sign;

- exp_diff = x->exponent - y->exponent;

- med = abs(exp_diff);

- if (exp_diff < 0) {

- x_largest = false;

- } else if (exp_diff > 0) {

- x_largest = true;

- } else {

- /* EXPONENTS EQUAL SO COMPARE SIGNS, THEN FRACTIONS IF NEC. */

- if (sign_prod < 0) {

- /* Signs are not equal. find out which mantissa is larger. */

- int j;

- for (j = 0; j < MP_T; j++) {

- int i = x->fraction[j] - y->fraction[j];

- if (i == 0)

- continue;

- if (i < 0)

- x_largest = false;

- else if (i > 0)

- x_largest = true;

- break;

- }

-

- /* Both mantissas equal, so result is zero. */

- if (j >= MP_T) {

- mp_set_from_integer(0, z);

- return;

- }

- }

- }

-

- mp_set_from_mp(x, &x_copy);

- mp_set_from_mp(y, &y_copy);

- mp_set_from_integer(0, z);

-

- if (x_largest) {

- z->sign = x_copy.sign;

- z->exponent = x_copy.exponent;

- big_fraction = x_copy.fraction;

- small_fraction = y_copy.fraction;

- } else {

- z->sign = y_sign;

- z->exponent = y_copy.exponent;

- big_fraction = y_copy.fraction;

- small_fraction = x_copy.fraction;

- }

-

- /* CLEAR GUARD DIGITS TO RIGHT OF X DIGITS */

- for(i = 3; i >= med; i--)

- z->fraction[MP_T + i] = 0;

-

- if (sign_prod >= 0) {

- /* This is probably insufficient overflow detection, but it makes us

- * not crash at least.

- */

- if (MP_T + 3 < med) {

- mperr(_("Overflow: the result couldn't be calculated"));

- mp_set_from_integer(0, z);

- return;

- }

-

- /* HERE DO ADDITION, EXPONENT(Y) >= EXPONENT(X) */

- for (i = MP_T + 3; i >= MP_T; i--)

- z->fraction[i] = small_fraction[i - med];

-

- c = 0;

- for (; i >= med; i--) {

- c = big_fraction[i] + small_fraction[i - med] + c;

-

- if (c < MP_BASE) {

- /* NO CARRY GENERATED HERE */

- z->fraction[i] = c;

- c = 0;

- } else {

- /* CARRY GENERATED HERE */

- z->fraction[i] = c - MP_BASE;

- c = 1;

- }

- }

-

- for (; i >= 0; i--)

- {

- c = big_fraction[i] + c;

- if (c < MP_BASE) {

- z->fraction[i] = c;

- i--;

-

- /* NO CARRY POSSIBLE HERE */

- for (; i >= 0; i--)

- z->fraction[i] = big_fraction[i];

-

- c = 0;

- break;

- }

-

- z->fraction[i] = 0;

- c = 1;

- }

-

- /* MUST SHIFT RIGHT HERE AS CARRY OFF END */

- if (c != 0) {

- for (i = MP_T + 3; i > 0; i--)

- z->fraction[i] = z->fraction[i - 1];

- z->fraction[0] = 1;

- z->exponent++;

- }

- }

- else {

- c = 0;

- for (i = MP_T + med - 1; i >= MP_T; i--) {

- /* HERE DO SUBTRACTION, ABS(Y) > ABS(X) */

- z->fraction[i] = c - small_fraction[i - med];

- c = 0;

-

- /* BORROW GENERATED HERE */

- if (z->fraction[i] < 0) {

- c = -1;

- z->fraction[i] += MP_BASE;

- }

- }

-

- for(; i >= med; i--) {

- c = big_fraction[i] + c - small_fraction[i - med];

- if (c >= 0) {

- /* NO BORROW GENERATED HERE */

- z->fraction[i] = c;

- c = 0;

- } else {

- /* BORROW GENERATED HERE */

- z->fraction[i] = c + MP_BASE;

- c = -1;

- }

- }

-

- for (; i >= 0; i--) {

- c = big_fraction[i] + c;

-

- if (c >= 0) {

- z->fraction[i] = c;

- i--;

-

- /* NO CARRY POSSIBLE HERE */

- for (; i >= 0; i--)

- z->fraction[i] = big_fraction[i];

-

- break;

- }

-

- z->fraction[i] = c + MP_BASE;

- c = -1;

- }

- }

-

- mp_normalize(z);

-}

-

-

-static void

-mp_add_with_sign(const MPNumber *x, int y_sign, const MPNumber *y, MPNumber *z)

-{

- if (mp_is_complex(x) || mp_is_complex(y)) {

- MPNumber real_x, real_y, im_x, im_y, real_z, im_z;

-

- mp_real_component(x, &real_x);

- mp_imaginary_component(x, &im_x);

- mp_real_component(y, &real_y);

- mp_imaginary_component(y, &im_y);

-

- mp_add_real(&real_x, y_sign * y->sign, &real_y, &real_z);

- mp_add_real(&im_x, y_sign * y->im_sign, &im_y, &im_z);

-

- mp_set_from_complex(&real_z, &im_z, z);

- }

- else

- mp_add_real(x, y_sign * y->sign, y, z);

-}

-

-

- mp_add_with_sign(x, 1, y, z);

- MPNumber t;

- mp_set_from_integer(y, &t);

- mp_add(x, &t, z);

-

-void

-mp_add_fraction(const MPNumber *x, int64_t i, int64_t j, MPNumber *y)

-{

- MPNumber t;

- mp_set_from_fraction(i, j, &t);

- mp_add(x, &t, y);

-}

-

-

- mp_add_with_sign(x, -1, y, z);

-

- if (mp_is_zero(x))

- mp_set_from_integer(0, z);

- else if (mp_is_negative(x))

- mp_set_from_integer(-1, z);

- else

- mp_set_from_integer(1, z);

- int i;

-

- /* Clear fraction */

- mp_set_from_mp(x, z);

- for (i = z->exponent; i < MP_SIZE; i++)

- z->fraction[i] = 0;

- z->im_sign = 0;

- z->im_exponent = 0;

- memset(z->im_fraction, 0, sizeof(int) * MP_SIZE);

- int i, shift;

-

- /* Fractional component of zero is 0 */

- if (mp_is_zero(x)) {

- mp_set_from_integer(0, z);

- return;

- }

-

- /* All fractional */

- if (x->exponent <= 0) {

- mp_set_from_mp(x, z);

- return;

- }

-

- /* Shift fractional component */

- shift = x->exponent;

- for (i = shift; i < MP_SIZE && x->fraction[i] == 0; i++)

- shift++;

- z->sign = x->sign;

- z->exponent = x->exponent - shift;

- for (i = 0; i < MP_SIZE; i++) {

- if (i + shift >= MP_SIZE)

- z->fraction[i] = 0;

- else

- z->fraction[i] = x->fraction[i + shift];

- }

- if (z->fraction[0] == 0)

- z->sign = 0;

-

- z->im_sign = 0;

- z->im_exponent = 0;

- memset(z->im_fraction, 0, sizeof(int) * MP_SIZE);

- MPNumber f;

- int i;

- bool have_fraction = false, is_negative;

-

- /* Integer component of zero = 0 */

- if (mp_is_zero(x)) {

- mp_set_from_mp(x, z);

- return;

- }

-

- /* If all fractional then no integer component */

- if (x->exponent <= 0) {

- mp_set_from_integer(0, z);

- return;

- }

-

- is_negative = mp_is_negative(x);

-

- /* Clear fraction */

- mp_set_from_mp(x, z);

- for (i = z->exponent; i < MP_SIZE; i++) {

- if (z->fraction[i])

- have_fraction = true;

- z->fraction[i] = 0;

- }

- z->im_sign = 0;

- z->im_exponent = 0;

- memset(z->im_fraction, 0, sizeof(int) * MP_SIZE);

-

- if (have_fraction && is_negative)

- mp_add_integer(z, -1, z);

- MPNumber f;

-

- mp_floor(x, z);

- mp_fractional_component(x, &f);

- if (mp_is_zero(&f))

- return;

- mp_add_integer(z, 1, z);

- MPNumber f, one;

- bool do_floor;

-

- do_floor = !mp_is_negative(x);

-

- mp_fractional_component(x, &f);

- mp_multiply_integer(&f, 2, &f);

- mp_abs(&f, &f);

- mp_set_from_integer(1, &one);

- if (mp_is_greater_equal(&f, &one))

- do_floor = !do_floor;

-

- if (do_floor)

- mp_floor(x, z);

- else

- mp_ceiling(x, z);

-mp_compare_mp_to_mp(const MPNumber *x, const MPNumber *y)

- int i;

-

- if (x->sign != y->sign) {

- if (x->sign > y->sign)

- return 1;

- else

- return -1;

- }

-

- /* x = y = 0 */

- if (mp_is_zero(x))

- return 0;

-

- /* See if numbers are of different magnitude */

- if (x->exponent != y->exponent) {

- if (x->exponent > y->exponent)

- return x->sign;

- else

- return -x->sign;

- }

-

- /* Compare fractions */

- for (i = 0; i < MP_SIZE; i++) {

- if (x->fraction[i] == y->fraction[i])

- continue;

-

- if (x->fraction[i] > y->fraction[i])

- return x->sign;

- else

- return -x->sign;

- }

-

- /* x = y */

- return 0;

-

- int i, ie;

- MPNumber t;

-

- /* x/0 */

- if (mp_is_zero(y)) {

- /* Translators: Error displayed attempted to divide by zero */

- mperr(_("Division by zero is undefined"));

- mp_set_from_integer(0, z);

- return;

- }

-

- /* 0/y = 0 */

- if (mp_is_zero(x)) {

- mp_set_from_integer(0, z);

- return;

- }

-

- /* z = x × y⁻¹ */

- /* FIXME: Set exponent to zero to avoid overflow in mp_multiply??? */

- mp_reciprocal(y, &t);

- ie = t.exponent;

- t.exponent = 0;

- i = t.fraction[0];

- mp_multiply(x, &t, z);

- mp_ext(i, z->fraction[0], z);

- z->exponent += ie;

-}

-

-

-static void

-mp_divide_integer_real(const MPNumber *x, int64_t y, MPNumber *z)

-{

- int c, i, k, b2, c2, j1, j2;

- MPNumber x_copy;

-

- /* x/0 */

- if (y == 0) {

-

- /* 0/y = 0 */

- if (mp_is_zero(x)) {

- mp_set_from_integer(0, z);

- return;

- }

-

- /* Division by -1 or 1 just changes sign */

- if (y == 1 || y == -1) {

- if (y < 0)

- mp_invert_sign(x, z);

- else

- mp_set_from_mp(x, z);

- return;

- }

-

- /* Copy x as z may also refer to x */

- mp_set_from_mp(x, &x_copy);

- mp_set_from_integer(0, z);

-

- if (y < 0) {

- y = -y;

- z->sign = -x_copy.sign;

- }

- else

- z->sign = x_copy.sign;

- z->exponent = x_copy.exponent;

-

- c = 0;

- i = 0;

-

- /* IF y*B NOT REPRESENTABLE AS AN INTEGER HAVE TO SIMULATE

- * LONG DIVISION. ASSUME AT LEAST 16-BIT WORD.

- */

-

- /* Computing MAX */

- b2 = max(MP_BASE << 3, 32767 / MP_BASE);

- if (y < b2) {

- int kh, r1;

-

- /* LOOK FOR FIRST NONZERO DIGIT IN QUOTIENT */

- do {

- c = MP_BASE * c;

- if (i < MP_T)

- c += x_copy.fraction[i];

- i++;

- r1 = c / y;

- if (r1 < 0)

- goto L210;

- } while (r1 == 0);

-

- /* ADJUST EXPONENT AND GET T+4 DIGITS IN QUOTIENT */

- z->exponent += 1 - i;

- z->fraction[0] = r1;

- c = MP_BASE * (c - y * r1);

- kh = 1;

- if (i < MP_T) {

- kh = MP_T + 1 - i;

- for (k = 1; k < kh; k++) {

- c += x_copy.fraction[i];

- z->fraction[k] = c / y;

- c = MP_BASE * (c - y * z->fraction[k]);

- i++;

- }

- if (c < 0)

- goto L210;

- }

-

- for (k = kh; k < MP_T + 4; k++) {

- z->fraction[k] = c / y;

- c = MP_BASE * (c - y * z->fraction[k]);

- }

- if (c < 0)

- goto L210;

-

- mp_normalize(z);

- return;

- }

-

- /* HERE NEED SIMULATED DOUBLE-PRECISION DIVISION */

- j1 = y / MP_BASE;

- j2 = y - j1 * MP_BASE;

-

- /* LOOK FOR FIRST NONZERO DIGIT */

- c2 = 0;

- do {

- c = MP_BASE * c + c2;

- c2 = i < MP_T ? x_copy.fraction[i] : 0;

- i++;

- } while (c < j1 || (c == j1 && c2 < j2));

-

- /* COMPUTE T+4 QUOTIENT DIGITS */

- z->exponent += 1 - i;

- i--;

-

- /* MAIN LOOP FOR LARGE ABS(y) CASE */

- for (k = 1; k <= MP_T + 4; k++) {

- int ir, iq, iqj;

-

- /* GET APPROXIMATE QUOTIENT FIRST */

- ir = c / (j1 + 1);

-

- /* NOW REDUCE SO OVERFLOW DOES NOT OCCUR */

- iq = c - ir * j1;

- if (iq >= b2) {

- /* HERE IQ*B WOULD POSSIBLY OVERFLOW SO INCREASE IR */

- ++ir;

- iq -= j1;

- }

-

- iq = iq * MP_BASE - ir * j2;

- if (iq < 0) {

- /* HERE IQ NEGATIVE SO IR WAS TOO LARGE */

- ir--;

- iq += y;

- }

-

- if (i < MP_T)

- iq += x_copy.fraction[i];

- i++;

- iqj = iq / y;

-

- /* R(K) = QUOTIENT, C = REMAINDER */

- z->fraction[k - 1] = iqj + ir;

- c = iq - y * iqj;

-

- if (c < 0)

- goto L210;

- }

-

- mp_normalize(z);

-

-L210:

- /* CARRY NEGATIVE SO OVERFLOW MUST HAVE OCCURRED */

- mperr("*** INTEGER OVERFLOW IN MP_DIVIDE_INTEGER, B TOO LARGE ***");

- mp_set_from_integer(0, z);

-

- if (mp_is_complex(x)) {

- MPNumber re_z, im_z;

-

- mp_real_component(x, &re_z);

- mp_imaginary_component(x, &im_z);

- mp_divide_integer_real(&re_z, y, &re_z);

- mp_divide_integer_real(&im_z, y, &im_z);

- mp_set_from_complex(&re_z, &im_z, z);

- }

- else

- mp_divide_integer_real(x, y, z);

-}

- MPNumber t1, t2, t3;

-

- /* This fix is required for 1/3 repiprocal not being detected as an integer */

- /* Multiplication and division by 10000 is used to get around a

- * limitation to the "fix" for Sun bugtraq bug #4006391 in the

- * mp_floor() routine in mp.c, when the exponent is less than 1.

- */

- mp_set_from_integer(10000, &t3);

- mp_multiply(x, &t3, &t1);

- mp_divide(&t1, &t3, &t1);

- mp_floor(&t1, &t2);

- return mp_is_equal(&t1, &t2);

-

- /* Correct way to check for integer */

- /*int i;

-

- // Zero is an integer

- if (mp_is_zero(x))

- return true;

-

- // Fractional

- if (x->exponent <= 0)

- return false;

-

- // Look for fractional components

- for (i = x->exponent; i < MP_SIZE; i++) {

- if (x->fraction[i] != 0)

- return false;

- }

-

- return true;*/

- return x->sign >= 0 && mp_is_integer(x);

- return x->sign > 0 && mp_is_integer(x);

- return x->im_sign != 0;

- return mp_compare_mp_to_mp(x, y) == 0;

-}

-

-

-/* Return e^x for |x| < 1 USING AN O(SQRT(T).M(T)) ALGORITHM

- * DESCRIBED IN - R. P. BRENT, THE COMPLEXITY OF MULTIPLE-

- * PRECISION ARITHMETIC (IN COMPLEXITY OF COMPUTATIONAL PROBLEM

- * SOLVING, UNIV. OF QUEENSLAND PRESS, BRISBANE, 1976, 126-165).

- * ASYMPTOTICALLY FASTER METHODS EXIST, BUT ARE NOT USEFUL

- * UNLESS T IS VERY LARGE. SEE COMMENTS TO MP_ATAN AND MPPIGL.

- */

-static void

-mp_exp(const MPNumber *x, MPNumber *z)

-{

- int i, q;

- float rlb;

- MPNumber t1, t2;

-

- /* e^0 = 1 */

- if (mp_is_zero(x)) {

- mp_set_from_integer(1, z);

- return;

- }

-

- /* Only defined for |x| < 1 */

- if (x->exponent > 0) {

- mperr("*** ABS(X) NOT LESS THAN 1 IN CALL TO MP_EXP ***");

- mp_set_from_integer(0, z);

- return;

- }

-

- mp_set_from_mp(x, &t1);

- rlb = log((float)MP_BASE);

-

- /* Compute approximately optimal q (and divide x by 2^q) */

- q = (int)(sqrt((float)MP_T * 0.48f * rlb) + (float) x->exponent * 1.44f * rlb);

-

- /* HALVE Q TIMES */

- if (q > 0) {

- int ib, ic;

-

- ib = MP_BASE << 2;

- ic = 1;

- for (i = 1; i <= q; ++i) {

- ic *= 2;

- if (ic < ib && ic != MP_BASE && i < q)

- continue;

- mp_divide_integer(&t1, ic, &t1);

- ic = 1;

- }

- }

-

- if (mp_is_zero(&t1)) {

- mp_set_from_integer(0, z);

- return;

- }

-

- /* Sum series, reducing t where possible */

- mp_set_from_mp(&t1, z);

- mp_set_from_mp(&t1, &t2);

- for (i = 2; MP_T + t2.exponent - z->exponent > 0; i++) {

- mp_multiply(&t1, &t2, &t2);

- mp_divide_integer(&t2, i, &t2);

- mp_add(&t2, z, z);

- if (mp_is_zero(&t2))

- break;

- }

-

- /* Apply (x+1)^2 - 1 = x(2 + x) for q iterations */

- for (i = 1; i <= q; ++i) {

- mp_add_integer(z, 2, &t1);

- mp_multiply(&t1, z, z);

- }

-

- mp_add_integer(z, 1, z);

-}

-

-

-static void

-mp_epowy_real(const MPNumber *x, MPNumber *z)

-{

- float r__1;

- int i, ix, xs, tss;

- float rx, rz;

- MPNumber t1, t2;

-

- /* e^0 = 1 */

- if (mp_is_zero(x)) {

- mp_set_from_integer(1, z);

- return;

- }

-

- /* If |x| < 1 use mp_exp */

- if (x->exponent <= 0) {

- mp_exp(x, z);

- return;

- }

-

- /* NOW SAFE TO CONVERT X TO REAL */

- rx = mp_cast_to_float(x);

-

- /* SAVE SIGN AND WORK WITH ABS(X) */

- xs = x->sign;

- mp_abs(x, &t2);

-

- /* GET FRACTIONAL AND INTEGER PARTS OF ABS(X) */

- ix = mp_cast_to_int(&t2);

- mp_fractional_component(&t2, &t2);

-

- /* ATTACH SIGN TO FRACTIONAL PART AND COMPUTE EXP OF IT */

- t2.sign *= xs;

- mp_exp(&t2, z);

-

- /* COMPUTE E-2 OR 1/E USING TWO EXTRA DIGITS IN CASE ABS(X) LARGE

- * (BUT ONLY ONE EXTRA DIGIT IF T < 4)

- */

- if (MP_T < 4)

- tss = MP_T + 1;

- else

- tss = MP_T + 2;

-

- /* LOOP FOR E COMPUTATION. DECREASE T IF POSSIBLE. */

- /* Computing MIN */

- mp_set_from_integer(xs, &t1);

-

- t2.sign = 0;

- for (i = 2 ; ; i++) {

- if (min(tss, tss + 2 + t1.exponent) <= 2)

- break;

-

- mp_divide_integer(&t1, i * xs, &t1);

- mp_add(&t2, &t1, &t2);

- if (mp_is_zero(&t1))

- break;

- }

-

- /* RAISE E OR 1/E TO POWER IX */

- if (xs > 0)

- mp_add_integer(&t2, 2, &t2);

- mp_xpowy_integer(&t2, ix, &t2);

-

- /* MULTIPLY EXPS OF INTEGER AND FRACTIONAL PARTS */

- mp_multiply(z, &t2, z);

-

- /* CHECK THAT RELATIVE ERROR LESS THAN 0.01 UNLESS ABS(X) LARGE

- * (WHEN EXP MIGHT OVERFLOW OR UNDERFLOW)

- */

- if (fabs(rx) > 10.0f)

- return;

-

- rz = mp_cast_to_float(z);

- r__1 = rz - exp(rx);

- if (fabs(r__1) < rz * 0.01f)

- return;

-

- /* THE FOLLOWING MESSAGE MAY INDICATE THAT

- * B**(T-1) IS TOO SMALL, OR THAT M IS TOO SMALL SO THE

- * RESULT UNDERFLOWED.

- */

- mperr("*** ERROR OCCURRED IN MP_EPOWY, RESULT INCORRECT ***");

- /* e^0 = 1 */

- if (mp_is_zero(x)) {

- mp_set_from_integer(1, z);

- return;

- }

-

- if (mp_is_complex(x)) {

- MPNumber x_real, r, theta;

-

- mp_real_component(x, &x_real);

- mp_imaginary_component(x, &theta);

-

- mp_epowy_real(&x_real, &r);

- mp_set_from_polar(&r, MP_RADIANS, &theta, z);

- }

- else

- mp_epowy_real(x, z);

-}

-

-

-/* RETURNS K = K/GCD AND L = L/GCD, WHERE GCD IS THE

- * GREATEST COMMON DIVISOR OF K AND L.

- * SAVE INPUT PARAMETERS IN LOCAL VARIABLES

- */

-void

-mp_gcd(int64_t *k, int64_t *l)

-{

- int64_t i, j;

-

- i = abs(*k);

- j = abs(*l);

- if (j == 0) {

- /* IF J = 0 RETURN (1, 0) UNLESS I = 0, THEN (0, 0) */

- *k = 1;

- *l = 0;

- if (i == 0)

- *k = 0;

- return;

- }

-

- /* EUCLIDEAN ALGORITHM LOOP */

- do {

- i %= j;

- if (i == 0) {

- *k = *k / j;

- *l = *l / j;

- return;

- }

- j %= i;

- } while (j != 0);

-

- /* HERE J IS THE GCD OF K AND L */

- *k = *k / i;

- *l = *l / i;

-

-mp_is_zero(const MPNumber *x)

- return x->sign == 0 && x->im_sign == 0;

-

- return x->sign < 0;

- return mp_compare_mp_to_mp(x, y) >= 0;

- return mp_compare_mp_to_mp(x, y) > 0;

- return mp_compare_mp_to_mp(x, y) <= 0;

-}

-

-

-/* RETURNS MP Y = LN(1+X) IF X IS AN MP NUMBER SATISFYING THE

- * CONDITION ABS(X) < 1/B, ERROR OTHERWISE.

- * USES NEWTONS METHOD TO SOLVE THE EQUATION

- * EXP1(-Y) = X, THEN REVERSES SIGN OF Y.

- */

-static void

-mp_lns(const MPNumber *x, MPNumber *z)

-{

- int t, it0;

- MPNumber t1, t2, t3;

-

- /* ln(1+0) = 0 */

- if (mp_is_zero(x)) {

- mp_set_from_integer(0, z);

- return;

- }

-

- /* Get starting approximation -ln(1+x) ~= -x + x^2/2 - x^3/3 + x^4/4 */

- mp_set_from_mp(x, &t2);

- mp_divide_integer(x, 4, &t1);

- mp_add_fraction(&t1, -1, 3, &t1);

- mp_multiply(x, &t1, &t1);

- mp_add_fraction(&t1, 1, 2, &t1);

- mp_multiply(x, &t1, &t1);

- mp_add_integer(&t1, -1, &t1);

- mp_multiply(x, &t1, z);

-

- /* Solve using Newtons method */

- it0 = t = 5;

- while(1)

- {

- int ts2, ts3;

-

- /* t3 = (e^t3 - 1) */

- /* z = z - (t2 + t3 + (t2 * t3)) */

- mp_epowy(z, &t3);

- mp_add_integer(&t3, -1, &t3);

- mp_multiply(&t2, &t3, &t1);

- mp_add(&t3, &t1, &t3);

- mp_add(&t2, &t3, &t3);

- mp_subtract(z, &t3, z);

- if (t >= MP_T)

- break;

-

- /* FOLLOWING LOOP COMPUTES NEXT VALUE OF T TO USE.

- * BECAUSE NEWTONS METHOD HAS 2ND ORDER CONVERGENCE,

- * WE CAN ALMOST DOUBLE T EACH TIME.

- */

- ts3 = t;

- t = MP_T;

- do {

- ts2 = t;

- t = (t + it0) / 2;

- } while (t > ts3);

- t = ts2;

- }

-

- /* CHECK THAT NEWTON ITERATION WAS CONVERGING AS EXPECTED */

- if (t3.sign != 0 && t3.exponent << 1 > it0 - MP_T) {

- mperr("*** ERROR OCCURRED IN MP_LNS, NEWTON ITERATION NOT CONVERGING PROPERLY ***");

- }

-

- z->sign = -z->sign;

-

-static void

-mp_ln_real(const MPNumber *x, MPNumber *z)

- int e, k;

- double rx, rlx;

- MPNumber t1, t2;

-

- /* LOOP TO GET APPROXIMATE LN(X) USING SINGLE-PRECISION */

- mp_set_from_mp(x, &t1);

- mp_set_from_integer(0, z);

- for(k = 0; k < 10; k++)

- {

- /* COMPUTE FINAL CORRECTION ACCURATELY USING MP_LNS */

- mp_add_integer(&t1, -1, &t2);

- if (mp_is_zero(&t2) || t2.exponent + 1 <= 0) {

- mp_lns(&t2, &t2);

- mp_add(z, &t2, z);

- return;

- }

-

- /* REMOVE EXPONENT TO AVOID FLOATING-POINT OVERFLOW */

- e = t1.exponent;

- t1.exponent = 0;

- rx = mp_cast_to_double(&t1);

- t1.exponent = e;

- rlx = log(rx) + e * log(MP_BASE);

- mp_set_from_double(-(double)rlx, &t2);

-

- /* UPDATE Z AND COMPUTE ACCURATE EXP OF APPROXIMATE LOG */

- mp_subtract(z, &t2, z);

- mp_epowy(&t2, &t2);

-

- /* COMPUTE RESIDUAL WHOSE LOG IS STILL TO BE FOUND */

- mp_multiply(&t1, &t2, &t1);

- }

-

- mperr("*** ERROR IN MP_LN, ITERATION NOT CONVERGING ***");

-

- if (mp_is_zero(x)) {

- if (mp_is_complex(x) || mp_is_negative(x)) {

- MPNumber r, theta, z_real;

-

- /* ln(re^iθ) = e^(ln(r)+iθ) */

- mp_abs(x, &r);

- mp_arg(x, MP_RADIANS, &theta);

-

- mp_ln_real(&r, &z_real);

- mp_set_from_complex(&z_real, &theta, z);

- }

- else

- mp_ln_real(x, z);

- MPNumber t1, t2;

-

- if (mp_is_zero(x)) {

- mp_ln(x, &t2);

- mp_divide(&t2, &t1, z);

-}

-

-

-bool

-mp_is_less_than(const MPNumber *x, const MPNumber *y)

-{

- return mp_compare_mp_to_mp(x, y) < 0;

-

-static void

-mp_multiply_real(const MPNumber *x, const MPNumber *y, MPNumber *z)

-{

- int c, i, xi;

- MPNumber r;

-

- /* x*0 = 0*y = 0 */

- if (x->sign == 0 || y->sign == 0) {

- mp_set_from_integer(0, z);

- return;

- }

-

- z->sign = x->sign * y->sign;

- z->exponent = x->exponent + y->exponent;

- memset(&r, 0, sizeof(MPNumber));

-

- /* PERFORM MULTIPLICATION */

- c = 8;

- for (i = 0; i < MP_T; i++) {

- int j;

-

- xi = x->fraction[i];

-

- /* FOR SPEED, PUT THE NUMBER WITH MANY ZEROS FIRST */

- if (xi == 0)

- continue;

-

- /* Computing MIN */

- for (j = 0; j < min(MP_T, MP_T + 3 - i); j++)

- r.fraction[i+j+1] += xi * y->fraction[j];

- c--;

- if (c > 0)

- continue;

-

- /* CHECK FOR LEGAL BASE B DIGIT */

- if (xi < 0 || xi >= MP_BASE) {

- mperr("*** ILLEGAL BASE B DIGIT IN CALL TO MP_MULTIPLY, POSSIBLE OVERWRITING PROBLEM ***");

- mp_set_from_integer(0, z);

- return;

- }

-

- /* PROPAGATE CARRIES AT END AND EVERY EIGHTH TIME,

- * FASTER THAN DOING IT EVERY TIME.

- */

- for (j = MP_T + 3; j >= 0; j--) {

- int ri = r.fraction[j] + c;

- if (ri < 0) {

- mperr("*** INTEGER OVERFLOW IN MP_MULTIPLY, B TOO LARGE ***");

- mp_set_from_integer(0, z);

- return;

- }

- c = ri / MP_BASE;

- r.fraction[j] = ri - MP_BASE * c;

- }

- if (c != 0) {

- mperr("*** ILLEGAL BASE B DIGIT IN CALL TO MP_MULTIPLY, POSSIBLE OVERWRITING PROBLEM ***");

- mp_set_from_integer(0, z);

- return;

- }

- c = 8;

- }

-

- if (c != 8) {

- if (xi < 0 || xi >= MP_BASE) {

- mperr("*** ILLEGAL BASE B DIGIT IN CALL TO MP_MULTIPLY, POSSIBLE OVERWRITING PROBLEM ***");

- mp_set_from_integer(0, z);

- return;

- }

-

- c = 0;

- for (i = MP_T + 3; i >= 0; i--) {

- int ri = r.fraction[i] + c;

- if (ri < 0) {

- mperr("*** INTEGER OVERFLOW IN MP_MULTIPLY, B TOO LARGE ***");

- mp_set_from_integer(0, z);

- return;

- }

- c = ri / MP_BASE;

- r.fraction[i] = ri - MP_BASE * c;

- }

-

- if (c != 0) {

- mperr("*** ILLEGAL BASE B DIGIT IN CALL TO MP_MULTIPLY, POSSIBLE OVERWRITING PROBLEM ***");

- mp_set_from_integer(0, z);

- return;

- }

- }

-

- /* Clear complex part */

- z->im_sign = 0;

- z->im_exponent = 0;

- memset(z->im_fraction, 0, sizeof(int) * MP_SIZE);

-

- /* NORMALIZE AND ROUND RESULT */

- // FIXME: Use stack variable because of mp_normalize brokeness

- for (i = 0; i < MP_SIZE; i++)

- z->fraction[i] = r.fraction[i];

- mp_normalize(z);

-}

-

-

- /* x*0 = 0*y = 0 */

- if (mp_is_zero(x) || mp_is_zero(y)) {

- mp_set_from_integer(0, z);

- return;

- }

-

- /* (a+bi)(c+di) = (ac-bd)+(ad+bc)i */

- if (mp_is_complex(x) || mp_is_complex(y)) {

- MPNumber real_x, real_y, im_x, im_y, t1, t2, real_z, im_z;

-

- mp_real_component(x, &real_x);

- mp_imaginary_component(x, &im_x);

- mp_real_component(y, &real_y);

- mp_imaginary_component(y, &im_y);

-

- mp_multiply_real(&real_x, &real_y, &t1);

- mp_multiply_real(&im_x, &im_y, &t2);

- mp_subtract(&t1, &t2, &real_z);

-

- mp_multiply_real(&real_x, &im_y, &t1);

- mp_multiply_real(&im_x, &real_y, &t2);

- mp_add(&t1, &t2, &im_z);

-

- mp_set_from_complex(&real_z, &im_z, z);

- }

- else {

- mp_multiply_real(x, y, z);

- }

-}

-

-

-static void

-mp_multiply_integer_real(const MPNumber *x, int64_t y, MPNumber *z)

-{

- int c, i;

- MPNumber x_copy;

-

- /* x*0 = 0*y = 0 */

- if (mp_is_zero(x) || y == 0) {

- mp_set_from_integer(0, z);

- return;

- }

-

- /* x*1 = x, x*-1 = -x */

- // FIXME: Why is this not working? mp_ext is using this function to do a normalization

- /*if (y == 1 || y == -1) {

- if (y < 0)

- mp_invert_sign(x, z);

- else

- mp_set_from_mp(x, z);

- return;

- }*/

-

- /* Copy x as z may also refer to x */

- mp_set_from_mp(x, &x_copy);

- mp_set_from_integer(0, z);

-

- if (y < 0) {

- y = -y;

- z->sign = -x_copy.sign;

- }

- else

- z->sign = x_copy.sign;

- z->exponent = x_copy.exponent + 4;

-

- /* FORM PRODUCT IN ACCUMULATOR */

- c = 0;

-

- /* IF y*B NOT REPRESENTABLE AS AN INTEGER WE HAVE TO SIMULATE

- * DOUBLE-PRECISION MULTIPLICATION.

- */

-

- /* Computing MAX */

- if (y >= max(MP_BASE << 3, 32767 / MP_BASE)) {

- int64_t j1, j2;

-

- /* HERE J IS TOO LARGE FOR SINGLE-PRECISION MULTIPLICATION */

- j1 = y / MP_BASE;

- j2 = y - j1 * MP_BASE;

-

- /* FORM PRODUCT */

- for (i = MP_T + 3; i >= 0; i--) {

- int64_t c1, c2, is, ix, t;

-

- c1 = c / MP_BASE;

- c2 = c - MP_BASE * c1;

- ix = 0;

- if (i > 3)

- ix = x_copy.fraction[i - 4];

-

- t = j2 * ix + c2;

- is = t / MP_BASE;

- c = j1 * ix + c1 + is;

- z->fraction[i] = t - MP_BASE * is;

- }

- }

- else

- {

- int64_t ri = 0;

-

- for (i = MP_T + 3; i >= 4; i--) {

- ri = y * x_copy.fraction[i - 4] + c;

- c = ri / MP_BASE;

- z->fraction[i] = ri - MP_BASE * c;

- }

-

- /* CHECK FOR INTEGER OVERFLOW */

- if (ri < 0) {

- mperr("*** INTEGER OVERFLOW IN mp_multiply_integer, B TOO LARGE ***");

- mp_set_from_integer(0, z);

- return;

- }

-

- /* HAVE TO TREAT FIRST FOUR WORDS OF R SEPARATELY */

- for (i = 3; i >= 0; i--) {

- int t;

-

- t = c;

- c = t / MP_BASE;

- z->fraction[i] = t - MP_BASE * c;

- }

- }

-

- /* HAVE TO SHIFT RIGHT HERE AS CARRY OFF END */

- while (c != 0) {

- int64_t t;

-

- for (i = MP_T + 3; i >= 1; i--)

- z->fraction[i] = z->fraction[i - 1];

- t = c;

- c = t / MP_BASE;

- z->fraction[0] = t - MP_BASE * c;

- z->exponent++;

- }

-

- if (c < 0) {

- mperr("*** INTEGER OVERFLOW IN mp_multiply_integer, B TOO LARGE ***");

- mp_set_from_integer(0, z);

- return;

- }

-

- z->im_sign = 0;

- z->im_exponent = 0;

- memset(z->im_fraction, 0, sizeof(int) * MP_SIZE);

- mp_normalize(z);

- if (mp_is_complex(x)) {

- MPNumber re_z, im_z;

- mp_real_component(x, &re_z);

- mp_imaginary_component(x, &im_z);

- mp_multiply_integer_real(&re_z, y, &re_z);

- mp_multiply_integer_real(&im_z, y, &im_z);

- mp_set_from_complex(&re_z, &im_z, z);

- }

- else

- mp_multiply_integer_real(x, y, z);

-}

-

-

-void

-mp_multiply_fraction(const MPNumber *x, int64_t numerator, int64_t denominator, MPNumber *z)

-{

- if (denominator == 0) {

- mperr(_("Division by zero is undefined"));

- mp_set_from_integer(0, z);

- return;

- }

-

- if (numerator == 0) {

- mp_set_from_integer(0, z);

- return;

- }

-

- /* Reduce to lowest terms */

- mp_gcd(&numerator, &denominator);

- mp_divide_integer(x, denominator, z);

- mp_multiply_integer(z, numerator, z);

- mp_set_from_mp(x, z);

- z->sign = -z->sign;

- z->im_sign = -z->im_sign;

-}

-

-

-// FIXME: Is r->fraction large enough? It seems to be in practise but it may be MP_T+4 instead of MP_T

-// FIXME: There is some sort of stack corruption/use of unitialised variables here. Some functions are

-// using stack variables as x otherwise there are corruption errors. e.g. "Cos(45) - 1/Sqrt(2) = -0"

-// (try in scientific mode)

-void

-mp_normalize(MPNumber *x)

-{

- int start_index;

-

- /* Find first non-zero digit */

- for (start_index = 0; start_index < MP_SIZE && x->fraction[start_index] == 0; start_index++);

-

- /* Mark as zero */

- if (start_index >= MP_SIZE) {

- x->sign = 0;

- x->exponent = 0;

- return;

- }

-

- /* Shift left so first digit is non-zero */

- if (start_index > 0) {

- int i;

-

- x->exponent -= start_index;

- for (i = 0; (i + start_index) < MP_SIZE; i++)

- x->fraction[i] = x->fraction[i + start_index];

- for (; i < MP_SIZE; i++)

- x->fraction[i] = 0;

- }

-}

-

-

-static void

-mp_pwr(const MPNumber *x, const MPNumber *y, MPNumber *z)

-{

- MPNumber t;

-

- /* (-x)^y imaginary */

- /* FIXME: Make complex numbers optional */

- /*if (x->sign < 0) {

- mperr(_("The power of negative numbers is only defined for integer exponents"));

- mp_set_from_integer(0, z);

- return;

- }*/

-

- /* 0^y = 0, 0^-y undefined */

- if (mp_is_zero(x)) {

- mp_set_from_integer(0, z);

- if (y->sign < 0)

- mperr(_("The power of zero is undefined for a negative exponent"));

- return;

- }

-

- /* x^0 = 1 */

- if (mp_is_zero(y)) {

- mp_set_from_integer(1, z);

- return;

- }

-

- mp_ln(x, &t);

- mp_multiply(y, &t, z);

- mp_epowy(z, z);

-}

-

-

-static void

-mp_reciprocal_real(const MPNumber *x, MPNumber *z)

-{

- MPNumber t1, t2;

- int it0, t;

-

- /* 1/0 invalid */

- if (mp_is_zero(x)) {

- mperr(_("Reciprocal of zero is undefined"));

- mp_set_from_integer(0, z);

- return;

- }

-

- /* Start by approximating value using floating point */

- mp_set_from_mp(x, &t1);

- t1.exponent = 0;

- mp_set_from_float(1.0f / mp_cast_to_float(&t1), &t1);

- t1.exponent -= x->exponent;

-

- it0 = t = 3;

- while(1) {

- int ts2, ts3;

-

- /* t1 = t1 - (t1 * ((x * t1) - 1)) (2*t1 - t1^2*x) */

- mp_multiply(x, &t1, &t2);

- mp_add_integer(&t2, -1, &t2);

- mp_multiply(&t1, &t2, &t2);

- mp_subtract(&t1, &t2, &t1);

- if (t >= MP_T)

- break;

-

- /* FOLLOWING LOOP ALMOST DOUBLES T (POSSIBLE

- * BECAUSE NEWTONS METHOD HAS 2ND ORDER CONVERGENCE).

- */

- ts3 = t;

- t = MP_T;

- do {

- ts2 = t;

- t = (t + it0) / 2;

- } while (t > ts3);

- t = min(ts2, MP_T);

- }

-

- /* RETURN IF NEWTON ITERATION WAS CONVERGING */

- if (t2.sign != 0 && (t1.exponent - t2.exponent) << 1 < MP_T - it0) {

- /* THE FOLLOWING MESSAGE MAY INDICATE THAT B**(T-1) IS TOO SMALL,

- * OR THAT THE STARTING APPROXIMATION IS NOT ACCURATE ENOUGH.

- */

- mperr("*** ERROR OCCURRED IN MP_RECIPROCAL, NEWTON ITERATION NOT CONVERGING PROPERLY ***");

- }

-

- mp_set_from_mp(&t1, z);

- if (mp_is_complex(x)) {

- MPNumber t1, t2;

- MPNumber real_x, im_x;

-

- mp_real_component(x, &real_x);

- mp_imaginary_component(x, &im_x);

-

- /* 1/(a+bi) = (a-bi)/(a+bi)(a-bi) = (a-bi)/(a²+b²) */

- mp_multiply(&real_x, &real_x, &t1);

- mp_multiply(&im_x, &im_x, &t2);

- mp_add(&t1, &t2, &t1);

- mp_reciprocal_real(&t1, z);

- mp_conjugate(x, &t1);

- mp_multiply(&t1, z, z);

- }

- else

- mp_reciprocal_real(x, z);

-

-static void

-mp_root_real(const MPNumber *x, int64_t n, MPNumber *z)

- float approximation;

- int ex, np, it0, t;

- MPNumber t1, t2;

-

- /* x^(1/1) = x */

- if (n == 1) {

- mp_set_from_mp(x, z);

- return;

- }

- /* x^(1/0) invalid */

- if (n == 0) {

- mperr(_("Root must be non-zero"));

- mp_set_from_integer(0, z);

- return;

- }

-

- np = abs(n);

-

- /* LOSS OF ACCURACY IF NP LARGE, SO ONLY ALLOW NP <= MAX (B, 64) */

- if (np > max(MP_BASE, 64)) {

- mperr("*** ABS(N) TOO LARGE IN CALL TO MP_ROOT ***");

- mp_set_from_integer(0, z);

- return;

-

- /* 0^(1/n) = 0 for positive n */

- if (mp_is_zero(x)) {

- mp_set_from_integer(0, z);

- if (n <= 0)

- mperr(_("Negative root of zero is undefined"));

- return;

-

- // FIXME: Imaginary root

- if (x->sign < 0 && np % 2 == 0) {

- mperr(_("nth root of negative number is undefined for even n"));

-

- /* DIVIDE EXPONENT BY NP */

- ex = x->exponent / np;

-

- /* Get initial approximation */

- mp_set_from_mp(x, &t1);

- t1.exponent = 0;

- approximation = exp(((float)(np * ex - x->exponent) * log((float)MP_BASE) -

- log((fabs(mp_cast_to_float(&t1))))) / (float)np);

- mp_set_from_float(approximation, &t1);

- t1.sign = x->sign;

- t1.exponent -= ex;

-

- /* MAIN ITERATION LOOP */

- it0 = t = 3;

- while(1) {

- int ts2, ts3;

-

- /* t1 = t1 - ((t1 * ((x * t1^np) - 1)) / np) */

- mp_xpowy_integer(&t1, np, &t2);

- mp_multiply(x, &t2, &t2);

- mp_add_integer(&t2, -1, &t2);

- mp_multiply(&t1, &t2, &t2);

- mp_divide_integer(&t2, np, &t2);

- mp_subtract(&t1, &t2, &t1);

-

- /* FOLLOWING LOOP ALMOST DOUBLES T (POSSIBLE BECAUSE

- * NEWTONS METHOD HAS 2ND ORDER CONVERGENCE).

- */

- if (t >= MP_T)

- break;

-

- ts3 = t;

- t = MP_T;

- do {

- ts2 = t;

- t = (t + it0) / 2;

- } while (t > ts3);

- t = min(ts2, MP_T);

- }

-

- /* NOW R(I2) IS X**(-1/NP)

- * CHECK THAT NEWTON ITERATION WAS CONVERGING

- */

- if (t2.sign != 0 && (t1.exponent - t2.exponent) << 1 < MP_T - it0) {

- /* THE FOLLOWING MESSAGE MAY INDICATE THAT B**(T-1) IS TOO SMALL,

- * OR THAT THE INITIAL APPROXIMATION OBTAINED USING ALOG AND EXP

- * IS NOT ACCURATE ENOUGH.

- */

- mperr("*** ERROR OCCURRED IN MP_ROOT, NEWTON ITERATION NOT CONVERGING PROPERLY ***");

- }

-

- if (n >= 0) {

- mp_xpowy_integer(&t1, n - 1, &t1);

- mp_multiply(x, &t1, z);

- return;

- }

-

- mp_set_from_mp(&t1, z);

-}

-

-

-void

-mp_root(const MPNumber *x, int64_t n, MPNumber *z)

-{

- if (!mp_is_complex(x) && mp_is_negative(x) && n % 2 == 1) {

- mp_abs(x, z);

- mp_root_real(z, n, z);

- mp_invert_sign(z, z);

- }

- else if (mp_is_complex(x) || mp_is_negative(x)) {

- MPNumber r, theta;

-

- mp_abs(x, &r);

- mp_arg(x, MP_RADIANS, &theta);

-

- mp_root_real(&r, n, &r);

- mp_divide_integer(&theta, n, &theta);

- mp_set_from_polar(&r, MP_RADIANS, &theta, z);

- mp_root_real(x, n, z);

- if (mp_is_zero(x))

- mp_set_from_integer(0, z);

- /* FIXME: Make complex numbers optional */

- /*else if (x->sign < 0) {

- mperr(_("Square root is undefined for negative values"));

- mp_set_from_integer(0, z);

- }*/

- else {

- MPNumber t;

-

- mp_root(x, -2, &t);

- mp_multiply(x, &t, z);

- mp_ext(t.fraction[0], z->fraction[0], z);

- }

- int i, value;

-

- if (mp_is_zero(x)) {

- mp_set_from_integer(1, z);

- if (!mp_is_natural(x)) {

- /* Translators: Error displayed when attempted take the factorial of a fractional number */

- mperr(_("Factorial is only defined for natural numbers"));

- mp_set_from_integer(0, z);

- return;

-

- /* Convert to integer - if couldn't be converted then the factorial would be too big anyway */

- value = mp_cast_to_int(x);

- mp_set_from_mp(x, z);

- for (i = 2; i < value; i++)

- mp_multiply_integer(z, i, z);

- MPNumber t1, t2;

-

- if (!mp_is_integer(x) || !mp_is_integer(y)) {

- /* Translators: Error displayed when attemping to do a modulus division on non-integer numbers */

- mp_divide(x, y, &t1);

- mp_floor(&t1, &t1);

- mp_multiply(&t1, y, &t2);

- mp_subtract(x, &t2, z);

- if ((mp_is_less_than(y, &t1) && mp_is_greater_than(z, &t1)) || mp_is_less_than(z, &t1))

- if (mp_is_integer(y)) {

- mp_xpowy_integer(x, mp_cast_to_int(y), z);

- } else {

- MPNumber reciprocal;

- mp_root(x, mp_cast_to_int(&reciprocal), z);

- else

- mp_pwr(x, y, z);

- int i;

- MPNumber t;

-

- if (mp_is_zero(x) && n < 0) {

- /* Translators: Error displayed when attempted to raise 0 to a negative exponent */

- /* x^0 = 1 */

- if (n == 0) {

- mp_set_from_integer(1, z);

- return;

- }

-

- /* 0^n = 0 */

- if (mp_is_zero(x)) {

- mp_set_from_integer(0, z);

- return;

- }

-

- /* x^1 = x */

- if (n == 1) {

- mp_set_from_mp(x, z);

- return;

- }

-

- if (n < 0) {

- mp_reciprocal(x, &t);

- n = -n;

- }

- else

- mp_set_from_mp(x, &t);

-

- /* Multply x n times */

- // FIXME: Can do mp_multiply(z, z, z) until close to answer (each call doubles number of multiples) */

- mp_set_from_integer(1, z);

- for (i = 0; i < n; i++)

- mp_multiply(z, &t, z);

-

- MPNumber value, tmp, divisor, root;

- if (mp_is_zero(&value)) {

- if (mp_is_equal(&value, &tmp)) {

- while (TRUE) {

- if (mp_is_integer(&tmp)) {

- value = tmp;

- } else {

- break;

- while (mp_is_less_equal(&divisor, &root)) {

- if (mp_is_integer(&tmp)) {

- value = tmp;

- } else {

- divisor = tmp;

- if (mp_is_greater_than(&value, &tmp)) {

- } else {

- if (mp_is_negative(x)) {

-/* Size of the multiple precision values */

-#define MP_SIZE 1000

-/* Base for numbers */

-#define MP_BASE 10000

-/* Object for a high precision floating point number representation

- *

- * x = sign * (MP_BASE^(exponent-1) + MP_BASE^(exponent-2) + ...)

- */

- /* Sign (+1, -1) or 0 for the value zero */

- int sign, im_sign;

-

- /* Exponent (to base MP_BASE) */

- int exponent, im_exponent;

-

- /* Normalized fraction */

- int fraction[MP_SIZE], im_fraction[MP_SIZE];

-const char *mp_get_error(void);

-void mp_clear_error(void);

-int mp_compare_mp_to_mp(const MPNumber *x, const MPNumber *y);

-/* Sets z = x + numerator ÷ denominator */

-void mp_add_fraction(const MPNumber *x, int64_t numerator, int64_t denominator, MPNumber *z);

-

-/* Sets z = x × numerator ÷ denominator */

-void mp_multiply_fraction(const MPNumber *x, int64_t numerator, int64_t denominator, MPNumber *z);

-

-float mp_cast_to_float(const MPNumber *x);

-double mp_cast_to_double(const MPNumber *x);

-int64_t mp_cast_to_int(const MPNumber *x);

-uint64_t mp_cast_to_unsigned_int(const MPNumber *x);

- free(ans);

-

- free(val);

- MPNumber tmp;

- MPNumber* ans;

- ans = (MPNumber *) malloc(sizeof(MPNumber));

- free(ans);

- free(ans);

- free(ans);

- MPNumber tmp;

- MPNumber* ans;

- ans = (MPNumber *) malloc(sizeof(MPNumber));

- free(ans);

- free(ans);

- MPNumber value;

- MPNumber t;

- MPNumber* ans;

- ans = (MPNumber*) malloc(sizeof(MPNumber));

- MPNumber value;

- MPNumber t;

- MPNumber* ans;

- ans = (MPNumber*) malloc(sizeof(MPNumber));

- MPNumber* ans;

- ans = (MPNumber*) malloc(sizeof(MPNumber));

- MPNumber* tmp;

- MPNumber* ans;

- tmp = (MPNumber*) malloc(sizeof(MPNumber));

- ans = (MPNumber*) malloc(sizeof(MPNumber));

- free(tmp);

- free(ans);

- free(val);

- free(tmp);

- free(ans);

- free(tmp);

- free(ans);

- free(val);

- free(val);

- free(tmp);

- MPNumber* tmp;

- MPNumber* ans;

- tmp = (MPNumber*) malloc(sizeof(MPNumber));

- ans = (MPNumber*) malloc(sizeof(MPNumber));

- free(tmp);

- free(ans);

- free(tmp);

- free(ans);

- free(tmp);

- free(ans);

- free(val);

- free(val);

- free(tmp);

- MPNumber* ans;

- ans = (MPNumber*) malloc(sizeof(MPNumber));

- free(ans);

- free(val);

- MPNumber* ans;

- ans = (MPNumber*) malloc(sizeof(MPNumber));

- MPNumber* ans;

- ans = (MPNumber*) malloc(sizeof(MPNumber));

- free(ans);

- free(val);

- MPNumber* ans;

- ans = (MPNumber*) malloc(sizeof(MPNumber));

- free(ans);

- free(val);

- MPNumber* ans;

- ans = (MPNumber*) malloc(sizeof(MPNumber));

- free(ans);

- free(val);

- MPNumber* ans;

- ans = (MPNumber*) malloc(sizeof(MPNumber));

- free(ans);

- free(val);

- MPNumber* ans;

- ans = (MPNumber*) malloc(sizeof(MPNumber));

- free(ans);

- free(val);

- MPNumber* ans;

- ans = (MPNumber*) malloc(sizeof(MPNumber));

- free(ans);

- free(val);

- MPNumber* ans;

- ans = (MPNumber*) malloc(sizeof(MPNumber));

- free(ans);

- free(val);

- MPNumber* ans;

- ans = (MPNumber*) malloc(sizeof(MPNumber));

- free(val);

- free(pow);

- free(ans);

- free(val);

- free(pow);

- MPNumber* ans;

- ans = (MPNumber*) malloc(sizeof(MPNumber));

- free(ans);

- free(val);

- MPNumber* ans;

- ans = (MPNumber*) malloc(sizeof(MPNumber));

- free(ans);

- free(val);

- MPNumber* ans;

- ans = (MPNumber*) malloc(sizeof(MPNumber));

- free(ans);

- free(val);

- MPNumber* ans;

- ans = (MPNumber*) malloc(sizeof(MPNumber));

- free(left);

- free(right);

- free(ans);

- free(left);

- free(right);

- MPNumber* ans;

- ans = (MPNumber*) malloc(sizeof(MPNumber));

- free(left);

- free(right);

- free(left);

- free(right);

- MPNumber* ans;

- ans = (MPNumber*) malloc(sizeof(MPNumber));

- free(left);

- free(right);

- free(ans);

- free(left);

- free(right);

- MPNumber* ans;

- ans = (MPNumber*) malloc(sizeof(MPNumber));

- free(left);

- free (right);

- free(left);

- free(right);

- MPNumber* ans;

- ans = (MPNumber*) malloc(sizeof(MPNumber));

- free(left);

- free(right);

- free(left);

- free(right);

- MPNumber* ans;

- ans = (MPNumber*) malloc(sizeof(MPNumber));

- free(val);

- free(per);

- free(ans);

- free(val);

- free(per);

- MPNumber* ans;

- ans = (MPNumber*) malloc(sizeof(MPNumber));

- free(val);

- free(per);

- free(val);

- free(per);

- MPNumber* ans;

- ans = (MPNumber*) malloc(sizeof(MPNumber));

- free(ans);

- free(val);

- MPNumber* ans;

- ans = (MPNumber*) malloc(sizeof(MPNumber));

- free(ans);

- free(ans);

- free(val);

- free(val);

- MPNumber* ans;

- ans = (MPNumber*) malloc(sizeof(MPNumber));

- free(left);

- free(right);

- free(ans);

- free(left);

- free(right);

- MPNumber* ans;

- ans = (MPNumber*) malloc(sizeof(MPNumber));

- free(left);

- free(right);

- free(ans);

- free(left);

- free(right);

- MPNumber* ans;

- ans = (MPNumber*) malloc(sizeof(MPNumber));

- free(left);

- free(right);

- free(left);

- free(right);

- MPNumber* ans;

- ans = (MPNumber*) malloc(sizeof(MPNumber));

- free(ans);

- <property name="upper">15</property>

- MPNumber result;

- test("0.1!", "", PARSER_ERR_MP);

-#include "mp-private.h"

- int i, j;

-

- printf("sign=%d exponent=%d fraction=%d", x->sign, x->exponent, x->fraction[0]);

- for (i = 1; i < MP_SIZE; i++) {

- for (j = i; j < MP_SIZE && x->fraction[j] == 0; j++);

- if (j == MP_SIZE) {

- printf(",...");

- break;

- }

- printf(",%d", x->fraction[i]);

- }

- MPNumber t;

- MPNumber t;

- printf("base=%d\n", MP_BASE);

- MPNumber zero, one, minus_one;

-

- MPNumber t;

-

-

- MPNumber result;

- int i, j;

-

- printf("sign=%d exponent=%d fraction=%d", x->sign, x->exponent, x->fraction[0]);

-

- for (i = 1; i < MP_SIZE; i++)

- {

- for (j = i; j < MP_SIZE && x->fraction[j] == 0; j++)

- {

- /* nothing*/

- }

-

- if (j == MP_SIZE)

- {

- printf(",...");

- break;

- }

-

- printf(",%d", x->fraction[i]);

- }

- MPNumber t;

- MPNumber t;

-#include "mp-private.h"

- printf("base=%d\n", MP_BASE);

-

- MPNumber zero, one, minus_one;

-

- mp_set_from_integer(0, &zero);

- mp_set_from_integer(1, &one);

- mp_set_from_integer(-1, &minus_one);