src/mp.c

   1 /*  Copyright (c) 1987-2008 Sun Microsystems, Inc. All Rights Reserved.
   2  *  Copyright (c) 2008-2009 Robert Ancell
   3  *
   4  *  This program is free software; you can redistribute it and/or modify
   5  *  it under the terms of the GNU General Public License as published by
   6  *  the Free Software Foundation; either version 2, or (at your option)
   7  *  any later version.
   8  *
   9  *  This program is distributed in the hope that it will be useful, but
  10  *  WITHOUT ANY WARRANTY; without even the implied warranty of
  11  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  12  *  General Public License for more details.
  13  *
  14  *  You should have received a copy of the GNU General Public License
  15  *  along with this program; if not, write to the Free Software
  16  *  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
  17  *  02111-1307, USA.
  18  */
  19
  20 #include <stdlib.h>
  21 #include <stdio.h>
  22 #include <math.h>
  23 #include <errno.h>
  24
  25 #include "mp.h"
  26 #include "mp-private.h"
  27
  28 static MPNumber eulers_number;
  29 static gboolean have_eulers_number = FALSE;
  30
  31 // FIXME: Re-add overflow and underflow detection
  32
  33 char *mp_error = NULL;
  34
  35 /*  THIS ROUTINE IS CALLED WHEN AN ERROR CONDITION IS ENCOUNTERED, AND
  36  *  AFTER A MESSAGE HAS BEEN WRITTEN TO STDERR.
  37  */
  38 void
  39 mperr(const char *format, ...)
  40 {
  41     char text[1024];
  42     va_list args;
  43
  44     va_start(args, format);
  45     vsnprintf(text, 1024, format, args);
  46     va_end(args);
  47
  48     if (mp_error)
  49         free(mp_error);
  50     mp_error = strdup(text);
  51 }
  52
  53
  54 const char *
  55 mp_get_error()
  56 {
  57     return mp_error;
  58 }
  59
  60
  61 void mp_clear_error()
  62 {
  63     if (mp_error)
  64         free(mp_error);
  65     mp_error = NULL;
  66 }
  67
  68
  69 /*  ROUTINE CALLED BY MP_DIVIDE AND MP_SQRT TO ENSURE THAT
  70  *  RESULTS ARE REPRESENTED EXACTLY IN T-2 DIGITS IF THEY
  71  *  CAN BE.  X IS AN MP NUMBER, I AND J ARE INTEGERS.
  72  */
  73 static void
  74 mp_ext(int i, int j, MPNumber *x)
  75 {
  76     int q, s;
  77
  78     if (mp_is_zero(x) || MP_T <= 2 || i == 0)
  79         return;
  80
  81     /* COMPUTE MAXIMUM POSSIBLE ERROR IN THE LAST PLACE */
  82     q = (j + 1) / i + 1;
  83     s = MP_BASE * x->fraction[MP_T - 2] + x->fraction[MP_T - 1];
  84
  85     /* SET LAST TWO DIGITS TO ZERO */
  86     if (s <= q) {
  87         x->fraction[MP_T - 2] = 0;
  88         x->fraction[MP_T - 1] = 0;
  89         return;
  90     }
  91
  92     if (s + q < MP_BASE * MP_BASE)
  93         return;
  94
  95     /* ROUND UP HERE */
  96     x->fraction[MP_T - 2] = MP_BASE - 1;
  97     x->fraction[MP_T - 1] = MP_BASE;
  98
  99     /* NORMALIZE X (LAST DIGIT B IS OK IN MP_MULTIPLY_INTEGER) */
 100     mp_multiply_integer(x, 1, x);
 101 }
 102
 103
 104 void
 105 mp_get_eulers(MPNumber *z)
 106 {
 107     if (!have_eulers_number) {
 108         MPNumber t;
 109         mp_set_from_integer(1, &t);
 110         mp_epowy(&t, &eulers_number);
 111         have_eulers_number = TRUE;
 112     }
 113     mp_set_from_mp(&eulers_number, z);
 114 }
 115
 116
 117 void
 118 mp_get_i(MPNumber *z)
 119 {
 120     mp_set_from_integer(0, z);
 121     z->im_sign = 1;
 122     z->im_exponent = 1;
 123     z->im_fraction[0] = 1;
 124 }
 125
 126
 127 void
 128 mp_abs(const MPNumber *x, MPNumber *z)
 129 {
 130     if (mp_is_complex(x)){
 131         MPNumber x_real, x_im;
 132
 133         mp_real_component(x, &x_real);
 134         mp_imaginary_component(x, &x_im);
 135
 136         mp_multiply(&x_real, &x_real, &x_real);
 137         mp_multiply(&x_im, &x_im, &x_im);
 138         mp_add(&x_real, &x_im, z);
 139         mp_root(z, 2, z);
 140     }
 141     else {
 142         mp_set_from_mp(x, z);
 143         if (z->sign < 0)
 144             z->sign = -z->sign;
 145     }
 146 }
 147
 148
 149 void
 150 mp_arg(const MPNumber *x, MPAngleUnit unit, MPNumber *z)
 151 {
 152     MPNumber x_real, x_im, pi;
 153
 154     if (mp_is_zero(x)) {
 155         /* Translators: Error display when attempting to take argument of zero */
 156         mperr(_("Argument not defined for zero"));
 157         mp_set_from_integer(0, z);
 158         return;
 159     }
 160
 161     mp_real_component(x, &x_real);
 162     mp_imaginary_component(x, &x_im);
 163     mp_get_pi(&pi);
 164
 165     if (mp_is_zero(&x_im)) {
 166         if (mp_is_negative(&x_real))
 167             convert_from_radians(&pi, MP_RADIANS, z);
 168         else
 169             mp_set_from_integer(0, z);
 170     }
 171     else if (mp_is_zero(&x_real)) {
 172         mp_set_from_mp(&pi, z);
 173         if (mp_is_negative(&x_im))
 174             mp_divide_integer(z, -2, z);
 175         else
 176             mp_divide_integer(z, 2, z);
 177     }
 178     else if (mp_is_negative(&x_real)) {
 179         mp_divide(&x_im, &x_real, z);
 180         mp_atan(z, MP_RADIANS, z);
 181         if (mp_is_negative(&x_im))
 182             mp_subtract(z, &pi, z);
 183         else
 184             mp_add(z, &pi, z);
 185     }
 186     else {
 187         mp_divide(&x_im, &x_real, z);
 188         mp_atan(z, MP_RADIANS, z);
 189     }
 190
 191     convert_from_radians(z, unit, z);
 192 }
 193
 194
 195 void
 196 mp_conjugate(const MPNumber *x, MPNumber *z)
 197 {
 198     mp_set_from_mp(x, z);
 199     z->im_sign = -z->im_sign;
 200 }
 201
 202
 203 void
 204 mp_real_component(const MPNumber *x, MPNumber *z)
 205 {
 206     mp_set_from_mp(x, z);
 207
 208     /* Clear imaginary component */
 209     z->im_sign = 0;
 210     z->im_exponent = 0;
 211     memset(z->im_fraction, 0, sizeof(int) * MP_SIZE);
 212 }
 213
 214
 215 void
 216 mp_imaginary_component(const MPNumber *x, MPNumber *z)
 217 {
 218     /* Copy imaginary component to real component */
 219     z->sign = x->im_sign;
 220     z->exponent = x->im_exponent;
 221     memcpy(z->fraction, x->im_fraction, sizeof(int) * MP_SIZE);
 222
 223     /* Clear (old) imaginary component */
 224     z->im_sign = 0;
 225     z->im_exponent = 0;
 226     memset(z->im_fraction, 0, sizeof(int) * MP_SIZE);
 227 }
 228
 229
 230 static void
 231 mp_add_real(const MPNumber *x, int y_sign, const MPNumber *y, MPNumber *z)
 232 {
 233     int sign_prod, i, c;
 234     int exp_diff, med;
 235     bool x_largest = false;
 236     const int *big_fraction, *small_fraction;
 237     MPNumber x_copy, y_copy;
 238
 239     /* 0 + y = y */
 240     if (mp_is_zero(x)) {
 241         mp_set_from_mp(y, z);
 242         z->sign = y_sign;
 243         return;
 244     }
 245     /* x + 0 = x */
 246     else if (mp_is_zero(y)) {
 247         mp_set_from_mp(x, z);
 248         return;
 249     }
 250
 251     sign_prod = y_sign * x->sign;
 252     exp_diff = x->exponent - y->exponent;
 253     med = abs(exp_diff);
 254     if (exp_diff < 0) {
 255         x_largest = false;
 256     } else if (exp_diff > 0) {
 257         x_largest = true;
 258     } else {
 259         /* EXPONENTS EQUAL SO COMPARE SIGNS, THEN FRACTIONS IF NEC. */
 260         if (sign_prod < 0) {
 261             /* Signs are not equal.  find out which mantissa is larger. */
 262             int j;
 263             for (j = 0; j < MP_T; j++) {
 264                 int i = x->fraction[j] - y->fraction[j];
 265                 if (i == 0)
 266                     continue;
 267                 if (i < 0)
 268                     x_largest = false;
 269                 else if (i > 0)
 270                     x_largest = true;
 271                 break;
 272             }
 273
 274             /* Both mantissas equal, so result is zero. */
 275             if (j >= MP_T) {
 276                 mp_set_from_integer(0, z);
 277                 return;
 278             }
 279         }
 280     }
 281
 282     mp_set_from_mp(x, &x_copy);
 283     mp_set_from_mp(y, &y_copy);
 284     mp_set_from_integer(0, z);
 285
 286     if (x_largest) {
 287         z->sign = x_copy.sign;
 288         z->exponent = x_copy.exponent;
 289         big_fraction = x_copy.fraction;
 290         small_fraction = y_copy.fraction;
 291     } else {
 292         z->sign = y_sign;
 293         z->exponent = y_copy.exponent;
 294         big_fraction = y_copy.fraction;
 295         small_fraction = x_copy.fraction;
 296     }
 297
 298     /* CLEAR GUARD DIGITS TO RIGHT OF X DIGITS */
 299     for(i = 3; i >= med; i--)
 300         z->fraction[MP_T + i] = 0;
 301
 302     if (sign_prod >= 0) {
 303         /* This is probably insufficient overflow detection, but it makes us
 304          * not crash at least.
 305          */
 306         if (MP_T + 3 < med) {
 307             mperr(_("Overflow: the result couldn't be calculated"));
 308             mp_set_from_integer(0, z);
 309             return;
 310         }
 311
 312         /* HERE DO ADDITION, EXPONENT(Y) >= EXPONENT(X) */
 313         for (i = MP_T + 3; i >= MP_T; i--)
 314             z->fraction[i] = small_fraction[i - med];
 315
 316         c = 0;
 317         for (; i >= med; i--) {
 318             c = big_fraction[i] + small_fraction[i - med] + c;
 319
 320             if (c < MP_BASE) {
 321                 /* NO CARRY GENERATED HERE */
 322                 z->fraction[i] = c;
 323                 c = 0;
 324             } else {
 325                 /* CARRY GENERATED HERE */
 326                 z->fraction[i] = c - MP_BASE;
 327                 c = 1;
 328             }
 329         }
 330
 331         for (; i >= 0; i--)
 332         {
 333             c = big_fraction[i] + c;
 334             if (c < MP_BASE) {
 335                 z->fraction[i] = c;
 336                 i--;
 337
 338                 /* NO CARRY POSSIBLE HERE */
 339                 for (; i >= 0; i--)
 340                     z->fraction[i] = big_fraction[i];
 341
 342                 c = 0;
 343                 break;
 344             }
 345
 346             z->fraction[i] = 0;
 347             c = 1;
 348         }
 349
 350         /* MUST SHIFT RIGHT HERE AS CARRY OFF END */
 351         if (c != 0) {
 352             for (i = MP_T + 3; i > 0; i--)
 353                 z->fraction[i] = z->fraction[i - 1];
 354             z->fraction[0] = 1;
 355             z->exponent++;
 356         }
 357     }
 358     else  {
 359         c = 0;
 360         for (i = MP_T + med - 1; i >= MP_T; i--) {
 361             /* HERE DO SUBTRACTION, ABS(Y) > ABS(X) */
 362             z->fraction[i] = c - small_fraction[i - med];
 363             c = 0;
 364
 365             /* BORROW GENERATED HERE */
 366             if (z->fraction[i] < 0) {
 367                 c = -1;
 368                 z->fraction[i] += MP_BASE;
 369             }
 370         }
 371
 372         for(; i >= med; i--) {
 373             c = big_fraction[i] + c - small_fraction[i - med];
 374             if (c >= 0) {
 375                 /* NO BORROW GENERATED HERE */
 376                 z->fraction[i] = c;
 377                 c = 0;
 378             } else {
 379                 /* BORROW GENERATED HERE */
 380                 z->fraction[i] = c + MP_BASE;
 381                 c = -1;
 382             }
 383         }
 384
 385         for (; i >= 0; i--) {
 386             c = big_fraction[i] + c;
 387
 388             if (c >= 0) {
 389                 z->fraction[i] = c;
 390                 i--;
 391
 392                 /* NO CARRY POSSIBLE HERE */
 393                 for (; i >= 0; i--)
 394                     z->fraction[i] = big_fraction[i];
 395
 396                 break;
 397             }
 398
 399             z->fraction[i] = c + MP_BASE;
 400             c = -1;
 401         }
 402     }
 403
 404     mp_normalize(z);
 405 }
 406
 407
 408 static void
 409 mp_add_with_sign(const MPNumber *x, int y_sign, const MPNumber *y, MPNumber *z)
 410 {
 411     if (mp_is_complex(x) || mp_is_complex(y)) {
 412         MPNumber real_x, real_y, im_x, im_y, real_z, im_z;
 413
 414         mp_real_component(x, &real_x);
 415         mp_imaginary_component(x, &im_x);
 416         mp_real_component(y, &real_y);
 417         mp_imaginary_component(y, &im_y);
 418
 419         mp_add_real(&real_x, y_sign * y->sign, &real_y, &real_z);
 420         mp_add_real(&im_x, y_sign * y->im_sign, &im_y, &im_z);
 421
 422         mp_set_from_complex(&real_z, &im_z, z);
 423     }
 424     else
 425         mp_add_real(x, y_sign * y->sign, y, z);
 426 }
 427
 428
 429 void
 430 mp_add(const MPNumber *x, const MPNumber *y, MPNumber *z)
 431 {
 432     mp_add_with_sign(x, 1, y, z);
 433 }
 434
 435
 436 void
 437 mp_add_integer(const MPNumber *x, int64_t y, MPNumber *z)
 438 {
 439     MPNumber t;
 440     mp_set_from_integer(y, &t);
 441     mp_add(x, &t, z);
 442 }
 443
 444
 445 void
 446 mp_add_fraction(const MPNumber *x, int64_t i, int64_t j, MPNumber *y)
 447 {
 448     MPNumber t;
 449     mp_set_from_fraction(i, j, &t);
 450     mp_add(x, &t, y);
 451 }
 452
 453
 454 void
 455 mp_subtract(const MPNumber *x, const MPNumber *y, MPNumber *z)
 456 {
 457     mp_add_with_sign(x, -1, y, z);
 458 }
 459
 460
 461 void
 462 mp_sgn(const MPNumber *x, MPNumber *z)
 463 {
 464     if (mp_is_zero(x))
 465         mp_set_from_integer(0, z);
 466     else if (mp_is_negative(x))
 467         mp_set_from_integer(-1, z);
 468     else
 469         mp_set_from_integer(1, z);
 470 }
 471
 472 void
 473 mp_integer_component(const MPNumber *x, MPNumber *z)
 474 {
 475     int i;
 476
 477     /* Clear fraction */
 478     mp_set_from_mp(x, z);
 479     for (i = z->exponent; i < MP_SIZE; i++)
 480         z->fraction[i] = 0;
 481     z->im_sign = 0;
 482     z->im_exponent = 0;
 483     memset(z->im_fraction, 0, sizeof(int) * MP_SIZE);
 484 }
 485
 486 void
 487 mp_fractional_component(const MPNumber *x, MPNumber *z)
 488 {
 489     int i, shift;
 490
 491     /* Fractional component of zero is 0 */
 492     if (mp_is_zero(x)) {
 493         mp_set_from_integer(0, z);
 494         return;
 495     }
 496
 497     /* All fractional */
 498     if (x->exponent <= 0) {
 499         mp_set_from_mp(x, z);
 500         return;
 501     }
 502
 503     /* Shift fractional component */
 504     shift = x->exponent;
 505     for (i = shift; i < MP_SIZE && x->fraction[i] == 0; i++)
 506         shift++;
 507     z->sign = x->sign;
 508     z->exponent = x->exponent - shift;
 509     for (i = 0; i < MP_SIZE; i++) {
 510         if (i + shift >= MP_SIZE)
 511             z->fraction[i] = 0;
 512         else
 513             z->fraction[i] = x->fraction[i + shift];
 514     }
 515     if (z->fraction[0] == 0)
 516         z->sign = 0;
 517
 518     z->im_sign = 0;
 519     z->im_exponent = 0;
 520     memset(z->im_fraction, 0, sizeof(int) * MP_SIZE);
 521 }
 522
 523
 524 void
 525 mp_fractional_part(const MPNumber *x, MPNumber *z)
 526 {
 527     MPNumber f;
 528     mp_floor(x, &f);
 529     mp_subtract(x, &f, z);
 530 }
 531
 532
 533 void
 534 mp_floor(const MPNumber *x, MPNumber *z)
 535 {
 536     int i;
 537     bool have_fraction = false, is_negative;
 538
 539     /* Integer component of zero = 0 */
 540     if (mp_is_zero(x)) {
 541         mp_set_from_mp(x, z);
 542         return;
 543     }
 544
 545     /* If all fractional then no integer component */
 546     if (x->exponent <= 0) {
 547         mp_set_from_integer(0, z);
 548         return;
 549     }
 550
 551     is_negative = mp_is_negative(x);
 552
 553     /* Clear fraction */
 554     mp_set_from_mp(x, z);
 555     for (i = z->exponent; i < MP_SIZE; i++) {
 556         if (z->fraction[i])
 557             have_fraction = true;
 558         z->fraction[i] = 0;
 559     }
 560     z->im_sign = 0;
 561     z->im_exponent = 0;
 562     memset(z->im_fraction, 0, sizeof(int) * MP_SIZE);
 563
 564     if (have_fraction && is_negative)
 565        mp_add_integer(z, -1, z);
 566 }
 567
 568
 569 void
 570 mp_ceiling(const MPNumber *x, MPNumber *z)
 571 {
 572     MPNumber f;
 573
 574     mp_floor(x, z);
 575     mp_fractional_component(x, &f);
 576     if (mp_is_zero(&f))
 577         return;
 578     mp_add_integer(z, 1, z);
 579 }
 580
 581
 582 void
 583 mp_round(const MPNumber *x, MPNumber *z)
 584 {
 585     MPNumber f, one;
 586     bool do_floor;
 587
 588     do_floor = !mp_is_negative(x);
 589
 590     mp_fractional_component(x, &f);
 591     mp_multiply_integer(&f, 2, &f);
 592     mp_abs(&f, &f);
 593     mp_set_from_integer(1, &one);
 594     if (mp_is_greater_equal(&f, &one))
 595         do_floor = !do_floor;
 596
 597     if (do_floor)
 598         mp_floor(x, z);
 599     else
 600         mp_ceiling(x, z);
 601 }
 602
 603 int
 604 mp_compare_mp_to_mp(const MPNumber *x, const MPNumber *y)
 605 {
 606     int i;
 607
 608     if (x->sign != y->sign) {
 609         if (x->sign > y->sign)
 610             return 1;
 611         else
 612             return -1;
 613     }
 614
 615     /* x = y = 0 */
 616     if (mp_is_zero(x))
 617         return 0;
 618
 619     /* See if numbers are of different magnitude */
 620     if (x->exponent != y->exponent) {
 621         if (x->exponent > y->exponent)
 622             return x->sign;
 623         else
 624             return -x->sign;
 625     }
 626
 627     /* Compare fractions */
 628     for (i = 0; i < MP_SIZE; i++) {
 629         if (x->fraction[i] == y->fraction[i])
 630             continue;
 631
 632         if (x->fraction[i] > y->fraction[i])
 633             return x->sign;
 634         else
 635             return -x->sign;
 636     }
 637
 638     /* x = y */
 639     return 0;
 640 }
 641
 642
 643 void
 644 mp_divide(const MPNumber *x, const MPNumber *y, MPNumber *z)
 645 {
 646     int i, ie;
 647     MPNumber t;
 648
 649     /* x/0 */
 650     if (mp_is_zero(y)) {
 651         /* Translators: Error displayed attempted to divide by zero */
 652         mperr(_("Division by zero is undefined"));
 653         mp_set_from_integer(0, z);
 654         return;
 655     }
 656
 657     /* 0/y = 0 */
 658     if (mp_is_zero(x)) {
 659         mp_set_from_integer(0, z);
 660         return;
 661     }
 662
 663     /* z = x × y⁻¹ */
 664     /* FIXME: Set exponent to zero to avoid overflow in mp_multiply??? */
 665     mp_reciprocal(y, &t);
 666     ie = t.exponent;
 667     t.exponent = 0;
 668     i = t.fraction[0];
 669     mp_multiply(x, &t, z);
 670     mp_ext(i, z->fraction[0], z);
 671     z->exponent += ie;
 672 }
 673
 674
 675 static void
 676 mp_divide_integer_real(const MPNumber *x, int64_t y, MPNumber *z)
 677 {
 678     int c, i, k, b2, c2, j1, j2;
 679     MPNumber x_copy;
 680
 681     /* x/0 */
 682     if (y == 0) {
 683         /* Translators: Error displayed attempted to divide by zero */
 684         mperr(_("Division by zero is undefined"));
 685         mp_set_from_integer(0, z);
 686         return;
 687     }
 688
 689     /* 0/y = 0 */
 690     if (mp_is_zero(x)) {
 691         mp_set_from_integer(0, z);
 692         return;
 693     }
 694
 695     /* Division by -1 or 1 just changes sign */
 696     if (y == 1 || y == -1) {
 697         if (y < 0)
 698             mp_invert_sign(x, z);
 699         else
 700             mp_set_from_mp(x, z);
 701         return;
 702     }
 703
 704     /* Copy x as z may also refer to x */
 705     mp_set_from_mp(x, &x_copy);
 706     mp_set_from_integer(0, z);
 707
 708     if (y < 0) {
 709         y = -y;
 710         z->sign = -x_copy.sign;
 711     }
 712     else
 713         z->sign = x_copy.sign;
 714     z->exponent = x_copy.exponent;
 715
 716     c = 0;
 717     i = 0;
 718
 719     /*  IF y*B NOT REPRESENTABLE AS AN INTEGER HAVE TO SIMULATE
 720      *  LONG DIVISION.  ASSUME AT LEAST 16-BIT WORD.
 721      */
 722
 723     /* Computing MAX */
 724     b2 = max(MP_BASE << 3, 32767 / MP_BASE);
 725     if (y < b2) {
 726         int kh, r1;
 727
 728         /* LOOK FOR FIRST NONZERO DIGIT IN QUOTIENT */
 729         do {
 730             c = MP_BASE * c;
 731             if (i < MP_T)
 732                 c += x_copy.fraction[i];
 733             i++;
 734             r1 = c / y;
 735             if (r1 < 0)
 736                 goto L210;
 737         } while (r1 == 0);
 738
 739         /* ADJUST EXPONENT AND GET T+4 DIGITS IN QUOTIENT */
 740         z->exponent += 1 - i;
 741         z->fraction[0] = r1;
 742         c = MP_BASE * (c - y * r1);
 743         kh = 1;
 744         if (i < MP_T) {
 745             kh = MP_T + 1 - i;
 746             for (k = 1; k < kh; k++) {
 747                 c += x_copy.fraction[i];
 748                 z->fraction[k] = c / y;
 749                 c = MP_BASE * (c - y * z->fraction[k]);
 750                 i++;
 751             }
 752             if (c < 0)
 753                 goto L210;
 754         }
 755
 756         for (k = kh; k < MP_T + 4; k++) {
 757             z->fraction[k] = c / y;
 758             c = MP_BASE * (c - y * z->fraction[k]);
 759         }
 760         if (c < 0)
 761             goto L210;
 762
 763         mp_normalize(z);
 764         return;
 765     }
 766
 767     /* HERE NEED SIMULATED DOUBLE-PRECISION DIVISION */
 768     j1 = y / MP_BASE;
 769     j2 = y - j1 * MP_BASE;
 770
 771     /* LOOK FOR FIRST NONZERO DIGIT */
 772     c2 = 0;
 773     do {
 774         c = MP_BASE * c + c2;
 775         c2 = i < MP_T ? x_copy.fraction[i] : 0;
 776         i++;
 777     } while (c < j1 || (c == j1 && c2 < j2));
 778
 779     /* COMPUTE T+4 QUOTIENT DIGITS */
 780     z->exponent += 1 - i;
 781     i--;
 782
 783     /* MAIN LOOP FOR LARGE ABS(y) CASE */
 784     for (k = 1; k <= MP_T + 4; k++) {
 785         int ir, iq, iqj;
 786
 787         /* GET APPROXIMATE QUOTIENT FIRST */
 788         ir = c / (j1 + 1);
 789
 790         /* NOW REDUCE SO OVERFLOW DOES NOT OCCUR */
 791         iq = c - ir * j1;
 792         if (iq >= b2) {
 793             /* HERE IQ*B WOULD POSSIBLY OVERFLOW SO INCREASE IR */
 794             ++ir;
 795             iq -= j1;
 796         }
 797
 798         iq = iq * MP_BASE - ir * j2;
 799         if (iq < 0) {
 800             /* HERE IQ NEGATIVE SO IR WAS TOO LARGE */
 801             ir--;
 802             iq += y;
 803         }
 804
 805         if (i < MP_T)
 806             iq += x_copy.fraction[i];
 807         i++;
 808         iqj = iq / y;
 809
 810         /* R(K) = QUOTIENT, C = REMAINDER */
 811         z->fraction[k - 1] = iqj + ir;
 812         c = iq - y * iqj;
 813
 814         if (c < 0)
 815             goto L210;
 816     }
 817
 818     mp_normalize(z);
 819
 820 L210:
 821     /* CARRY NEGATIVE SO OVERFLOW MUST HAVE OCCURRED */
 822     mperr("*** INTEGER OVERFLOW IN MP_DIVIDE_INTEGER, B TOO LARGE ***");
 823     mp_set_from_integer(0, z);
 824 }
 825
 826
 827 void
 828 mp_divide_integer(const MPNumber *x, int64_t y, MPNumber *z)
 829 {
 830     if (mp_is_complex(x)) {
 831         MPNumber re_z, im_z;
 832
 833         mp_real_component(x, &re_z);
 834         mp_imaginary_component(x, &im_z);
 835         mp_divide_integer_real(&re_z, y, &re_z);
 836         mp_divide_integer_real(&im_z, y, &im_z);
 837         mp_set_from_complex(&re_z, &im_z, z);
 838     }
 839     else
 840         mp_divide_integer_real(x, y, z);
 841 }
 842
 843
 844 bool
 845 mp_is_integer(const MPNumber *x)
 846 {
 847     MPNumber t1, t2, t3;
 848
 849     if (mp_is_complex(x))
 850         return false;
 851
 852     /* This fix is required for 1/3 repiprocal not being detected as an integer */
 853     /* Multiplication and division by 10000 is used to get around a
 854      * limitation to the "fix" for Sun bugtraq bug #4006391 in the
 855      * mp_floor() routine in mp.c, when the exponent is less than 1.
 856      */
 857     mp_set_from_integer(10000, &t3);
 858     mp_multiply(x, &t3, &t1);
 859     mp_divide(&t1, &t3, &t1);
 860     mp_floor(&t1, &t2);
 861     return mp_is_equal(&t1, &t2);
 862
 863     /* Correct way to check for integer */
 864     /*int i;
 865
 866     // Zero is an integer
 867     if (mp_is_zero(x))
 868         return true;
 869
 870     // Fractional
 871     if (x->exponent <= 0)
 872         return false;
 873
 874     // Look for fractional components
 875     for (i = x->exponent; i < MP_SIZE; i++) {
 876         if (x->fraction[i] != 0)
 877             return false;
 878     }
 879
 880     return true;*/
 881 }
 882
 883
 884 bool
 885 mp_is_positive_integer(const MPNumber *x)
 886 {
 887     if (mp_is_complex(x))
 888         return false;
 889     else
 890         return x->sign >= 0 && mp_is_integer(x);
 891 }
 892
 893
 894 bool
 895 mp_is_natural(const MPNumber *x)
 896 {
 897     if (mp_is_complex(x))
 898         return false;
 899     else
 900         return x->sign > 0 && mp_is_integer(x);
 901 }
 902
 903
 904 bool
 905 mp_is_complex(const MPNumber *x)
 906 {
 907     return x->im_sign != 0;
 908 }
 909
 910
 911 bool
 912 mp_is_equal(const MPNumber *x, const MPNumber *y)
 913 {
 914     return mp_compare_mp_to_mp(x, y) == 0;
 915 }
 916
 917
 918 /*  Return e^x for |x| < 1 USING AN O(SQRT(T).M(T)) ALGORITHM
 919  *  DESCRIBED IN - R. P. BRENT, THE COMPLEXITY OF MULTIPLE-
 920  *  PRECISION ARITHMETIC (IN COMPLEXITY OF COMPUTATIONAL PROBLEM
 921  *  SOLVING, UNIV. OF QUEENSLAND PRESS, BRISBANE, 1976, 126-165).
 922  *  ASYMPTOTICALLY FASTER METHODS EXIST, BUT ARE NOT USEFUL
 923  *  UNLESS T IS VERY LARGE. SEE COMMENTS TO MP_ATAN AND MPPIGL.
 924  */
 925 static void
 926 mp_exp(const MPNumber *x, MPNumber *z)
 927 {
 928     int i, q;
 929     float rlb;
 930     MPNumber t1, t2;
 931
 932     /* e^0 = 1 */
 933     if (mp_is_zero(x)) {
 934         mp_set_from_integer(1, z);
 935         return;
 936     }
 937
 938     /* Only defined for |x| < 1 */
 939     if (x->exponent > 0) {
 940         mperr("*** ABS(X) NOT LESS THAN 1 IN CALL TO MP_EXP ***");
 941         mp_set_from_integer(0, z);
 942         return;
 943     }
 944
 945     mp_set_from_mp(x, &t1);
 946     rlb = log((float)MP_BASE);
 947
 948     /* Compute approximately optimal q (and divide x by 2^q) */
 949     q = (int)(sqrt((float)MP_T * 0.48f * rlb) + (float) x->exponent * 1.44f * rlb);
 950
 951     /* HALVE Q TIMES */
 952     if (q > 0) {
 953         int ib, ic;
 954
 955         ib = MP_BASE << 2;
 956         ic = 1;
 957         for (i = 1; i <= q; ++i) {
 958             ic *= 2;
 959             if (ic < ib && ic != MP_BASE && i < q)
 960                 continue;
 961             mp_divide_integer(&t1, ic, &t1);
 962             ic = 1;
 963         }
 964     }
 965
 966     if (mp_is_zero(&t1)) {
 967         mp_set_from_integer(0, z);
 968         return;
 969     }
 970
 971     /* Sum series, reducing t where possible */
 972     mp_set_from_mp(&t1, z);
 973     mp_set_from_mp(&t1, &t2);
 974     for (i = 2; MP_T + t2.exponent - z->exponent > 0; i++) {
 975         mp_multiply(&t1, &t2, &t2);
 976         mp_divide_integer(&t2, i, &t2);
 977         mp_add(&t2, z, z);
 978         if (mp_is_zero(&t2))
 979             break;
 980     }
 981
 982     /* Apply (x+1)^2 - 1 = x(2 + x) for q iterations */
 983     for (i = 1; i <= q; ++i) {
 984         mp_add_integer(z, 2, &t1);
 985         mp_multiply(&t1, z, z);
 986     }
 987
 988     mp_add_integer(z, 1, z);
 989 }
 990
 991
 992 static void
 993 mp_epowy_real(const MPNumber *x, MPNumber *z)
 994 {
 995     float r__1;
 996     int i, ix, xs, tss;
 997     float rx, rz, rlb;
 998     MPNumber t1, t2;
 999
1000     /* e^0 = 1 */
1001     if (mp_is_zero(x)) {
1002         mp_set_from_integer(1, z);
1003         return;
1004     }
1005
1006     /* If |x| < 1 use mp_exp */
1007     if (x->exponent <= 0) {
1008         mp_exp(x, z);
1009         return;
1010     }
1011
1012     /*  SEE IF ABS(X) SO LARGE THAT EXP(X) WILL CERTAINLY OVERFLOW
1013      *  OR UNDERFLOW.  1.01 IS TO ALLOW FOR ERRORS IN ALOG.
1014      */
1015     rlb = log((float)MP_BASE) * 1.01f;
1016
1017     /* NOW SAFE TO CONVERT X TO REAL */
1018     rx = mp_cast_to_float(x);
1019
1020     /* SAVE SIGN AND WORK WITH ABS(X) */
1021     xs = x->sign;
1022     mp_abs(x, &t2);
1023
1024     /* GET FRACTIONAL AND INTEGER PARTS OF ABS(X) */
1025     ix = mp_cast_to_int(&t2);
1026     mp_fractional_component(&t2, &t2);
1027
1028     /* ATTACH SIGN TO FRACTIONAL PART AND COMPUTE EXP OF IT */
1029     t2.sign *= xs;
1030     mp_exp(&t2, z);
1031
1032     /*  COMPUTE E-2 OR 1/E USING TWO EXTRA DIGITS IN CASE ABS(X) LARGE
1033      *  (BUT ONLY ONE EXTRA DIGIT IF T < 4)
1034      */
1035     if (MP_T < 4)
1036         tss = MP_T + 1;
1037     else
1038         tss = MP_T + 2;
1039
1040     /* LOOP FOR E COMPUTATION. DECREASE T IF POSSIBLE. */
1041     /* Computing MIN */
1042     mp_set_from_integer(xs, &t1);
1043
1044     t2.sign = 0;
1045     for (i = 2 ; ; i++) {
1046         if (min(tss, tss + 2 + t1.exponent) <= 2)
1047             break;
1048
1049         mp_divide_integer(&t1, i * xs, &t1);
1050         mp_add(&t2, &t1, &t2);
1051         if (mp_is_zero(&t1))
1052             break;
1053     }
1054
1055     /* RAISE E OR 1/E TO POWER IX */
1056     if (xs > 0)
1057         mp_add_integer(&t2, 2, &t2);
1058     mp_xpowy_integer(&t2, ix, &t2);
1059
1060     /* MULTIPLY EXPS OF INTEGER AND FRACTIONAL PARTS */
1061     mp_multiply(z, &t2, z);
1062
1063     /*  CHECK THAT RELATIVE ERROR LESS THAN 0.01 UNLESS ABS(X) LARGE
1064      *  (WHEN EXP MIGHT OVERFLOW OR UNDERFLOW)
1065      */
1066     if (fabs(rx) > 10.0f)
1067         return;
1068
1069     rz = mp_cast_to_float(z);
1070     r__1 = rz - exp(rx);
1071     if (fabs(r__1) < rz * 0.01f)
1072         return;
1073
1074     /*  THE FOLLOWING MESSAGE MAY INDICATE THAT
1075      *  B**(T-1) IS TOO SMALL, OR THAT M IS TOO SMALL SO THE
1076      *  RESULT UNDERFLOWED.
1077      */
1078     mperr("*** ERROR OCCURRED IN MP_EPOWY, RESULT INCORRECT ***");
1079 }
1080
1081
1082 void
1083 mp_epowy(const MPNumber *x, MPNumber *z)
1084 {
1085     /* e^0 = 1 */
1086     if (mp_is_zero(x)) {
1087         mp_set_from_integer(1, z);
1088         return;
1089     }
1090
1091     if (mp_is_complex(x)) {
1092         MPNumber x_real, r, theta;
1093
1094         mp_real_component(x, &x_real);
1095         mp_imaginary_component(x, &theta);
1096
1097         mp_epowy_real(&x_real, &r);
1098         mp_set_from_polar(&r, MP_RADIANS, &theta, z);
1099     }
1100     else
1101         mp_epowy_real(x, z);
1102 }
1103
1104
1105 /*  RETURNS K = K/GCD AND L = L/GCD, WHERE GCD IS THE
1106  *  GREATEST COMMON DIVISOR OF K AND L.
1107  *  SAVE INPUT PARAMETERS IN LOCAL VARIABLES
1108  */
1109 void
1110 mp_gcd(int64_t *k, int64_t *l)
1111 {
1112     int64_t i, j;
1113
1114     i = abs(*k);
1115     j = abs(*l);
1116     if (j == 0) {
1117         /* IF J = 0 RETURN (1, 0) UNLESS I = 0, THEN (0, 0) */
1118         *k = 1;
1119         *l = 0;
1120         if (i == 0)
1121             *k = 0;
1122         return;
1123     }
1124
1125     /* EUCLIDEAN ALGORITHM LOOP */
1126     do {
1127         i %= j;
1128         if (i == 0) {
1129             *k = *k / j;
1130             *l = *l / j;
1131             return;
1132         }
1133         j %= i;
1134     } while (j != 0);
1135
1136     /* HERE J IS THE GCD OF K AND L */
1137     *k = *k / i;
1138     *l = *l / i;
1139 }
1140
1141
1142 bool
1143 mp_is_zero(const MPNumber *x)
1144 {
1145     return x->sign == 0 && x->im_sign == 0;
1146 }
1147
1148
1149 bool
1150 mp_is_negative(const MPNumber *x)
1151 {
1152     return x->sign < 0;
1153 }
1154
1155
1156 bool
1157 mp_is_greater_equal(const MPNumber *x, const MPNumber *y)
1158 {
1159     return mp_compare_mp_to_mp(x, y) >= 0;
1160 }
1161
1162
1163 bool
1164 mp_is_greater_than(const MPNumber *x, const MPNumber *y)
1165 {
1166     return mp_compare_mp_to_mp(x, y) > 0;
1167 }
1168
1169
1170 bool
1171 mp_is_less_equal(const MPNumber *x, const MPNumber *y)
1172 {
1173     return mp_compare_mp_to_mp(x, y) <= 0;
1174 }
1175
1176
1177 /*  RETURNS MP Y = LN(1+X) IF X IS AN MP NUMBER SATISFYING THE
1178  *  CONDITION ABS(X) < 1/B, ERROR OTHERWISE.
1179  *  USES NEWTONS METHOD TO SOLVE THE EQUATION
1180  *  EXP1(-Y) = X, THEN REVERSES SIGN OF Y.
1181  */
1182 static void
1183 mp_lns(const MPNumber *x, MPNumber *z)
1184 {
1185     int t, it0;
1186     MPNumber t1, t2, t3;
1187
1188     /* ln(1+0) = 0 */
1189     if (mp_is_zero(x)) {
1190         mp_set_from_integer(0, z);
1191         return;
1192     }
1193
1194     /* Get starting approximation -ln(1+x) ~= -x + x^2/2 - x^3/3 + x^4/4 */
1195     mp_set_from_mp(x, &t2);
1196     mp_divide_integer(x, 4, &t1);
1197     mp_add_fraction(&t1, -1, 3, &t1);
1198     mp_multiply(x, &t1, &t1);
1199     mp_add_fraction(&t1, 1, 2, &t1);
1200     mp_multiply(x, &t1, &t1);
1201     mp_add_integer(&t1, -1, &t1);
1202     mp_multiply(x, &t1, z);
1203
1204     /* Solve using Newtons method */
1205     it0 = t = 5;
1206     while(1)
1207     {
1208         int ts2, ts3;
1209
1210         /* t3 = (e^t3 - 1) */
1211         /* z = z - (t2 + t3 + (t2 * t3)) */
1212         mp_epowy(z, &t3);
1213         mp_add_integer(&t3, -1, &t3);
1214         mp_multiply(&t2, &t3, &t1);
1215         mp_add(&t3, &t1, &t3);
1216         mp_add(&t2, &t3, &t3);
1217         mp_subtract(z, &t3, z);
1218         if (t >= MP_T)
1219             break;
1220
1221         /*  FOLLOWING LOOP COMPUTES NEXT VALUE OF T TO USE.
1222          *  BECAUSE NEWTONS METHOD HAS 2ND ORDER CONVERGENCE,
1223          *  WE CAN ALMOST DOUBLE T EACH TIME.
1224          */
1225         ts3 = t;
1226         t = MP_T;
1227         do {
1228             ts2 = t;
1229             t = (t + it0) / 2;
1230         } while (t > ts3);
1231         t = ts2;
1232     }
1233
1234     /* CHECK THAT NEWTON ITERATION WAS CONVERGING AS EXPECTED */
1235     if (t3.sign != 0 && t3.exponent << 1 > it0 - MP_T) {
1236         mperr("*** ERROR OCCURRED IN MP_LNS, NEWTON ITERATION NOT CONVERGING PROPERLY ***");
1237     }
1238
1239     z->sign = -z->sign;
1240 }
1241
1242
1243 static void
1244 mp_ln_real(const MPNumber *x, MPNumber *z)
1245 {
1246     int e, k;
1247     float rx, rlx;
1248     MPNumber t1, t2;
1249
1250     /* LOOP TO GET APPROXIMATE LN(X) USING SINGLE-PRECISION */
1251     mp_set_from_mp(x, &t1);
1252     mp_set_from_integer(0, z);
1253     for(k = 0; k < 10; k++)
1254     {
1255         /* COMPUTE FINAL CORRECTION ACCURATELY USING MP_LNS */
1256         mp_add_integer(&t1, -1, &t2);
1257         if (mp_is_zero(&t2) || t2.exponent + 1 <= 0) {
1258             mp_lns(&t2, &t2);
1259             mp_add(z, &t2, z);
1260             return;
1261         }
1262
1263         /* REMOVE EXPONENT TO AVOID FLOATING-POINT OVERFLOW */
1264         e = t1.exponent;
1265         t1.exponent = 0;
1266         rx = mp_cast_to_float(&t1);
1267         t1.exponent = e;
1268         rlx = log(rx) + (float)e * log((float)MP_BASE);
1269         mp_set_from_float(-(double)rlx, &t2);
1270
1271         /* UPDATE Z AND COMPUTE ACCURATE EXP OF APPROXIMATE LOG */
1272         mp_subtract(z, &t2, z);
1273         mp_epowy(&t2, &t2);
1274
1275         /* COMPUTE RESIDUAL WHOSE LOG IS STILL TO BE FOUND */
1276         mp_multiply(&t1, &t2, &t1);
1277     }
1278
1279     mperr("*** ERROR IN MP_LN, ITERATION NOT CONVERGING ***");
1280 }
1281
1282
1283 void
1284 mp_ln(const MPNumber *x, MPNumber *z)
1285 {
1286     /* ln(0) undefined */
1287     if (mp_is_zero(x)) {
1288         /* Translators: Error displayed when attempting to take logarithm of zero */
1289         mperr(_("Logarithm of zero is undefined"));
1290         mp_set_from_integer(0, z);
1291         return;
1292     }
1293
1294     /* ln(-x) complex */
1295     /* FIXME: Make complex numbers optional */
1296     /*if (mp_is_negative(x)) {
1297         // Translators: Error displayed attempted to take logarithm of negative value
1298         mperr(_("Logarithm of negative values is undefined"));
1299         mp_set_from_integer(0, z);
1300         return;
1301     }*/
1302
1303     if (mp_is_complex(x) || mp_is_negative(x)) {
1304         MPNumber r, theta, z_real;
1305
1306         /* ln(re^iθ) = e^(ln(r)+iθ) */
1307         mp_abs(x, &r);
1308         mp_arg(x, MP_RADIANS, &theta);
1309
1310         mp_ln_real(&r, &z_real);
1311         mp_set_from_complex(&z_real, &theta, z);
1312     }
1313     else
1314         mp_ln_real(x, z);
1315 }
1316
1317
1318 void
1319 mp_logarithm(int64_t n, const MPNumber *x, MPNumber *z)
1320 {
1321     MPNumber t1, t2;
1322
1323     /* log(0) undefined */
1324     if (mp_is_zero(x)) {
1325         /* Translators: Error displayed when attempting to take logarithm of zero */
1326         mperr(_("Logarithm of zero is undefined"));
1327         mp_set_from_integer(0, z);
1328         return;
1329     }
1330
1331     /* logn(x) = ln(x) / ln(n) */
1332     mp_set_from_integer(n, &t1);
1333     mp_ln(&t1, &t1);
1334     mp_ln(x, &t2);
1335     mp_divide(&t2, &t1, z);
1336 }
1337
1338
1339 bool
1340 mp_is_less_than(const MPNumber *x, const MPNumber *y)
1341 {
1342     return mp_compare_mp_to_mp(x, y) < 0;
1343 }
1344
1345
1346 static void
1347 mp_multiply_real(const MPNumber *x, const MPNumber *y, MPNumber *z)
1348 {
1349     int c, i, xi;
1350     MPNumber r;
1351
1352     /* x*0 = 0*y = 0 */
1353     if (x->sign == 0 || y->sign == 0) {
1354         mp_set_from_integer(0, z);
1355         return;
1356     }
1357
1358     z->sign = x->sign * y->sign;
1359     z->exponent = x->exponent + y->exponent;
1360     memset(&r, 0, sizeof(MPNumber));
1361
1362     /* PERFORM MULTIPLICATION */
1363     c = 8;
1364     for (i = 0; i < MP_T; i++) {
1365         int j;
1366
1367         xi = x->fraction[i];
1368
1369         /* FOR SPEED, PUT THE NUMBER WITH MANY ZEROS FIRST */
1370         if (xi == 0)
1371             continue;
1372
1373         /* Computing MIN */
1374         for (j = 0; j < min(MP_T, MP_T + 3 - i); j++)
1375             r.fraction[i+j+1] += xi * y->fraction[j];
1376         c--;
1377         if (c > 0)
1378             continue;
1379
1380         /* CHECK FOR LEGAL BASE B DIGIT */
1381         if (xi < 0 || xi >= MP_BASE) {
1382             mperr("*** ILLEGAL BASE B DIGIT IN CALL TO MP_MULTIPLY, POSSIBLE OVERWRITING PROBLEM ***");
1383             mp_set_from_integer(0, z);
1384             return;
1385         }
1386
1387         /*  PROPAGATE CARRIES AT END AND EVERY EIGHTH TIME,
1388          *  FASTER THAN DOING IT EVERY TIME.
1389          */
1390         for (j = MP_T + 3; j >= 0; j--) {
1391             int ri = r.fraction[j] + c;
1392             if (ri < 0) {
1393                 mperr("*** INTEGER OVERFLOW IN MP_MULTIPLY, B TOO LARGE ***");
1394                 mp_set_from_integer(0, z);
1395                 return;
1396             }
1397             c = ri / MP_BASE;
1398             r.fraction[j] = ri - MP_BASE * c;
1399         }
1400         if (c != 0) {
1401             mperr("*** ILLEGAL BASE B DIGIT IN CALL TO MP_MULTIPLY, POSSIBLE OVERWRITING PROBLEM ***");
1402             mp_set_from_integer(0, z);
1403             return;
1404         }
1405         c = 8;
1406     }
1407
1408     if (c != 8) {
1409         if (xi < 0 || xi >= MP_BASE) {
1410             mperr("*** ILLEGAL BASE B DIGIT IN CALL TO MP_MULTIPLY, POSSIBLE OVERWRITING PROBLEM ***");
1411             mp_set_from_integer(0, z);
1412             return;
1413         }
1414
1415         c = 0;
1416         for (i = MP_T + 3; i >= 0; i--) {
1417             int ri = r.fraction[i] + c;
1418             if (ri < 0) {
1419                 mperr("*** INTEGER OVERFLOW IN MP_MULTIPLY, B TOO LARGE ***");
1420                 mp_set_from_integer(0, z);
1421                 return;
1422             }
1423             c = ri / MP_BASE;
1424             r.fraction[i] = ri - MP_BASE * c;
1425         }
1426
1427         if (c != 0) {
1428             mperr("*** ILLEGAL BASE B DIGIT IN CALL TO MP_MULTIPLY, POSSIBLE OVERWRITING PROBLEM ***");
1429             mp_set_from_integer(0, z);
1430             return;
1431         }
1432     }
1433
1434     /* Clear complex part */
1435     z->im_sign = 0;
1436     z->im_exponent = 0;
1437     memset(z->im_fraction, 0, sizeof(int) * MP_SIZE);
1438
1439     /* NORMALIZE AND ROUND RESULT */
1440     // FIXME: Use stack variable because of mp_normalize brokeness
1441     for (i = 0; i < MP_SIZE; i++)
1442         z->fraction[i] = r.fraction[i];
1443     mp_normalize(z);
1444 }
1445
1446
1447 void
1448 mp_multiply(const MPNumber *x, const MPNumber *y, MPNumber *z)
1449 {
1450     /* x*0 = 0*y = 0 */
1451     if (mp_is_zero(x) || mp_is_zero(y)) {
1452         mp_set_from_integer(0, z);
1453         return;
1454     }
1455
1456     /* (a+bi)(c+di) = (ac-bd)+(ad+bc)i */
1457     if (mp_is_complex(x) || mp_is_complex(y)) {
1458         MPNumber real_x, real_y, im_x, im_y, t1, t2, real_z, im_z;
1459
1460         mp_real_component(x, &real_x);
1461         mp_imaginary_component(x, &im_x);
1462         mp_real_component(y, &real_y);
1463         mp_imaginary_component(y, &im_y);
1464
1465         mp_multiply_real(&real_x, &real_y, &t1);
1466         mp_multiply_real(&im_x, &im_y, &t2);
1467         mp_subtract(&t1, &t2, &real_z);
1468
1469         mp_multiply_real(&real_x, &im_y, &t1);
1470         mp_multiply_real(&im_x, &real_y, &t2);
1471         mp_add(&t1, &t2, &im_z);
1472
1473         mp_set_from_complex(&real_z, &im_z, z);
1474     }
1475     else {
1476         mp_multiply_real(x, y, z);
1477     }
1478 }
1479
1480
1481 static void
1482 mp_multiply_integer_real(const MPNumber *x, int64_t y, MPNumber *z)
1483 {
1484     int c, i;
1485     MPNumber x_copy;
1486
1487     /* x*0 = 0*y = 0 */
1488     if (mp_is_zero(x) || y == 0) {
1489         mp_set_from_integer(0, z);
1490         return;
1491     }
1492
1493     /* x*1 = x, x*-1 = -x */
1494     // FIXME: Why is this not working? mp_ext is using this function to do a normalization
1495     /*if (y == 1 || y == -1) {
1496         if (y < 0)
1497             mp_invert_sign(x, z);
1498         else
1499             mp_set_from_mp(x, z);
1500         return;
1501     }*/
1502
1503     /* Copy x as z may also refer to x */
1504     mp_set_from_mp(x, &x_copy);
1505     mp_set_from_integer(0, z);
1506
1507     if (y < 0) {
1508         y = -y;
1509         z->sign = -x_copy.sign;
1510     }
1511     else
1512         z->sign = x_copy.sign;
1513     z->exponent = x_copy.exponent + 4;
1514
1515     /* FORM PRODUCT IN ACCUMULATOR */
1516     c = 0;
1517
1518     /*  IF y*B NOT REPRESENTABLE AS AN INTEGER WE HAVE TO SIMULATE
1519      *  DOUBLE-PRECISION MULTIPLICATION.
1520      */
1521
1522     /* Computing MAX */
1523     if (y >= max(MP_BASE << 3, 32767 / MP_BASE)) {
1524         int64_t j1, j2;
1525
1526         /* HERE J IS TOO LARGE FOR SINGLE-PRECISION MULTIPLICATION */
1527         j1 = y / MP_BASE;
1528         j2 = y - j1 * MP_BASE;
1529
1530         /* FORM PRODUCT */
1531         for (i = MP_T + 3; i >= 0; i--) {
1532             int64_t c1, c2, is, ix, t;
1533
1534             c1 = c / MP_BASE;
1535             c2 = c - MP_BASE * c1;
1536             ix = 0;
1537             if (i > 3)
1538                 ix = x_copy.fraction[i - 4];
1539
1540             t = j2 * ix + c2;
1541             is = t / MP_BASE;
1542             c = j1 * ix + c1 + is;
1543             z->fraction[i] = t - MP_BASE * is;
1544         }
1545     }
1546     else
1547     {
1548         int64_t ri = 0;
1549
1550         for (i = MP_T + 3; i >= 4; i--) {
1551             ri = y * x_copy.fraction[i - 4] + c;
1552             c = ri / MP_BASE;
1553             z->fraction[i] = ri - MP_BASE * c;
1554         }
1555
1556         /* CHECK FOR INTEGER OVERFLOW */
1557         if (ri < 0) {
1558             mperr("*** INTEGER OVERFLOW IN mp_multiply_integer, B TOO LARGE ***");
1559             mp_set_from_integer(0, z);
1560             return;
1561         }
1562
1563         /* HAVE TO TREAT FIRST FOUR WORDS OF R SEPARATELY */
1564         for (i = 3; i >= 0; i--) {
1565             int t;
1566
1567             t = c;
1568             c = t / MP_BASE;
1569             z->fraction[i] = t - MP_BASE * c;
1570         }
1571     }
1572
1573     /* HAVE TO SHIFT RIGHT HERE AS CARRY OFF END */
1574     while (c != 0) {
1575         int64_t t;
1576
1577         for (i = MP_T + 3; i >= 1; i--)
1578             z->fraction[i] = z->fraction[i - 1];
1579         t = c;
1580         c = t / MP_BASE;
1581         z->fraction[0] = t - MP_BASE * c;
1582         z->exponent++;
1583     }
1584
1585     if (c < 0) {
1586         mperr("*** INTEGER OVERFLOW IN mp_multiply_integer, B TOO LARGE ***");
1587         mp_set_from_integer(0, z);
1588         return;
1589     }
1590
1591     z->im_sign = 0;
1592     z->im_exponent = 0;
1593     memset(z->im_fraction, 0, sizeof(int) * MP_SIZE);
1594     mp_normalize(z);
1595 }
1596
1597
1598 void
1599 mp_multiply_integer(const MPNumber *x, int64_t y, MPNumber *z)
1600 {
1601     if (mp_is_complex(x)) {
1602         MPNumber re_z, im_z;
1603         mp_real_component(x, &re_z);
1604         mp_imaginary_component(x, &im_z);
1605         mp_multiply_integer_real(&re_z, y, &re_z);
1606         mp_multiply_integer_real(&im_z, y, &im_z);
1607         mp_set_from_complex(&re_z, &im_z, z);
1608     }
1609     else
1610         mp_multiply_integer_real(x, y, z);
1611 }
1612
1613
1614 void
1615 mp_multiply_fraction(const MPNumber *x, int64_t numerator, int64_t denominator, MPNumber *z)
1616 {
1617     if (denominator == 0) {
1618         mperr(_("Division by zero is undefined"));
1619         mp_set_from_integer(0, z);
1620         return;
1621     }
1622
1623     if (numerator == 0) {
1624         mp_set_from_integer(0, z);
1625         return;
1626     }
1627
1628     /* Reduce to lowest terms */
1629     mp_gcd(&numerator, &denominator);
1630     mp_divide_integer(x, denominator, z);
1631     mp_multiply_integer(z, numerator, z);
1632 }
1633
1634
1635 void
1636 mp_invert_sign(const MPNumber *x, MPNumber *z)
1637 {
1638     mp_set_from_mp(x, z);
1639     z->sign = -z->sign;
1640     z->im_sign = -z->im_sign;
1641 }
1642
1643
1644 // FIXME: Is r->fraction large enough?  It seems to be in practise but it may be MP_T+4 instead of MP_T
1645 // FIXME: There is some sort of stack corruption/use of unitialised variables here.  Some functions are
1646 // using stack variables as x otherwise there are corruption errors. e.g. "Cos(45) - 1/Sqrt(2) = -0"
1647 // (try in scientific mode)
1648 void
1649 mp_normalize(MPNumber *x)
1650 {
1651     int start_index;
1652
1653     /* Find first non-zero digit */
1654     for (start_index = 0; start_index < MP_SIZE && x->fraction[start_index] == 0; start_index++);
1655
1656     /* Mark as zero */
1657     if (start_index >= MP_SIZE) {
1658         x->sign = 0;
1659         x->exponent = 0;
1660         return;
1661     }
1662
1663     /* Shift left so first digit is non-zero */
1664     if (start_index > 0) {
1665         int i;
1666
1667         x->exponent -= start_index;
1668         for (i = 0; (i + start_index) < MP_SIZE; i++)
1669             x->fraction[i] = x->fraction[i + start_index];
1670         for (; i < MP_SIZE; i++)
1671             x->fraction[i] = 0;
1672     }
1673 }
1674
1675
1676 static void
1677 mp_pwr(const MPNumber *x, const MPNumber *y, MPNumber *z)
1678 {
1679     MPNumber t;
1680
1681     /* (-x)^y imaginary */
1682     /* FIXME: Make complex numbers optional */
1683     /*if (x->sign < 0) {
1684         mperr(_("The power of negative numbers is only defined for integer exponents"));
1685         mp_set_from_integer(0, z);
1686         return;
1687     }*/
1688
1689     /* 0^-y illegal */
1690     if (mp_is_zero(x) && y->sign < 0) {
1691         mperr(_("The power of zero is undefined for a negative exponent"));
1692         mp_set_from_integer(0, z);
1693         return;
1694     }
1695
1696     /* x^0 = 1 */
1697     if (mp_is_zero(y)) {
1698         mp_set_from_integer(1, z);
1699         return;
1700     }
1701
1702     mp_ln(x, &t);
1703     mp_multiply(y, &t, z);
1704     mp_epowy(z, z);
1705 }
1706
1707
1708 static void
1709 mp_reciprocal_real(const MPNumber *x, MPNumber *z)
1710 {
1711     MPNumber t1, t2;
1712     int it0, t;
1713
1714     /* 1/0 invalid */
1715     if (mp_is_zero(x)) {
1716         mperr(_("Reciprocal of zero is undefined"));
1717         mp_set_from_integer(0, z);
1718         return;
1719     }
1720
1721     /* Start by approximating value using floating point */
1722     mp_set_from_mp(x, &t1);
1723     t1.exponent = 0;
1724     mp_set_from_float(1.0f / mp_cast_to_float(&t1), &t1);
1725     t1.exponent -= x->exponent;
1726
1727     it0 = t = 3;
1728     while(1) {
1729         int ts2, ts3;
1730
1731         /* t1 = t1 - (t1 * ((x * t1) - 1))   (2*t1 - t1^2*x) */
1732         mp_multiply(x, &t1, &t2);
1733         mp_add_integer(&t2, -1, &t2);
1734         mp_multiply(&t1, &t2, &t2);
1735         mp_subtract(&t1, &t2, &t1);
1736         if (t >= MP_T)
1737             break;
1738
1739         /*  FOLLOWING LOOP ALMOST DOUBLES T (POSSIBLE
1740          *  BECAUSE NEWTONS METHOD HAS 2ND ORDER CONVERGENCE).
1741          */
1742         ts3 = t;
1743         t = MP_T;
1744         do {
1745             ts2 = t;
1746             t = (t + it0) / 2;
1747         } while (t > ts3);
1748         t = min(ts2, MP_T);
1749     }
1750
1751     /* RETURN IF NEWTON ITERATION WAS CONVERGING */
1752     if (t2.sign != 0 && (t1.exponent - t2.exponent) << 1 < MP_T - it0) {
1753         /*  THE FOLLOWING MESSAGE MAY INDICATE THAT B**(T-1) IS TOO SMALL,
1754          *  OR THAT THE STARTING APPROXIMATION IS NOT ACCURATE ENOUGH.
1755          */
1756         mperr("*** ERROR OCCURRED IN MP_RECIPROCAL, NEWTON ITERATION NOT CONVERGING PROPERLY ***");
1757     }
1758
1759     mp_set_from_mp(&t1, z);
1760 }
1761
1762
1763 void
1764 mp_reciprocal(const MPNumber *x, MPNumber *z)
1765 {
1766     if (mp_is_complex(x)) {
1767         MPNumber t1, t2;
1768         MPNumber real_x, im_x;
1769
1770         mp_real_component(x, &real_x);
1771         mp_imaginary_component(x, &im_x);
1772
1773         /* 1/(a+bi) = (a-bi)/(a+bi)(a-bi) = (a-bi)/(a²+b²) */
1774         mp_multiply(&real_x, &real_x, &t1);
1775         mp_multiply(&im_x, &im_x, &t2);
1776         mp_add(&t1, &t2, &t1);
1777         mp_reciprocal_real(&t1, z);
1778         mp_conjugate(x, &t1);
1779         mp_multiply(&t1, z, z);
1780     }
1781     else
1782         mp_reciprocal_real(x, z);
1783 }
1784
1785
1786 static void
1787 mp_root_real(const MPNumber *x, int64_t n, MPNumber *z)
1788 {
1789     float approximation;
1790     int ex, np, it0, t;
1791     MPNumber t1, t2;
1792
1793     /* x^(1/1) = x */
1794     if (n == 1) {
1795         mp_set_from_mp(x, z);
1796         return;
1797     }
1798
1799     /* x^(1/0) invalid */
1800     if (n == 0) {
1801         mperr(_("Root must be non-zero"));
1802         mp_set_from_integer(0, z);
1803         return;
1804     }
1805
1806     np = abs(n);
1807
1808     /* LOSS OF ACCURACY IF NP LARGE, SO ONLY ALLOW NP <= MAX (B, 64) */
1809     if (np > max(MP_BASE, 64)) {
1810         mperr("*** ABS(N) TOO LARGE IN CALL TO MP_ROOT ***");
1811         mp_set_from_integer(0, z);
1812         return;
1813     }
1814
1815     /* 0^(1/n) = 0 for positive n */
1816     if (mp_is_zero(x)) {
1817         mp_set_from_integer(0, z);
1818         if (n <= 0)
1819             mperr(_("Negative root of zero is undefined"));
1820         return;
1821     }
1822
1823     // FIXME: Imaginary root
1824     if (x->sign < 0 && np % 2 == 0) {
1825         mperr(_("nth root of negative number is undefined for even n"));
1826         mp_set_from_integer(0, z);
1827         return;
1828     }
1829
1830     /* DIVIDE EXPONENT BY NP */
1831     ex = x->exponent / np;
1832
1833     /* Get initial approximation */
1834     mp_set_from_mp(x, &t1);
1835     t1.exponent = 0;
1836     approximation = exp(((float)(np * ex - x->exponent) * log((float)MP_BASE) -
1837                          log((fabs(mp_cast_to_float(&t1))))) / (float)np);
1838     mp_set_from_float(approximation, &t1);
1839     t1.sign = x->sign;
1840     t1.exponent -= ex;
1841
1842     /* MAIN ITERATION LOOP */
1843     it0 = t = 3;
1844     while(1) {
1845         int ts2, ts3;
1846
1847         /* t1 = t1 - ((t1 * ((x * t1^np) - 1)) / np) */
1848         mp_xpowy_integer(&t1, np, &t2);
1849         mp_multiply(x, &t2, &t2);
1850         mp_add_integer(&t2, -1, &t2);
1851         mp_multiply(&t1, &t2, &t2);
1852         mp_divide_integer(&t2, np, &t2);
1853         mp_subtract(&t1, &t2, &t1);
1854
1855         /*  FOLLOWING LOOP ALMOST DOUBLES T (POSSIBLE BECAUSE
1856          *  NEWTONS METHOD HAS 2ND ORDER CONVERGENCE).
1857          */
1858         if (t >= MP_T)
1859             break;
1860
1861         ts3 = t;
1862         t = MP_T;
1863         do {
1864             ts2 = t;
1865             t = (t + it0) / 2;
1866         } while (t > ts3);
1867         t = min(ts2, MP_T);
1868     }
1869
1870     /*  NOW R(I2) IS X**(-1/NP)
1871      *  CHECK THAT NEWTON ITERATION WAS CONVERGING
1872      */
1873     if (t2.sign != 0 && (t1.exponent - t2.exponent) << 1 < MP_T - it0) {
1874         /*  THE FOLLOWING MESSAGE MAY INDICATE THAT B**(T-1) IS TOO SMALL,
1875          *  OR THAT THE INITIAL APPROXIMATION OBTAINED USING ALOG AND EXP
1876          *  IS NOT ACCURATE ENOUGH.
1877          */
1878         mperr("*** ERROR OCCURRED IN MP_ROOT, NEWTON ITERATION NOT CONVERGING PROPERLY ***");
1879     }
1880
1881     if (n >= 0) {
1882         mp_xpowy_integer(&t1, n - 1, &t1);
1883         mp_multiply(x, &t1, z);
1884         return;
1885     }
1886
1887     mp_set_from_mp(&t1, z);
1888 }
1889
1890
1891 void
1892 mp_root(const MPNumber *x, int64_t n, MPNumber *z)
1893 {
1894     if (!mp_is_complex(x) && mp_is_negative(x) && n % 2 == 1) {
1895         mp_abs(x, z);
1896         mp_root_real(z, n, z);
1897         mp_invert_sign(z, z);
1898     }
1899     else if (mp_is_complex(x) || mp_is_negative(x)) {
1900         MPNumber r, theta;
1901
1902         mp_abs(x, &r);
1903         mp_arg(x, MP_RADIANS, &theta);
1904
1905         mp_root_real(&r, n, &r);
1906         mp_divide_integer(&theta, n, &theta);
1907         mp_set_from_polar(&r, MP_RADIANS, &theta, z);
1908     }
1909     else
1910         mp_root_real(x, n, z);
1911 }
1912
1913
1914 void
1915 mp_sqrt(const MPNumber *x, MPNumber *z)
1916 {
1917     if (mp_is_zero(x))
1918        mp_set_from_integer(0, z);
1919     /* FIXME: Make complex numbers optional */
1920     /*else if (x->sign < 0) {
1921         mperr(_("Square root is undefined for negative values"));
1922         mp_set_from_integer(0, z);
1923     }*/
1924     else {
1925         MPNumber t;
1926
1927         mp_root(x, -2, &t);
1928         mp_multiply(x, &t, z);
1929         mp_ext(t.fraction[0], z->fraction[0], z);
1930     }
1931 }
1932
1933
1934 void
1935 mp_factorial(const MPNumber *x, MPNumber *z)
1936 {
1937     int i, value;
1938
1939     /* 0! == 1 */
1940     if (mp_is_zero(x)) {
1941         mp_set_from_integer(1, z);
1942         return;
1943     }
1944     if (!mp_is_natural(x)) {
1945         /* Translators: Error displayed when attempted take the factorial of a fractional number */
1946         mperr(_("Factorial is only defined for natural numbers"));
1947         mp_set_from_integer(0, z);
1948         return;
1949     }
1950
1951     /* Convert to integer - if couldn't be converted then the factorial would be too big anyway */
1952     value = mp_cast_to_int(x);
1953     mp_set_from_mp(x, z);
1954     for (i = 2; i < value; i++)
1955         mp_multiply_integer(z, i, z);
1956 }
1957
1958
1959 void
1960 mp_modulus_divide(const MPNumber *x, const MPNumber *y, MPNumber *z)
1961 {
1962     MPNumber t1, t2;
1963
1964     if (!mp_is_integer(x) || !mp_is_integer(y)) {
1965         /* Translators: Error displayed when attemping to do a modulus division on non-integer numbers */
1966         mperr(_("Modulus division is only defined for integers"));
1967         mp_set_from_integer(0, z);
1968     }
1969
1970     mp_divide(x, y, &t1);
1971     mp_floor(&t1, &t1);
1972     mp_multiply(&t1, y, &t2);
1973     mp_subtract(x, &t2, z);
1974
1975     mp_set_from_integer(0, &t1);
1976     if ((mp_is_less_than(y, &t1) && mp_is_greater_than(z, &t1)) || mp_is_less_than(z, &t1))
1977         mp_add(z, y, z);
1978 }
1979
1980
1981 void
1982 mp_xpowy(const MPNumber *x, const MPNumber *y, MPNumber *z)
1983 {
1984     if (mp_is_integer(y)) {
1985         mp_xpowy_integer(x, mp_cast_to_int(y), z);
1986     } else {
1987         MPNumber reciprocal;
1988         mp_reciprocal(y, &reciprocal);
1989         if (mp_is_integer(&reciprocal))
1990             mp_root(x, mp_cast_to_int(&reciprocal), z);
1991         else
1992             mp_pwr(x, y, z);
1993     }
1994 }
1995
1996
1997 void
1998 mp_xpowy_integer(const MPNumber *x, int64_t n, MPNumber *z)
1999 {
2000     int i;
2001     MPNumber t;
2002
2003     /* 0^-n invalid */
2004     if (mp_is_zero(x) && n < 0) {
2005         /* Translators: Error displayed when attempted to raise 0 to a negative exponent */
2006         mperr(_("The power of zero is undefined for a negative exponent"));
2007         mp_set_from_integer(0, z);
2008         return;
2009     }
2010
2011     /* x^0 = 1 */
2012     if (n == 0) {
2013         mp_set_from_integer(1, z);
2014         return;
2015     }
2016
2017     /* 0^n = 0 */
2018     if (mp_is_zero(x)) {
2019         mp_set_from_integer(0, z);
2020         return;
2021     }
2022
2023     /* x^1 = x */
2024     if (n == 1) {
2025         mp_set_from_mp(x, z);
2026         return;
2027     }
2028
2029     if (n < 0) {
2030         mp_reciprocal(x, &t);
2031         n = -n;
2032     }
2033     else
2034         mp_set_from_mp(x, &t);
2035
2036     /* Multply x n times */
2037     // FIXME: Can do mp_multiply(z, z, z) until close to answer (each call doubles number of multiples) */
2038     mp_set_from_integer(1, z);
2039     for (i = 0; i < n; i++)
2040         mp_multiply(z, &t, z);
2041 }
2042
2043
2044 GList*
2045 mp_factorize(const MPNumber *x)
2046 {
2047     GList *list = NULL;
2048     MPNumber *factor = g_slice_alloc0(sizeof(MPNumber));
2049     MPNumber value, tmp, divisor, root;
2050
2051     mp_abs(x, &value);
2052
2053     if (mp_is_zero(&value)) {
2054         mp_set_from_mp(&value, factor);
2055         list = g_list_append(list, factor);
2056         return list;
2057     }
2058
2059     mp_set_from_integer(1, &tmp);
2060     if (mp_is_equal(&value, &tmp)) {
2061         mp_set_from_mp(x, factor);
2062         list = g_list_append(list, factor);
2063         return list;
2064     }
2065
2066     mp_set_from_integer(2, &divisor);
2067     while (TRUE) {
2068         mp_divide(&value, &divisor, &tmp);
2069         if (mp_is_integer(&tmp)) {
2070             value = tmp;
2071             mp_set_from_mp(&divisor, factor);
2072             list = g_list_append(list, factor);
2073             factor = g_slice_alloc0(sizeof(MPNumber));
2074         } else {
2075             break;
2076         }
2077     }
2078
2079     mp_set_from_integer(3, &divisor);
2080     mp_sqrt(&value, &root);
2081     while (mp_is_less_equal(&divisor, &root)) {
2082         mp_divide(&value, &divisor, &tmp);
2083         if (mp_is_integer(&tmp)) {
2084             value = tmp;
2085             mp_sqrt(&value, &root);
2086             mp_set_from_mp(&divisor, factor);
2087             list = g_list_append(list, factor);
2088             factor = g_slice_alloc0(sizeof(MPNumber));
2089         } else {
2090             mp_add_integer(&divisor, 2, &tmp);
2091             divisor = tmp;
2092         }
2093     }
2094
2095     mp_set_from_integer(1, &tmp);
2096     if (mp_is_greater_than(&value, &tmp)) {
2097         mp_set_from_mp(&value, factor);
2098         list = g_list_append(list, factor);
2099     } else {
2100         g_slice_free(MPNumber, factor);
2101     }
2102
2103     if (mp_is_negative(x)) {
2104         mp_invert_sign(list->data, list->data);
2105     }
2106
2107     return list;
2108 }