gcalctool/mpmath.c

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