src/mp.c

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