mp.c

Bug Summary

File:	mp.c
Warning:	line 365, column 17 Access of 'x_copy.fraction' at an overflowing index, while it holds only 1000 'int' elements
Annotated Source Code

Press '?' to see keyboard shortcuts
Show analyzer invocation
clang -cc1 -cc1 -triple x86_64-pc-linux-gnu -analyze -disable-free -clear-ast-before-backend -disable-llvm-verifier -discard-value-names -main-file-name mp.c -analyzer-checker=core -analyzer-checker=apiModeling -analyzer-checker=unix -analyzer-checker=deadcode -analyzer-checker=security.insecureAPI.UncheckedReturn -analyzer-checker=security.insecureAPI.getpw -analyzer-checker=security.insecureAPI.gets -analyzer-checker=security.insecureAPI.mktemp -analyzer-checker=security.insecureAPI.mkstemp -analyzer-checker=security.insecureAPI.vfork -analyzer-checker=nullability.NullPassedToNonnull -analyzer-checker=nullability.NullReturnedFromNonnull -analyzer-output plist -w -setup-static-analyzer -mrelocation-model pic -pic-level 2 -pic-is-pie -mframe-pointer=all -fmath-errno -ffp-contract=on -fno-rounding-math -mconstructor-aliases -funwind-tables=2 -target-cpu x86-64 -tune-cpu generic -debugger-tuning=gdb -fdebug-compilation-dir=/rootdir/src -fcoverage-compilation-dir=/rootdir/src -resource-dir /usr/lib/llvm-21/lib/clang/21 -D HAVE_CONFIG_H -I . -I .. -D VERSION="2.1.0" -D LOCALE_DIR="/usr/share/locale" -D GETTEXT_PACKAGE="cafe-calc" -I /usr/include/ctk-3.0 -I /usr/include/pango-1.0 -I /usr/include/cairo -I /usr/include/cdk-pixbuf-2.0 -I /usr/include/at-spi2-atk/2.0 -I /usr/include/at-spi-2.0 -I /usr/include/atk-1.0 -I /usr/include/dbus-1.0 -I /usr/lib/x86_64-linux-gnu/dbus-1.0/include -I /usr/include/fribidi -I /usr/include/pixman-1 -I /usr/include/harfbuzz -I /usr/include/freetype2 -I /usr/include/libpng16 -I /usr/include/gio-unix-2.0 -I /usr/include/libmount -I /usr/include/blkid -I /usr/include/libxml2 -I /usr/include/glib-2.0 -I /usr/lib/x86_64-linux-gnu/glib-2.0/include -I /usr/include/sysprof-6 -internal-isystem /usr/lib/llvm-21/lib/clang/21/include -internal-isystem /usr/local/include -internal-isystem /usr/lib/gcc/x86_64-linux-gnu/15/../../../../x86_64-linux-gnu/include -internal-externc-isystem /usr/include/x86_64-linux-gnu -internal-externc-isystem /include -internal-externc-isystem /usr/include -ferror-limit 19 -fgnuc-version=4.2.1 -fskip-odr-check-in-gmf -analyzer-checker deadcode.DeadStores -analyzer-checker security.ArrayBound -analyzer-checker unix.cstring.NotNullTerminated -analyzer-checker alpha.deadcode.UnreachableCode -analyzer-checker alpha.core.CastToStruct -analyzer-checker alpha.security.ReturnPtrRange -analyzer-checker alpha.unix.SimpleStream -analyzer-checker alpha.unix.cstring.BufferOverlap -analyzer-checker alpha.unix.cstring.OutOfBounds -analyzer-checker alpha.core.FixedAddr -analyzer-output=html -faddrsig -D__GCC_HAVE_DWARF2_CFI_ASM=1 -o /rootdir/html-report/2026-08-16-122727-15300-1 -x c mp.c
1/*
* Copyright (C) 1987-2008 Sun Microsystems, Inc. All Rights Reserved.
* Copyright (C) 2008-2011 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.
*/

12#include <stdlib.h>
13#include <stdio.h>
14#include <math.h>
15#include <errno(*__errno_location ()).h>

17#include "mp.h"
18#include "mp-private.h"

20static MPNumber eulers_number;
21static gboolean have_eulers_number = FALSE(0);

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

25char *mp_error = NULL((void*)0);

27/*  THIS ROUTINE IS CALLED WHEN AN ERROR CONDITION IS ENCOUNTERED, AND
*  AFTER A MESSAGE HAS BEEN WRITTEN TO STDERR.
*/
30void
31mperr(const char *format, ...)
32{
  char text[1024];
  va_list args;

  va_start(args, format)__builtin_va_start(args, format);
  vsnprintf(text, 1024, format, args);
  va_end(args)__builtin_va_end(args);

  if (mp_error)
      free(mp_error);
  mp_error = strdup(text);
43}


46const char *
47mp_get_error()
48{
  return mp_error;
50}


53void mp_clear_error()
54{
  if (mp_error)
      free(mp_error);
  mp_error = NULL((void*)0);
58}


61/*  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.
*/
65static void
66mp_ext(int i, int j, MPNumber *x)
67{
  int q, s;

  if (mp_is_zero(x) || MP_T100 <= 2 || i == 0)
      return;

  /* COMPUTE MAXIMUM POSSIBLE ERROR IN THE LAST PLACE */
  q = (j + 1) / i + 1;
  s = MP_BASE10000 * x->fraction[MP_T100 - 2] + x->fraction[MP_T100 - 1];

  /* SET LAST TWO DIGITS TO ZERO */
  if (s <= q) {
      x->fraction[MP_T100 - 2] = 0;
      x->fraction[MP_T100 - 1] = 0;
      return;
  }

  if (s + q < MP_BASE10000 * MP_BASE10000)
      return;

  /* ROUND UP HERE */
  x->fraction[MP_T100 - 2] = MP_BASE10000 - 1;
  x->fraction[MP_T100 - 1] = MP_BASE10000;

  /* NORMALIZE X (LAST DIGIT B IS OK IN MP_MULTIPLY_INTEGER) */
  mp_multiply_integer(x, 1, x);
93}


96void
97mp_get_eulers(MPNumber *z)
98{
  if (!have_eulers_number) {
      MPNumber t;
      mp_set_from_integer(1, &t);
      mp_epowy(&t, &eulers_number);
      have_eulers_number = TRUE(!(0));
  }
  mp_set_from_mp(&eulers_number, z);
106}


109void
110mp_get_i(MPNumber *z)
111{
  mp_set_from_integer(0, z);
  z->im_sign = 1;
  z->im_exponent = 1;
  z->im_fraction[0] = 1;
116}


119void
120mp_abs(const MPNumber *x, MPNumber *z)
121{
  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;
  }
138}


141void
142mp_arg(const MPNumber *x, MPAngleUnit unit, MPNumber *z)
143{
  MPNumber x_real, x_im, pi;

  if (mp_is_zero(x)) {
      /* Translators: Error display when attempting to take argument of zero */
      mperr(_("Argument not defined for zero")gettext ("Argument not defined for zero"));
      mp_set_from_integer(0, z);
      return;
  }

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

  convert_from_radians(z, unit, z);
184}


187void
188mp_conjugate(const MPNumber *x, MPNumber *z)
189{
  mp_set_from_mp(x, z);
  z->im_sign = -z->im_sign;
192}


195void
196mp_real_component(const MPNumber *x, MPNumber *z)
197{
  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_SIZE1000);
204}


207void
208mp_imaginary_component(const MPNumber *x, MPNumber *z)
209{
  /* 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_SIZE1000);

  /* Clear (old) imaginary component */
  z->im_sign = 0;
  z->im_exponent = 0;
  memset(z->im_fraction, 0, sizeof(int) * MP_SIZE1000);
219}


222static void
223mp_add_real(const MPNumber *x, int y_sign, const MPNumber *y, MPNumber *z)
224{
  int sign_prod, i, c;
  int exp_diff, med;
  bool_Bool x_largest = false0;
  const int *big_fraction, *small_fraction;
  MPNumber x_copy, y_copy;

  /* 0 + y = y */
  if (mp_is_zero(x)) {
12
←
Assuming the condition is false→
13
←
Taking false branch→
      mp_set_from_mp(y, z);
      z->sign = y_sign;
      return;
  }
  /* x + 0 = x */
  else if (mp_is_zero(y)) {
14
←
Assuming the condition is false→
15
←
Taking false branch→
      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) {
16
←
Assuming 'exp_diff' is >= 0→
17
←
Taking false branch→
      x_largest = false0;
  } else if (exp_diff > 0) {
18
←
Assuming 'exp_diff' is <= 0→
19
←
Taking false branch→
      x_largest = true1;
  } else {
      /* EXPONENTS EQUAL SO COMPARE SIGNS, THEN FRACTIONS IF NEC. */
      if (sign_prod < 0) {
20
←
Assuming 'sign_prod' is < 0→
21
←
Taking true branch→
          /* Signs are not equal.  find out which mantissa is larger. */
          int j;
          for (j = 0; j < MP_T100; j++) {
22
←
Loop condition is true.  Entering loop body→
              int i = x->fraction[j] - y->fraction[j];
              if (i == 0)
23
←
Assuming 'i' is not equal to 0→
24
←
Taking false branch→
                  continue;
              if (i < 0)
25
←
Assuming 'i' is >= 0→
26
←
Taking false branch→
                  x_largest = false0;
              else if (i26.1
'i' is > 0
 > 0)
27
←
Taking true branch→
                  x_largest = true1;
              break;
          }

          /* Both mantissas equal, so result is zero. */
          if (j28.1
'j' is < MP_T
 >= MP_T100) {
28
←
 Execution continues on line 267→
29
←
Taking false branch→
              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_largest29.1
'x_largest' is true
) {
30
←
Taking true branch→
      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--)
31
←
Assuming 'i' is < 'med'→
32
←
Loop condition is false. Execution continues on line 294→
      z->fraction[MP_T100 + i] = 0;

  if (sign_prod32.1
'sign_prod' is < 0
 >= 0) {
33
←
Taking false branch→
      /* This is probably insufficient overflow detection, but it makes us
       * not crash at least.
       */
      if (MP_T100 + 3 < med) {
          mperr(_("Overflow: the result couldn't be calculated")gettext ("Overflow: the result couldn't be calculated"));
          mp_set_from_integer(0, z);
          return;
      }

      /* HERE DO ADDITION, EXPONENT(Y) >= EXPONENT(X) */
      for (i = MP_T100 + 3; i >= MP_T100; 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_BASE10000) {
              /* NO CARRY GENERATED HERE */
              z->fraction[i] = c;
              c = 0;
          } else {
              /* CARRY GENERATED HERE */
              z->fraction[i] = c - MP_BASE10000;
              c = 1;
          }
      }

      for (; i >= 0; i--)
      {
          c = big_fraction[i] + c;
          if (c < MP_BASE10000) {
              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_T100 + 3; i > 0; i--)
              z->fraction[i] = z->fraction[i - 1];
          z->fraction[0] = 1;
          z->exponent++;
      }
  }
  else  {
      c = 0;
      for (i = MP_T100 + med - 1; i >= MP_T100; i--) {
34
←
Assuming 'i' is < MP_T→
35
←
Loop condition is false. Execution continues on line 364→
          /* 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_BASE10000;
          }
      }

      for(; i >= med; i--) {
36
←
Assuming 'i' is >= 'med'→
37
←
Loop condition is true.  Entering loop body→
          c = big_fraction[i] + c - small_fraction[i - med];
38
←
Access of 'x_copy.fraction' at an overflowing index, while it holds only 1000 'int' elements
          if (c >= 0) {
              /* NO BORROW GENERATED HERE */
              z->fraction[i] = c;
              c = 0;
          } else {
              /* BORROW GENERATED HERE */
              z->fraction[i] = c + MP_BASE10000;
              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_BASE10000;
          c = -1;
      }
  }

  mp_normalize(z);
397}


400static void
401mp_add_with_sign(const MPNumber *x, int y_sign, const MPNumber *y, MPNumber *z)
402{
  if (mp_is_complex(x) || mp_is_complex(y)) {
10
←
Taking false branch→
      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);
11
←
Calling 'mp_add_real'→
418}


421void
422mp_add(const MPNumber *x, const MPNumber *y, MPNumber *z)
423{
  mp_add_with_sign(x, 1, y, z);
425}


428void
429mp_add_integer(const MPNumber *x, int64_t y, MPNumber *z)
430{
  MPNumber t;
  mp_set_from_integer(y, &t);
  mp_add(x, &t, z);
434}


437void
438mp_add_fraction(const MPNumber *x, int64_t i, int64_t j, MPNumber *y)
439{
  MPNumber t;
  mp_set_from_fraction(i, j, &t);
  mp_add(x, &t, y);
443}


446void
447mp_subtract(const MPNumber *x, const MPNumber *y, MPNumber *z)
448{
  mp_add_with_sign(x, -1, y, z);
9
←
Calling 'mp_add_with_sign'→
450}


453void
454mp_sgn(const MPNumber *x, MPNumber *z)
455{
  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);
462}

464void
465mp_integer_component(const MPNumber *x, MPNumber *z)
466{
  int i;

  /* Clear fraction */
  mp_set_from_mp(x, z);
  for (i = z->exponent; i < MP_SIZE1000; i++)
      z->fraction[i] = 0;
  z->im_sign = 0;
  z->im_exponent = 0;
  memset(z->im_fraction, 0, sizeof(int) * MP_SIZE1000);
476}

478void
479mp_fractional_component(const MPNumber *x, MPNumber *z)
480{
  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_SIZE1000 && x->fraction[i] == 0; i++)
      shift++;
  z->sign = x->sign;
  z->exponent = x->exponent - shift;
  for (i = 0; i < MP_SIZE1000; i++) {
      if (i + shift >= MP_SIZE1000)
          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_SIZE1000);
513}


516void
517mp_fractional_part(const MPNumber *x, MPNumber *z)
518{
  MPNumber f;
  mp_floor(x, &f);
  mp_subtract(x, &f, z);
522}


525void
526mp_floor(const MPNumber *x, MPNumber *z)
527{
  int i;
  bool_Bool have_fraction = false0, 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_SIZE1000; i++) {
      if (z->fraction[i])
          have_fraction = true1;
      z->fraction[i] = 0;
  }
  z->im_sign = 0;
  z->im_exponent = 0;
  memset(z->im_fraction, 0, sizeof(int) * MP_SIZE1000);

  if (have_fraction && is_negative)
     mp_add_integer(z, -1, z);
558}


561void
562mp_ceiling(const MPNumber *x, MPNumber *z)
563{
  MPNumber f;

  mp_floor(x, z);
  mp_fractional_component(x, &f);
  if (mp_is_zero(&f))
      return;
  mp_add_integer(z, 1, z);
571}


574void
575mp_round(const MPNumber *x, MPNumber *z)
576{
  MPNumber f, one;
  bool_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);
593}

595int
596mp_compare_mp_to_mp(const MPNumber *x, const MPNumber *y)
597{
  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_SIZE1000; 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;
632}


635void
636mp_divide(const MPNumber *x, const MPNumber *y, MPNumber *z)
637{
  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")gettext ("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;
664}


667static void
668mp_divide_integer_real(const MPNumber *x, int64_t y, MPNumber *z)
669{
  int c, i, k, b2, c2, j1, j2;
  MPNumber x_copy;

  /* x/0 */
  if (y == 0) {
      /* Translators: Error displayed attempted to divide by zero */
      mperr(_("Division by zero is undefined")gettext ("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;
  }

  /* 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)((10000 << 3) >= (32767 / 10000) ? (10000 << 3
) : (32767 / 10000));
  if (y < b2) {
      int kh, r1;

      /* LOOK FOR FIRST NONZERO DIGIT IN QUOTIENT */
      do {
          c = MP_BASE10000 * c;
          if (i < MP_T100)
              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_BASE10000 * (c - y * r1);
      kh = 1;
      if (i < MP_T100) {
          kh = MP_T100 + 1 - i;
          for (k = 1; k < kh; k++) {
              c += x_copy.fraction[i];
              z->fraction[k] = c / y;
              c = MP_BASE10000 * (c - y * z->fraction[k]);
              i++;
          }
          if (c < 0)
              goto L210;
      }

      for (k = kh; k < MP_T100 + 4; k++) {
          z->fraction[k] = c / y;
          c = MP_BASE10000 * (c - y * z->fraction[k]);
      }
      if (c < 0)
          goto L210;

      mp_normalize(z);
      return;
  }

  /* HERE NEED SIMULATED DOUBLE-PRECISION DIVISION */
  j1 = y / MP_BASE10000;
  j2 = y - j1 * MP_BASE10000;

  /* LOOK FOR FIRST NONZERO DIGIT */
  c2 = 0;
  do {
      c = MP_BASE10000 * c + c2;
      c2 = i < MP_T100 ? 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_T100 + 4; k++) {
      int ir, iq, iqj;

      /* GET APPROXICAFE 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_BASE10000 - ir * j2;
      if (iq < 0) {
          /* HERE IQ NEGATIVE SO IR WAS TOO LARGE */
          ir--;
          iq += y;
      }

      if (i < MP_T100)
          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);

812L210:
  /* CARRY NEGATIVE SO OVERFLOW MUST HAVE OCCURRED */
  mperr("*** INTEGER OVERFLOW IN MP_DIVIDE_INTEGER, B TOO LARGE ***");
  mp_set_from_integer(0, z);
816}


819void
820mp_divide_integer(const MPNumber *x, int64_t y, MPNumber *z)
821{
  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);
833}


836bool_Bool
837mp_is_integer(const MPNumber *x)
838{
  MPNumber t1, t2, t3;

  if (mp_is_complex(x))
5
←
Taking false branch→
      return false0;

  /* 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);
6
←
Calling 'mp_multiply'→
  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;*/
873}


876bool_Bool
877mp_is_positive_integer(const MPNumber *x)
878{
  if (mp_is_complex(x))
      return false0;
  else
      return x->sign >= 0 && mp_is_integer(x);
883}


886bool_Bool
887mp_is_natural(const MPNumber *x)
888{
  if (mp_is_complex(x))
      return false0;
  else
      return x->sign > 0 && mp_is_integer(x);
893}


896bool_Bool
897mp_is_complex(const MPNumber *x)
898{
  return x->im_sign != 0;
900}


903bool_Bool
904mp_is_equal(const MPNumber *x, const MPNumber *y)
905{
  return mp_compare_mp_to_mp(x, y) == 0;
907}


910/*  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.
*/
917static void
918mp_exp(const MPNumber *x, MPNumber *z)
919{
  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_BASE10000);

  /* Compute approxicafely optimal q (and divide x by 2^q) */
  q = (int)(sqrt((float)MP_T100 * 0.48f * rlb) + (float) x->exponent * 1.44f * rlb);

  /* HALVE Q TIMES */
  if (q > 0) {
      int ib, ic;

      ib = MP_BASE10000 << 2;
      ic = 1;
      for (i = 1; i <= q; ++i) {
          ic *= 2;
          if (ic < ib && ic != MP_BASE10000 && 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_T100 + 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);
981}


984static void
985mp_epowy_real(const MPNumber *x, MPNumber *z)
986{
  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_T100 < 4)
      tss = MP_T100 + 1;
  else
      tss = MP_T100 + 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)((tss) <= (tss + 2 + t1.exponent) ? (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 ***");
1066}


1069void
1070mp_epowy(const MPNumber *x, MPNumber *z)
1071{
  /* 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);
1089}


1092/*  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
*/
1096void
1097mp_gcd(int64_t *k, int64_t *l)
1098{
  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;
1126}


1129bool_Bool
1130mp_is_zero(const MPNumber *x)
1131{
  return x->sign == 0 && x->im_sign == 0;
1133}


1136bool_Bool
1137mp_is_negative(const MPNumber *x)
1138{
  return x->sign < 0;
1140}


1143bool_Bool
1144mp_is_greater_equal(const MPNumber *x, const MPNumber *y)
1145{
  return mp_compare_mp_to_mp(x, y) >= 0;
1147}


1150bool_Bool
1151mp_is_greater_than(const MPNumber *x, const MPNumber *y)
1152{
  return mp_compare_mp_to_mp(x, y) > 0;
1154}


1157bool_Bool
1158mp_is_less_equal(const MPNumber *x, const MPNumber *y)
1159{
  return mp_compare_mp_to_mp(x, y) <= 0;
1161}


1164/*  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.
*/
1169static void
1170mp_lns(const MPNumber *x, MPNumber *z)
1171{
  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_T100)
          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_T100;
      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_T100) {
      mperr("*** ERROR OCCURRED IN MP_LNS, NEWTON ITERATION NOT CONVERGING PROPERLY ***");
  }

  z->sign = -z->sign;
1227}


1230static void
1231mp_ln_real(const MPNumber *x, MPNumber *z)
1232{
  int k;
  MPNumber t1, t2;

  /* LOOP TO GET APPROXICAFE LN(X) USING SINGLE-PRECISION */
  mp_set_from_mp(x, &t1);
  mp_set_from_integer(0, z);
  for(k = 0; k < 10; k++)
  {
      int e;
      double rx, rlx;

      /* 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_BASE10000);
      mp_set_from_double(-(double)rlx, &t2);

      /* UPDATE Z AND COMPUTE ACCURATE EXP OF APPROXICAFE 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 ***");
1269}


1272void
1273mp_ln(const MPNumber *x, MPNumber *z)
1274{
  /* ln(0) undefined */
  if (mp_is_zero(x)) {
      /* Translators: Error displayed when attempting to take logarithm of zero */
      mperr(_("Logarithm of zero is undefined")gettext ("Logarithm of zero is undefined"));
      mp_set_from_integer(0, z);
      return;
  }

  /* ln(-x) complex */
  /* FIXME: Make complex numbers optional */
  /*if (mp_is_negative(x)) {
      // Translators: Error displayed attempted to take logarithm of negative value
      mperr(_("Logarithm of negative values is undefined"));
      mp_set_from_integer(0, z);
      return;
  }*/

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


1307void
1308mp_logarithm(int64_t n, const MPNumber *x, MPNumber *z)
1309{
  MPNumber t1, t2;

  /* log(0) undefined */
  if (mp_is_zero(x)) {
      /* Translators: Error displayed when attempting to take logarithm of zero */
      mperr(_("Logarithm of zero is undefined")gettext ("Logarithm of zero is undefined"));
      mp_set_from_integer(0, z);
      return;
  }

  /* logn(x) = ln(x) / ln(n) */
  mp_set_from_integer(n, &t1);
  mp_ln(&t1, &t1);
  mp_ln(x, &t2);
  mp_divide(&t2, &t1, z);
1325}


1328bool_Bool
1329mp_is_less_than(const MPNumber *x, const MPNumber *y)
1330{
  return mp_compare_mp_to_mp(x, y) < 0;
1332}


1335static void
1336mp_multiply_real(const MPNumber *x, const MPNumber *y, MPNumber *z)
1337{
  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_T100; 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)((100) <= (100 + 3 - i) ? (100) : (100 + 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_BASE10000) {
          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_T100 + 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_BASE10000;
          r.fraction[j] = ri - MP_BASE10000 * 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_BASE10000) {
          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_T100 + 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_BASE10000;
          r.fraction[i] = ri - MP_BASE10000 * 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_SIZE1000);

  /* NORMALIZE AND ROUND RESULT */
  // FIXME: Use stack variable because of mp_normalize brokeness
  for (i = 0; i < MP_SIZE1000; i++)
      z->fraction[i] = r.fraction[i];
  mp_normalize(z);
1433}


1436void
1437mp_multiply(const MPNumber *x, const MPNumber *y, MPNumber *z)
1438{
  /* 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)) {
7
←
Taking true branch→
      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);
8
←
Calling 'mp_subtract'→

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


1470static void
1471mp_multiply_integer_real(const MPNumber *x, int64_t y, MPNumber *z)
1472{
  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)((10000 << 3) >= (32767 / 10000) ? (10000 << 3
) : (32767 / 10000))) {
      int64_t j1, j2;

      /* HERE J IS TOO LARGE FOR SINGLE-PRECISION MULTIPLICATION */
      j1 = y / MP_BASE10000;
      j2 = y - j1 * MP_BASE10000;

      /* FORM PRODUCT */
      for (i = MP_T100 + 3; i >= 0; i--) {
          int64_t c1, c2, is, ix, t;

          c1 = c / MP_BASE10000;
          c2 = c - MP_BASE10000 * c1;
          ix = 0;
          if (i > 3)
              ix = x_copy.fraction[i - 4];

          t = j2 * ix + c2;
          is = t / MP_BASE10000;
          c = j1 * ix + c1 + is;
          z->fraction[i] = t - MP_BASE10000 * is;
      }
  }
  else
  {
      int64_t ri = 0;

      for (i = MP_T100 + 3; i >= 4; i--) {
          ri = y * x_copy.fraction[i - 4] + c;
          c = ri / MP_BASE10000;
          z->fraction[i] = ri - MP_BASE10000 * 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_BASE10000;
          z->fraction[i] = t - MP_BASE10000 * c;
      }
  }

  /* HAVE TO SHIFT RIGHT HERE AS CARRY OFF END */
  while (c != 0) {
      int64_t t;

      for (i = MP_T100 + 3; i >= 1; i--)
          z->fraction[i] = z->fraction[i - 1];
      t = c;
      c = t / MP_BASE10000;
      z->fraction[0] = t - MP_BASE10000 * 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_SIZE1000);
  mp_normalize(z);
1584}


1587void
1588mp_multiply_integer(const MPNumber *x, int64_t y, MPNumber *z)
1589{
  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);
1600}


1603void
1604mp_multiply_fraction(const MPNumber *x, int64_t numerator, int64_t denominator, MPNumber *z)
1605{
  if (denominator == 0) {
      mperr(_("Division by zero is undefined")gettext ("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);
1621}


1624void
1625mp_invert_sign(const MPNumber *x, MPNumber *z)
1626{
  mp_set_from_mp(x, z);
  z->sign = -z->sign;
  z->im_sign = -z->im_sign;
1630}


1633// FIXME: Is r->fraction large enough?  It seems to be in practise but it may be MP_T+4 instead of MP_T
1634// FIXME: There is some sort of stack corruption/use of unitialised variables here.  Some functions are
1635// using stack variables as x otherwise there are corruption errors. e.g. "Cos(45) - 1/Sqrt(2) = -0"
1636// (try in scientific mode)
1637void
1638mp_normalize(MPNumber *x)
1639{
  int start_index;

  /* Find first non-zero digit */
  for (start_index = 0; start_index < MP_SIZE1000 && x->fraction[start_index] == 0; start_index++);

  /* Mark as zero */
  if (start_index >= MP_SIZE1000) {
      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_SIZE1000; i++)
          x->fraction[i] = x->fraction[i + start_index];
      for (; i < MP_SIZE1000; i++)
          x->fraction[i] = 0;
  }
1662}


1665static void
1666mp_pwr(const MPNumber *x, const MPNumber *y, MPNumber *z)
1667{
  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")gettext ("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);
1695}


1698static void
1699mp_reciprocal_real(const MPNumber *x, MPNumber *z)
1700{
  MPNumber t1, t2;
  int it0, t;

  /* 1/0 invalid */
  if (mp_is_zero(x)) {
      mperr(_("Reciprocal of zero is undefined")gettext ("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_T100)
          break;

      /*  FOLLOWING LOOP ALMOST DOUBLES T (POSSIBLE
       *  BECAUSE NEWTONS METHOD HAS 2ND ORDER CONVERGENCE).
       */
      ts3 = t;
      t = MP_T100;
      do {
          ts2 = t;
          t = (t + it0) / 2;
      } while (t > ts3);
      t = min(ts2, MP_T)((ts2) <= (100) ? (ts2) : (100));
  }

  /* RETURN IF NEWTON ITERATION WAS CONVERGING */
  if (t2.sign != 0 && (t1.exponent - t2.exponent) << 1 < MP_T100 - 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);
1750}


1753void
1754mp_reciprocal(const MPNumber *x, MPNumber *z)
1755{
  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);
1773}


1776static void
1777mp_root_real(const MPNumber *x, int64_t n, MPNumber *z)
1778{
  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")gettext ("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)((10000) >= (64) ? (10000) : (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")gettext ("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")gettext ("nth root of negative number is undefined for even n"
));
      mp_set_from_integer(0, z);
      return;
  }

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

      ts3 = t;
      t = MP_T100;
      do {
          ts2 = t;
          t = (t + it0) / 2;
      } while (t > ts3);
      t = min(ts2, MP_T)((ts2) <= (100) ? (ts2) : (100));
  }

  /*  NOW R(I2) IS X**(-1/NP)
   *  CHECK THAT NEWTON ITERATION WAS CONVERGING
   */
  if (t2.sign != 0 && (t1.exponent - t2.exponent) << 1 < MP_T100 - 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);
1878}


1881void
1882mp_root(const MPNumber *x, int64_t n, MPNumber *z)
1883{
  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);
  }
  else
      mp_root_real(x, n, z);
1901}


1904void
1905mp_sqrt(const MPNumber *x, MPNumber *z)
1906{
  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);
  }
1921}


1924void
1925mp_factorial(const MPNumber *x, MPNumber *z)
1926{
  int i, value;

  /* 0! == 1 */
  if (mp_is_zero(x)) {
      mp_set_from_integer(1, z);
      return;
  }
  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")gettext ("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);
1946}


1949void
1950mp_modulus_divide(const MPNumber *x, const MPNumber *y, MPNumber *z)
1951{
  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 */
      mperr(_("Modulus division is only defined for integers")gettext ("Modulus division is only defined for integers"));
      mp_set_from_integer(0, z);
  }

  mp_divide(x, y, &t1);
  mp_floor(&t1, &t1);
  mp_multiply(&t1, y, &t2);
  mp_subtract(x, &t2, z);

  mp_set_from_integer(0, &t1);
  if ((mp_is_less_than(y, &t1) && mp_is_greater_than(z, &t1)) || mp_is_less_than(z, &t1))
      mp_add(z, y, z);
1968}


1971void
1972mp_xpowy(const MPNumber *x, const MPNumber *y, MPNumber *z)
1973{
  if (mp_is_integer(y)) {
      mp_xpowy_integer(x, mp_cast_to_int(y), z);
  } else {
      MPNumber reciprocal;
      mp_reciprocal(y, &reciprocal);
      if (mp_is_integer(&reciprocal))
          mp_root(x, mp_cast_to_int(&reciprocal), z);
      else
          mp_pwr(x, y, z);
  }
1984}


1987void
1988mp_xpowy_integer(const MPNumber *x, int64_t n, MPNumber *z)
1989{
  int i;
  MPNumber t;

  /* 0^-n invalid */
  if (mp_is_zero(x) && n < 0) {
      /* Translators: Error displayed when attempted to raise 0 to a negative exponent */
      mperr(_("The power of zero is undefined for a negative exponent")gettext ("The power of zero is undefined for a negative exponent"
));
      mp_set_from_integer(0, z);
      return;
  }

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


2034GList*
2035mp_factorize(const MPNumber *x)
2036{
  GList *list = NULL((void*)0);
  MPNumber *factor = g_slice_alloc0(sizeof(MPNumber));
  MPNumber value, tmp, divisor, root;

  mp_abs(x, &value);

  if (mp_is_zero(&value)) {
1
Taking false branch→
      mp_set_from_mp(&value, factor);
      list = g_list_append(list, factor);
      return list;
  }

  mp_set_from_integer(1, &tmp);
  if (mp_is_equal(&value, &tmp)) {
2
←
Taking false branch→
      mp_set_from_mp(x, factor);
      list = g_list_append(list, factor);
      return list;
  }

  mp_set_from_integer(2, &divisor);
  while (TRUE(!(0))) {
3
←
Loop condition is true.  Entering loop body→
      mp_divide(&value, &divisor, &tmp);
      if (mp_is_integer(&tmp)) {
4
←
Calling 'mp_is_integer'→
          value = tmp;
          mp_set_from_mp(&divisor, factor);
          list = g_list_append(list, factor);
          factor = g_slice_alloc0(sizeof(MPNumber));
      } else {
          break;
      }
  }

  mp_set_from_integer(3, &divisor);
  mp_sqrt(&value, &root);
  while (mp_is_less_equal(&divisor, &root)) {
      mp_divide(&value, &divisor, &tmp);
      if (mp_is_integer(&tmp)) {
          value = tmp;
          mp_sqrt(&value, &root);
          mp_set_from_mp(&divisor, factor);
          list = g_list_append(list, factor);
          factor = g_slice_alloc0(sizeof(MPNumber));
      } else {
          mp_add_integer(&divisor, 2, &tmp);
          divisor = tmp;
      }
  }

  mp_set_from_integer(1, &tmp);
  if (mp_is_greater_than(&value, &tmp)) {
      mp_set_from_mp(&value, factor);
      list = g_list_append(list, factor);
  } else {
      g_slice_free(MPNumber, factor)do { if (1) g_slice_free1 (sizeof (MPNumber), (factor)); else
 (void) ((MPNumber*) 0 == (factor)); } while (0);
  }

  if (mp_is_negative(x)) {
      mp_invert_sign(list->data, list->data);
  }

  return list;
2098}