gcalctool/mpmath.c

   1
   2 /*  Copyright (c) 1987-2006 Sun Microsystems, Inc. All Rights Reserved.
   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 <assert.h>
  21 #include <errno.h>
  22
  23 #include "calctool.h"    /* FIXME: include only needed stuff. */
  24 #include "mpmath.h"
  25 #include "extern.h"
  26
  27 static void mpacos(int *, int *);
  28 static void mpacosh(int *, int *);
  29 static void mpasinh(int *, int *);
  30 static void mpatanh(int *, int *);
  31
  32
  33 BOOLEAN
  34 ibool(double x)
  35 {
  36     BOOLEAN p = (BOOLEAN) x;
  37
  38     return(p);
  39 }
  40
  41
  42 double
  43 setbool(BOOLEAN p)
  44 {
  45     BOOLEAN q;
  46     double val;
  47
  48     q = p & 0x80000000;
  49     p &= 0x7fffffff;
  50     val = p;
  51     if (q) {
  52         val += 2147483648.0;
  53     }
  54
  55     return(val);
  56 }
  57
  58
  59 void
  60 calc_and(int t[MP_SIZE], int s1[MP_SIZE], int s2[MP_SIZE])
  61 {
  62     double dres, dval;
  63
  64     mpcmd(s1, &dres);
  65     mpcmd(s2, &dval);
  66     dres = setbool(ibool(dres) & ibool(dval));
  67     mpcdm(&dres, t);
  68 }
  69
  70
  71 void
  72 calc_or(int t[MP_SIZE], int s1[MP_SIZE], int s2[MP_SIZE])
  73 {
  74     double dres, dval;
  75
  76     mpcmd(s1, &dres);
  77     mpcmd(s2, &dval);
  78     dres = setbool(ibool(dres) | ibool(dval));
  79     mpcdm(&dres, t);
  80 }
  81
  82
  83 void
  84 calc_xor(int t[MP_SIZE], int s1[MP_SIZE], int s2[MP_SIZE])
  85 {
  86     double dres, dval;
  87
  88     mpcmd(s1, &dres);
  89     mpcmd(s2, &dval);
  90     dres = setbool(ibool(dres) ^ ibool(dval));
  91     mpcdm(&dres, t);
  92 }
  93
  94
  95 void
  96 calc_xnor(int t[MP_SIZE], int s1[MP_SIZE], int s2[MP_SIZE])
  97 {
  98     double dres, dval;
  99
 100     mpcmd(s1, &dres);
 101     mpcmd(s2, &dval);
 102     dres = setbool(~ibool(dres) ^ ibool(dval));
 103     mpcdm(&dres, t);
 104 }
 105
 106
 107 void
 108 calc_not(int s1[MP_SIZE], int t[MP_SIZE])
 109 {
 110     double dval;
 111
 112     mpcmd(s1, &dval);
 113     dval = setbool(~ibool(dval));
 114     mpcdm(&dval, t);
 115 }
 116
 117
 118 void
 119 calc_rand(int t[MP_SIZE])
 120 {
 121     double dval = drand48();
 122
 123     mpcdm(&dval, t);
 124 }
 125
 126
 127 void
 128 calc_u32(int s1[MP_SIZE], int t1[MP_SIZE])
 129 {
 130     double dval;
 131
 132     mpcmd(s1, &dval);
 133     dval = setbool(ibool(dval));
 134     mpcdm(&dval, t1);
 135 }
 136
 137
 138 void
 139 calc_u16(int s1[MP_SIZE], int t1[MP_SIZE])
 140 {
 141     double dval;
 142
 143     mpcmd(s1, &dval);
 144     dval = setbool(ibool(dval) & 0xffff);
 145     mpcdm(&dval, t1);
 146 }
 147
 148
 149 void
 150 calc_inv(int s1[MP_SIZE], int t1[MP_SIZE])     /* Calculate 1/x */
 151 {
 152     int MP1[MP_SIZE];
 153     int MP2[MP_SIZE];
 154     int i = 1;
 155
 156     mpcim(&i, MP1);
 157     mpstr(s1, MP2);
 158     mpdiv(MP1, MP2, t1);
 159 }
 160
 161
 162 void
 163 calc_tenpowx(int s1[MP_SIZE], int t1[MP_SIZE])   /* Calculate 10^x */
 164 {
 165     int MP1[MP_SIZE];
 166     int i = 10;
 167
 168     mpcim(&i, MP1);
 169     mppwr2(MP1, s1, t1);
 170 }
 171
 172
 173 void
 174 calc_xpowy(int MPx[MP_SIZE], int MPy[MP_SIZE], int MPres[MP_SIZE]) /* Do x^y */
 175 {
 176     int MP0[MP_SIZE], val;
 177
 178     val = 0;
 179     mpcim(&val, MP0);
 180     if (mplt(MPx, MP0)) {          /* Is x < 0 ? */
 181         int MPtmp[MP_SIZE];
 182
 183         mpcmim(MPy, MPtmp);
 184         if (mpeq(MPtmp, MPy)) {   /* Is y == int(y) ? */
 185             int y;
 186
 187             mpcmi(MPy, &y);
 188             mppwr(MPx, &y, MPres);
 189         } else {        /* y != int(y). Force mppwr2 to generate an error. */
 190             mppwr2(MPx, MPy, MPres);
 191         }
 192     } else {
 193         mppwr2(MPx, MPy, MPres);
 194     }
 195 }
 196
 197
 198 void
 199 calc_xtimestenpowx(int s1[MP_SIZE], int s2[MP_SIZE], int t1[MP_SIZE])
 200 {
 201     int MP1[MP_SIZE];
 202
 203     calc_tenpowx(s2, MP1);
 204     mpmul(s1, MP1, t1);
 205 }
 206
 207 int
 208 calc_modulus(int op1[MP_SIZE],
 209              int op2[MP_SIZE],
 210              int result[MP_SIZE])
 211 {
 212     int MP1[MP_SIZE], MP2[MP_SIZE];
 213
 214     mpcmim(op1, MP1);
 215     mpcmim(op2, MP2);
 216     if (!mpeq(op1, MP1) || !mpeq(op2, MP2)) {
 217         return -EINVAL;
 218     }
 219
 220     mpdiv(op1, op2, MP1);
 221     mpcmim(MP1, MP1);
 222     mpmul(MP1, op2, MP2);
 223     mpsub(op1, MP2, result);
 224     return 0;
 225 }
 226
 227 void
 228 calc_percent(int s1[MP_SIZE], int s2[MP_SIZE], int t1[MP_SIZE])
 229 {
 230     int MP1[MP_SIZE];
 231     int MP2[MP_SIZE];
 232
 233     mpmul(s1, s2, MP1);
 234     MPstr_to_num("0.01", DEC, MP2);
 235     mpmul(MP1, MP2, t1);
 236 }
 237
 238
 239 void
 240 do_zero(int t1[MP_SIZE])
 241 {
 242     int i = 0;
 243
 244     mpcim(&i, t1);
 245 }
 246
 247
 248 void
 249 do_e(int t1[MP_SIZE])
 250 {
 251     double e = 2.71828182846;
 252
 253     mpcdm(&e, t1);
 254 }
 255
 256
 257 static void
 258 mptan(int s1[MP_SIZE], int t1[MP_SIZE])
 259 {
 260     int MPcos[MP_SIZE];
 261     int MPsin[MP_SIZE];
 262     double cval;
 263
 264     mpsin(s1, MPsin);
 265     mpcos(s1, MPcos);
 266     mpcmd(MPcos, &cval);
 267     if (cval == 0.0) {
 268         doerr(_("Error, cannot calculate cosine"));
 269     }
 270     mpdiv(MPsin, MPcos, t1);
 271 }
 272
 273
 274 /* Change type to radian */
 275
 276 static void
 277 to_rad(int s1[MP_SIZE], int t1[MP_SIZE])
 278 {
 279     int i, MP1[MP_SIZE], MP2[MP_SIZE];
 280
 281     if (v->ttype == DEG) {
 282         mppi(MP1);
 283         mpmul(s1, MP1, MP2);
 284         i = 180;
 285         mpcim(&i, MP1);
 286         mpdiv(MP2, MP1, t1);
 287     } else if (v->ttype == GRAD) {
 288         mppi(MP1);
 289         mpmul(s1, MP1, MP2);
 290         i = 200;
 291         mpcim(&i, MP1);
 292         mpdiv(MP2, MP1, t1);
 293     } else {
 294         mpstr(s1, t1);
 295     }
 296 }
 297
 298
 299 static void
 300 do_trig_typeconv(enum trig_type ttype, int s1[MP_SIZE], int t1[MP_SIZE])
 301 {
 302     int i, MP1[MP_SIZE], MP2[MP_SIZE];
 303
 304     switch (ttype) {
 305
 306         case DEG:
 307             i = 180;
 308             mpcim(&i, MP1);
 309             mpmul(s1, MP1, MP2);
 310             mppi(MP1);
 311             mpdiv(MP2, MP1, t1);
 312             break;
 313
 314         case RAD:
 315             mpstr(s1, t1);
 316             break;
 317
 318         case GRAD:
 319             i = 200;
 320             mpcim(&i, MP1);
 321             mpmul(s1, MP1, MP2);
 322             mppi(MP1);
 323             mpdiv(MP2, MP1, t1);
 324             break;
 325
 326         default:
 327             assert(0);
 328             break;
 329     }
 330 }
 331
 332
 333 /* Calculate selected trigonometric function */
 334
 335 int
 336 calc_trigfunc(enum trigfunc_type type, int s1[MP_SIZE], int t1[MP_SIZE])
 337 {
 338     switch (type) {
 339         case sin_t:
 340             to_rad(s1, s1);
 341             mpsin(s1, t1);
 342             break;
 343
 344         case cos_t:
 345             to_rad(s1, s1);
 346             mpcos(s1, t1);
 347             break;
 348
 349         case tan_t:
 350             to_rad(s1, s1);
 351             mptan(s1, t1);
 352             break;
 353
 354         case sinh_t:
 355             mpsinh(s1, t1);
 356             break;
 357
 358         case cosh_t:
 359             mpcosh(s1, t1);
 360             break;
 361
 362         case tanh_t:
 363             mptanh(s1, t1);
 364             break;
 365
 366         case asin_t:
 367             mpasin(s1, t1);
 368             do_trig_typeconv(v->ttype, t1, t1);
 369             break;
 370
 371         case acos_t:
 372             mpacos(s1, t1);
 373             do_trig_typeconv(v->ttype, t1, t1);
 374             break;
 375
 376         case atan_t:
 377             mpatan(s1, t1);
 378             do_trig_typeconv(v->ttype, t1, t1);
 379             break;
 380
 381         case asinh_t:
 382             mpasinh(s1, t1);
 383             break;
 384
 385         case acosh_t:
 386             mpacosh(s1, t1);
 387             break;
 388
 389         case atanh_t:
 390             mpatanh(s1, t1);
 391             break;
 392     }
 393
 394     return(0);
 395 }
 396
 397
 398 /*  The following MP routines were not in the Brent FORTRAN package. They are
 399  *  derived here, in terms of the existing routines.
 400  */
 401
 402 /*  MP precision arc cosine.
 403  *
 404  *  1. If (x < -1.0  or x > 1.0) then report DOMAIN error and return 0.0.
 405  *
 406  *  2. If (x = 0.0) then acos(x) = PI/2.
 407  *
 408  *  3. If (x = 1.0) then acos(x) = 0.0
 409  *
 410  *  4. If (x = -1.0) then acos(x) = PI.
 411  *
 412  *  5. If (0.0 < x < 1.0) then  acos(x) = atan(sqrt(1-(x**2)) / x)
 413  *
 414  *  6. If (-1.0 < x < 0.0) then acos(x) = atan(sqrt(1-(x**2)) / x) + PI
 415  */
 416
 417 static void
 418 mpacos(int *MPx, int *MPretval)
 419 {
 420     int MP0[MP_SIZE],  MP1[MP_SIZE],  MP2[MP_SIZE];
 421     int MPn1[MP_SIZE], MPpi[MP_SIZE], MPy[MP_SIZE], val;
 422
 423     mppi(MPpi);
 424     val = 0;
 425     mpcim(&val, MP0);
 426     val = 1;
 427     mpcim(&val, MP1);
 428     val = -1;
 429     mpcim(&val, MPn1);
 430
 431     if (mpgt(MPx, MP1) || mplt(MPx, MPn1)) {
 432         doerr(_("Error"));
 433         mpstr(MP0, MPretval);
 434     } else if (mpeq(MPx, MP0)) {
 435         val = 2;
 436         mpdivi(MPpi, &val, MPretval);
 437     } else if (mpeq(MPx, MP1)) {
 438         mpstr(MP0, MPretval);
 439     } else if (mpeq(MPx, MPn1)) {
 440         mpstr(MPpi, MPretval);
 441     } else {
 442         mpmul(MPx, MPx, MP2);
 443         mpsub(MP1, MP2, MP2);
 444         mpsqrt(MP2, MP2);
 445         mpdiv(MP2, MPx, MP2);
 446         mpatan(MP2, MPy);
 447         if (mpgt(MPx, MP0)) {
 448             mpstr(MPy, MPretval);
 449         } else {
 450             mpadd(MPy, MPpi, MPretval);
 451         }
 452     }
 453 }
 454
 455
 456 /*  MP precision hyperbolic arc cosine.
 457  *
 458  *  1. If (x < 1.0) then report DOMAIN error and return 0.0.
 459  *
 460  *  2. acosh(x) = log(x + sqrt(x**2 - 1))
 461  */
 462
 463 static void
 464 mpacosh(int *MPx, int *MPretval)
 465 {
 466     int MP1[MP_SIZE], val;
 467
 468     val = 1;
 469     mpcim(&val, MP1);
 470     if (mplt(MPx, MP1)) {
 471         doerr(_("Error"));
 472         val = 0;
 473         mpcim(&val, MPretval);
 474     } else {
 475         mpmul(MPx, MPx, MP1);
 476         val = -1;
 477         mpaddi(MP1, &val, MP1);
 478         mpsqrt(MP1, MP1);
 479         mpadd(MPx, MP1, MP1);
 480         mpln(MP1, MPretval);
 481     }
 482 }
 483
 484
 485 /*  MP precision hyperbolic arc sine.
 486  *
 487  *  1. asinh(x) = log(x + sqrt(x**2 + 1))
 488  */
 489
 490 static void
 491 mpasinh(int *MPx, int *MPretval)
 492 {
 493     int MP1[MP_SIZE], val;
 494
 495     mpmul(MPx, MPx, MP1);
 496     val = 1;
 497     mpaddi(MP1, &val, MP1);
 498     mpsqrt(MP1, MP1);
 499     mpadd(MPx, MP1, MP1);
 500     mpln(MP1, MPretval);
 501 }
 502
 503
 504 /*  MP precision hyperbolic arc tangent.
 505  *
 506  *  1. If (x <= -1.0 or x >= 1.0) then report a DOMAIn error and return 0.0.
 507  *
 508  *  2. atanh(x) = 0.5 * log((1 + x) / (1 - x))
 509  */
 510
 511 static void
 512 mpatanh(int *MPx, int *MPretval)
 513 {
 514     int MP0[MP_SIZE], MP1[MP_SIZE], MP2[MP_SIZE];
 515     int MP3[MP_SIZE], MPn1[MP_SIZE];
 516     int val;
 517
 518     val = 0;
 519     mpcim(&val, MP0);
 520     val = 1;
 521     mpcim(&val, MP1);
 522     val = -1;
 523     mpcim(&val, MPn1);
 524
 525     if (mpge(MPx, MP1) || mple(MPx, MPn1)) {
 526         doerr(_("Error"));
 527         mpstr(MP0, MPretval);
 528     } else {
 529         mpadd(MP1, MPx, MP2);
 530         mpsub(MP1, MPx, MP3);
 531         mpdiv(MP2, MP3, MP3);
 532         mpln(MP3, MP3);
 533         MPstr_to_num("0.5", DEC, MP1);
 534         mpmul(MP1, MP3, MPretval);
 535     }
 536 }
 537
 538
 539 /*  MP precision common log.
 540  *
 541  *  1. log10(x) = log10(e) * log(x)
 542  */
 543
 544 void
 545 mplog10(int *MPx, int *MPretval)
 546 {
 547     int MP1[MP_SIZE], MP2[MP_SIZE], n;
 548
 549     n = 10;
 550     mpcim(&n, MP1);
 551     mpln(MP1, MP1);
 552     mpln(MPx, MP2);
 553     mpdiv(MP2, MP1, MPretval);
 554 }
 555
 556
 557 void
 558 calc_ctrm(int t[MP_SIZE])
 559 {
 560
 561 /*  Cterm - MEM0 = int (periodic interest rate).
 562  *          MEM1 = fv  (future value).
 563  *          MEM2 = pv  (present value).
 564  *
 565  *          RESULT = log(MEM1 / MEM2) / log(1 + MEM0)
 566  */
 567
 568     int val;
 569     int MP1[MP_SIZE], MP2[MP_SIZE], MP3[MP_SIZE], MP4[MP_SIZE];
 570
 571     mpdiv(v->MPmvals[1], v->MPmvals[2], MP1);
 572     mpln(MP1, MP2);
 573     val = 1;
 574     mpaddi(v->MPmvals[0], &val, MP3);
 575     mpln(MP3, MP4);
 576     mpdiv(MP2, MP4, t);
 577 }
 578
 579
 580 void
 581 calc_ddb(int t[MP_SIZE])
 582 {
 583
 584 /*  Ddb   - MEM0 = cost    (amount paid for asset).
 585  *          MEM1 = salvage (value of asset at end of its life).
 586  *          MEM2 = life    (useful life of the asset).
 587  *          MEM3 = period  (time period for depreciation allowance).
 588  *
 589  *          bv = 0.0;
 590  *          for (i = 0; i < MEM3; i++)
 591  *            {
 592  *              VAL = ((MEM0 - bv) * 2) / MEM2
 593  *              bv += VAL
 594  *            }
 595  *          RESULT = VAL
 596  */
 597
 598     int i;
 599     int len;
 600     int val;
 601     int MPbv[MP_SIZE], MP1[MP_SIZE], MP2[MP_SIZE];
 602
 603     i = 0;
 604     mpcim(&i, MPbv);
 605     mpcmi(v->MPmvals[3], &len);
 606     for (i = 0; i < len; i++) {
 607         mpsub(v->MPmvals[0], MPbv, MP1);
 608         val = 2;
 609         mpmuli(MP1, &val, MP2);
 610         mpdiv(MP2, v->MPmvals[2], t);
 611         mpstr(MPbv, MP1);
 612         mpadd(MP1, t, MPbv); /* TODO: why result is MPbv, for next loop? */
 613     }
 614 }
 615
 616
 617 void
 618 calc_fv(int t[MP_SIZE])
 619 {
 620
 621 /*  Fv    - MEM0 = pmt (periodic payment).
 622  *          MEM1 = int (periodic interest rate).
 623  *          MEM2 = n   (number of periods).
 624  *
 625  *          RESULT = MEM0 * (pow(1 + MEM1, MEM2) - 1) / MEM1
 626  */
 627
 628     int val;
 629     int MP1[MP_SIZE], MP2[MP_SIZE], MP3[MP_SIZE], MP4[MP_SIZE];
 630
 631     val = 1;
 632     mpaddi(v->MPmvals[1], &val, MP1);
 633     mppwr2(MP1, v->MPmvals[2], MP2);
 634     val = -1;
 635     mpaddi(MP2, &val, MP3);
 636     mpmul(v->MPmvals[0], MP3, MP4);
 637     mpdiv(MP4, v->MPmvals[1], t);
 638 }
 639
 640
 641 void
 642 calc_pmt(int t[MP_SIZE])
 643 {
 644
 645 /*  Pmt   - MEM0 = prin (principal).
 646  *          MEM1 = int  (periodic interest rate).
 647  *          MEM2 = n    (term).
 648  *
 649  *          RESULT = MEM0 * (MEM1 / (1 - pow(MEM1 + 1, -1 * MEM2)))
 650  */
 651
 652     int val;
 653     int MP1[MP_SIZE], MP2[MP_SIZE], MP3[MP_SIZE], MP4[MP_SIZE];
 654
 655     val = 1;
 656     mpaddi(v->MPmvals[1], &val, MP1);
 657     val = -1;
 658     mpmuli(v->MPmvals[2], &val, MP2);
 659     mppwr2(MP1, MP2, MP3);
 660     val = -1;
 661     mpmuli(MP3, &val, MP4);
 662     val = 1;
 663     mpaddi(MP4, &val, MP1);
 664     mpdiv(v->MPmvals[1], MP1, MP2);
 665     mpmul(v->MPmvals[0], MP2, t);
 666 }
 667
 668
 669 void
 670 calc_pv(int t[MP_SIZE])
 671 {
 672
 673 /*  Pv    - MEM0 = pmt (periodic payment).
 674  *          MEM1 = int (periodic interest rate).
 675  *          MEM2 = n   (term).
 676  *
 677  *          RESULT = MEM0 * (1 - pow(1 + MEM1, -1 * MEM2)) / MEM1
 678  */
 679
 680     int val;
 681     int MP1[MP_SIZE], MP2[MP_SIZE], MP3[MP_SIZE], MP4[MP_SIZE];
 682
 683     val = 1;
 684     mpaddi(v->MPmvals[1], &val, MP1);
 685     val = -1;
 686     mpmuli(v->MPmvals[2], &val, MP2);
 687     mppwr2(MP1, MP2, MP3);
 688     val = -1;
 689     mpmuli(MP3, &val, MP4);
 690     val = 1;
 691     mpaddi(MP4, &val, MP1);
 692     mpdiv(MP1, v->MPmvals[1], MP2);
 693     mpmul(v->MPmvals[0], MP2, t);
 694 }
 695
 696
 697 void
 698 calc_rate(int t[MP_SIZE])
 699 {
 700
 701 /*  Rate  - MEM0 = fv (future value).
 702  *          MEM1 = pv (present value).
 703  *          MEM2 = n  (term).
 704  *
 705  *          RESULT = pow(MEM0 / MEM1, 1 / MEM2) - 1
 706  */
 707
 708     int val;
 709     int MP1[MP_SIZE], MP2[MP_SIZE], MP3[MP_SIZE], MP4[MP_SIZE];
 710
 711     mpdiv(v->MPmvals[0], v->MPmvals[1], MP1);
 712     val = 1;
 713     mpcim(&val, MP2);
 714     mpdiv(MP2, v->MPmvals[2], MP3);
 715     mppwr2(MP1, MP3, MP4);
 716     val = -1;
 717     mpaddi(MP4, &val, t);
 718 }
 719
 720
 721 void
 722 calc_sln(int t[MP_SIZE])
 723 {
 724
 725 /*  Sln   - MEM0 = cost    (cost of the asset).
 726  *          MEM1 = salvage (salvage value of the asset).
 727  *          MEM2 = life    (useful life of the asset).
 728  *
 729  *          RESULT = (MEM0 - MEM1) / MEM2
 730  */
 731
 732     int MP1[MP_SIZE];
 733
 734     mpsub(v->MPmvals[0], v->MPmvals[1], MP1);
 735     mpdiv(MP1, v->MPmvals[2], t);
 736 }
 737
 738
 739 void
 740 calc_syd(int t[MP_SIZE])
 741 {
 742
 743 /*  Syd   - MEM0 = cost    (cost of the asset).
 744  *          MEM1 = salvage (salvage value of the asset).
 745  *          MEM2 = life    (useful life of the asset).
 746  *          MEM3 = period  (period for which depreciation is computed).
 747  *
 748  *          RESULT = ((MEM0 - MEM1) * (MEM2 - MEM3 + 1)) /
 749  *                   (MEM2 * (MEM2 + 1) / 2)
 750  */
 751
 752     int val;
 753     int MP1[MP_SIZE], MP2[MP_SIZE], MP3[MP_SIZE], MP4[MP_SIZE];
 754
 755     mpsub(v->MPmvals[2], v->MPmvals[3], MP2);
 756     val = 1;
 757     mpaddi(MP2, &val, MP3);
 758     mpaddi(v->MPmvals[2], &val, MP2);
 759     mpmul(v->MPmvals[2], MP2, MP4);
 760     val = 2;
 761     mpcim(&val, MP2);
 762     mpdiv(MP4, MP2, MP1);
 763     mpdiv(MP3, MP1, MP2);
 764     mpsub(v->MPmvals[0], v->MPmvals[1], MP1);
 765     mpmul(MP1, MP2, t);
 766 }
 767
 768
 769 void
 770 calc_term(int t[MP_SIZE])
 771 {
 772
 773 /*  Term  - MEM0 = pmt (periodic payment).
 774  *          MEM1 = fv  (future value).
 775  *          MEM2 = int (periodic interest rate).
 776  *
 777  *          RESULT = log(1 + (MEM1 * MEM2 / MEM0)) / log(1 + MEM2)
 778  */
 779
 780     int val;
 781     int MP1[MP_SIZE], MP2[MP_SIZE], MP3[MP_SIZE], MP4[MP_SIZE];
 782
 783     val = 1;
 784     mpaddi(v->MPmvals[2], &val, MP1);
 785     mpln(MP1, MP2);
 786     mpmul(v->MPmvals[1], v->MPmvals[2], MP1);
 787     mpdiv(MP1, v->MPmvals[0], MP3);
 788     val = 1;
 789     mpaddi(MP3, &val, MP4);
 790     mpln(MP4, MP1);
 791     mpdiv(MP1, MP2, t);
 792 }
 793
 794
 795 void
 796 calc_rshift(int s[MP_SIZE], int t[MP_SIZE], int times, enum shiftd dir)
 797 {
 798   /* Implementation derived from old code.
 799    * Using BOOLEAN is strange at least. Assumed that
 800    * boolean means BINARY representation
 801    */
 802
 803     double dval;
 804     BOOLEAN temp;
 805     mpcmd(s, &dval);
 806     temp = ibool(dval);
 807
 808     /* There is a reason to do shift like this. Reason is that
 809      * processors define shift only in a certain range. i386 uses only 5
 810      * bits to describe shiftable amount. So, shift 32 times gives original
 811      * number. That can cause very strange results (and bugs).
 812      */
 813
 814     assert(times >= 0);
 815
 816     while (times--) {
 817         temp = (dir == right) ? temp >> 1 : temp << 1;
 818     }
 819
 820     dval = setbool(temp);
 821     mpcdm(&dval, t);
 822 }
 823
 824
 825 int
 826 is_integer(int MPnum[MP_SIZE])
 827 {
 828     int MP1[MP_SIZE];
 829     mpcmim(MPnum, MP1);
 830     return mpeq(MPnum, MP1);
 831 }
 832
 833 int
 834 is_natural(int MPnum[MP_SIZE])
 835 {
 836     int MP1[MP_SIZE];
 837     if (!is_integer(MPnum)) {
 838         return 0;
 839     }
 840     mpabs(MPnum, MP1);
 841     return mpeq(MPnum, MP1);
 842 }
 843
 844 void
 845 calc_epowy(int s[MP_SIZE], int t[MP_SIZE])
 846 {
 847     int MP1[MP_SIZE];
 848
 849     mpstr(s, MP1);
 850     mpexp(MP1, t);
 851 }