home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Professional / OS2PRO194.ISO / os2 / prgramer / bigcalc / bigmath.c < prev    next >
Encoding:
C/C++ Source or Header  |  1989-02-04  |  87.6 KB  |  2,985 lines
  1. /*
  2.  *    **************************************************
  3.  *    *                                                *
  4.  *    *                   BIGMATH.C                    *
  5.  *    *                                                *
  6.  *    *          Extended Precision Calculator         *
  7.  *    *                                                *
  8.  *    *        Extended Precision Math Routines        *
  9.  *    *                                                *
  10.  *    *              Version 4.3 02-04-89              *
  11.  *    *                                                *
  12.  *    *              Judson D. McClendon               *
  13.  *    *              329 37th Court N.E.               *
  14.  *    *              Birmingham, AL 35215              *
  15.  *    *                 205-853-8440                   *
  16.  *    *            Compuserve [74415,1003]             *
  17.  *    *                                                *
  18.  *    **************************************************
  19.  */
  20.  
  21.  
  22.  
  23. /*
  24.  *    **************************************************
  25.  *    *                                                *
  26.  *    *                   Includes                     *
  27.  *    *                                                *
  28.  *    **************************************************
  29.  */
  30.  
  31. #include <ctype.h>
  32. #include <stdio.h>
  33. #include <stdlib.h>
  34. #include <string.h>
  35. #include "bigcalc.h"
  36. #include "biggvar.h"
  37.  
  38.  
  39.  
  40. /*
  41.  *    **************************************************
  42.  *    *                                                *
  43.  *    *            Extended Math Routines              *
  44.  *    *                                                *
  45.  *    **************************************************
  46.  */
  47.  
  48.  
  49.  
  50. /*
  51.  *    **************************************************
  52.  *    *                                                *
  53.  *    *                 Extended Add                   *
  54.  *    *                                                *
  55.  *    *         work[2] = work[1] + work[0]            *
  56.  *    *                                                *
  57.  *    **************************************************
  58.  */
  59.  
  60. extern int ExtendedAdd(void)
  61.  
  62. {
  63.  
  64.    int a, b;                       /* operand indexes           */
  65.    long shift;                     /* digit allignment shift    */
  66.    WORKDIGIT *ax, *ar, *cx, *ct;   /* ar = a-right, ct = c-temp */
  67.  
  68.  
  69.          /* Special Cases */
  70.  
  71.    if (! work[0].digits) {             /* work[0] zero, return work[1] */
  72.       MoveWorkWork(1, 2);
  73.       return(TRUE);
  74.       }
  75.  
  76.    if (! work[1].digits) {             /* work[1] zero, return work[0] */
  77.       MoveWorkWork(0, 2);
  78.       return(TRUE);
  79.       }
  80.  
  81.    if (work[0].sign != work[1].sign) {       /* If signs are different, */
  82.       work[0].sign = FlipSign(work[0].sign); /*  reverse sign           */
  83.       return(ExtendedSubtract());            /*  and subtract           */
  84.       }
  85.  
  86.    if (work[0].exp <= work[1].exp) {   /* a->smallest, b->largest */
  87.       a = 0;
  88.       b = 1;
  89.       }
  90.    else {
  91.       a = 1;
  92.       b = 0;
  93.       }
  94.  
  95.    shift = work[b].exp - work[a].exp;
  96.    if (shift > (long)compprec) {       /* a operand is insignificant */
  97.       MoveWorkWork(b, 2);              /*  return b operand          */
  98.       return(TRUE);
  99.       }
  100.  
  101.                                        /* work[2] = b shifted 1 digit */
  102.    work[2].man[0] = 0;
  103.    memcpy(&work[2].man[1], work[b].man, (work[b].digits * sizeof(WORKDIGIT)));
  104.    memset(&(work[2].man[work[b].digits + 1]), 0, ((workprec - work[b].digits - 1) * sizeof(WORKDIGIT)));
  105.  
  106.    work[2].exp    = work[b].exp + 1L;
  107.    work[2].sign   = work[b].sign;
  108.  
  109.    if ((work[a].digits + shift + 1) > (work[b].digits + 1))
  110.       work[2].digits = work[a].digits + (int)shift + 1;
  111.    else
  112.       work[2].digits = work[b].digits + 1;
  113.  
  114.    ax = &work[a].man[0];
  115.    ar = &work[a].man[work[a].digits - 1];
  116.    cx = &work[2].man[shift + 1];
  117.  
  118.  
  119.    /*
  120.     *    **************************************************
  121.     *    *                                                *
  122.     *    *  ADDITION: By Shifted Column Add & Carry       *
  123.     *    *                                                *
  124.     *    **************************************************
  125.     */
  126.  
  127.    do {
  128.       if ((*(cx++) += *(ax++)) > 9) {
  129.          ct = cx - 1;                  /* Carry */
  130.          do {
  131.             *ct -= 10;
  132.             } while (++(*(--ct)) > 9);
  133.          }
  134.  
  135.       } while (ax <= ar);
  136.  
  137.    return(Normalize(2));
  138.  
  139. }
  140.  
  141.  
  142.  
  143.  
  144. /*
  145.  *    **************************************************
  146.  *    *                                                *
  147.  *    *               Extended Subtract                *
  148.  *    *                                                *
  149.  *    *         work[2] = work[1] - work[0]            *
  150.  *    *                                                *
  151.  *    **************************************************
  152.  */
  153.  
  154. extern int ExtendedSubtract(void)
  155.  
  156. {
  157.  
  158.    int a, b;                       /* operand indexes               */
  159.    int digits, result;             /* work variables                */
  160.    long shift;                     /* digit allignment shift        */
  161.    WORKDIGIT *al, *ax, *cx, *ct;   /* aleft, aindex, cindex, ctempx */
  162.  
  163.  
  164.          /* Special Cases */
  165.  
  166.    if (! work[0].digits) {             /* work[0] zero, return work[1] */
  167.       MoveWorkWork(1, 2);
  168.       return(TRUE);
  169.       }
  170.  
  171.    if (! work[1].digits) {             /* work[1] zero, return -work[0] */
  172.       MoveWorkWork(0, 2);
  173.       work[2].sign = FlipSign(work[2].sign);
  174.       return(TRUE);
  175.       }
  176.  
  177.    if (work[0].sign != work[1].sign) {       /* If signs different, */
  178.       work[0].sign = FlipSign(work[0].sign); /*  reverse sign       */
  179.       return(ExtendedAdd());                 /*  and add            */
  180.       }
  181.  
  182.    if (work[0].exp < work[1].exp) {    /* compare work[0] : work[1] */
  183.  
  184.       a = 0;                           /* magnitude work[0] < work[1] */
  185.       b = 1;
  186.       }
  187.  
  188.    else if (work[0].exp > work[1].exp) {
  189.  
  190.       a = 1;                           /* Magnitude work[1] > work[0] */
  191.       b = 0;
  192.       work[b].sign = FlipSign(work[b].sign);
  193.       }
  194.  
  195.    else {                                    /* Can't tell by exponents, */
  196.                                              /*  so compare mantissae    */
  197.       if (work[0].digits >= work[1].digits)
  198.          digits = work[0].digits;
  199.       else
  200.          digits = work[1].digits;
  201.  
  202.       result = memcmp(work[0].man, work[1].man, (digits * sizeof(WORKDIGIT)));
  203.  
  204.       if (result < 0) {
  205.          a = 0;                        /* work[0] < work[1] */
  206.          b = 1;
  207.          }
  208.  
  209.       else if (result > 0) {
  210.          a = 1;                        /* work[0] > work[1] */
  211.          b = 0;
  212.          work[b].sign = FlipSign(work[b].sign);
  213.          }
  214.  
  215.       else {
  216.          ClearWork(2);                 /* work[0] == work[1], return 0 */
  217.          return(TRUE);
  218.          }
  219.  
  220.       }  /* End compare */
  221.  
  222.    MoveWorkWork(b, 2);                 /* Put larger number in work[2] */
  223.                                        /*  so borrow always successful */
  224.    shift = work[b].exp - work[a].exp;
  225.    if (shift > (long)compprec) {       /* a operand is insignificant   */
  226.       return(TRUE);                    /*  return larger nbr as answer */
  227.       }
  228.  
  229.    if ((work[a].digits + (int)shift) > work[b].digits)
  230.       work[2].digits = work[a].digits + (int)shift;
  231.    else
  232.       work[2].digits = work[b].digits;
  233.  
  234.    al = &work[a].man[0];
  235.    ax = &work[a].man[work[a].digits - 1];
  236.    cx = &work[2].man[(int)shift + work[a].digits - 1];
  237.  
  238.  
  239.    /*
  240.     *    **************************************************
  241.     *    *                                                *
  242.     *    *  SUBTRACTION: By Shifted Column Sub & Borrow   *
  243.     *    *                                                *
  244.     *    **************************************************
  245.     */
  246.  
  247.    do {
  248.       if ((*(cx--) -= *(ax--)) < 0) {
  249.          ct = cx;                         /* Borrow */
  250.          do {
  251.             *(ct + 1) += 10;              /* Borrow always successful */
  252.             } while (--(*(ct--)) < 0);    /*  because a < c           */
  253.          }
  254.  
  255.       } while (ax >= al);
  256.  
  257.    return(Normalize(2));
  258.  
  259. }
  260.  
  261.  
  262.  
  263.  
  264. /*
  265.  *    **************************************************
  266.  *    *                                                *
  267.  *    *               Extended Multiply                *
  268.  *    *                                                *
  269.  *    *         work[2] = work[1] * work[0]            *
  270.  *    *                                                *
  271.  *    **************************************************
  272.  */
  273.  
  274. extern int ExtendedMultiply(void)
  275.  
  276. {
  277.  
  278.    WORKDIGIT *ax, *bx, *cx;        /* indexes          */
  279.    WORKDIGIT *al, *bl, *cl;        /* left pointers    */
  280.    WORKDIGIT *ar, *br, *cr;        /* right pointers   */
  281.    WORKDIGIT ad;                   /* a digit          */
  282.    long exponent;                  /* work exponent    */
  283.  
  284.  
  285.          /* Special Cases */
  286.  
  287.    ClearWork(2);
  288.  
  289.    if ( (! work[0].digits)             /* If either is zero, */
  290.          ||
  291.         (! work[1].digits) ) {
  292.       return(TRUE);                    /*  return zero       */
  293.       }
  294.  
  295.    exponent = work[0].exp + work[1].exp;
  296.  
  297.    if (exponent < MINEXP + 2L) {       /* If underflow, */
  298.       return(TRUE);                    /*  return zero  */
  299.       }
  300.  
  301.    if (exponent > (MAXEXP + 2L)) {
  302.       Overflow();
  303.       return(FALSE);
  304.       }
  305.  
  306.    work[2].exp = exponent;
  307.  
  308.    if (work[0].sign == work[1].sign)
  309.       work[2].sign = '+';
  310.    else
  311.       work[2].sign = '-';
  312.  
  313.    if (work[0].digits <= work[1].digits) {
  314.       work[0].digits = (work[0].digits + 2) / 3;   /* Round up to mod 3  */
  315.       work[0].digits *= 3;
  316.       al = &work[0].man[0];                        /* Multiplier left    */
  317.       ar = &work[0].man[work[0].digits - 1];       /* Multiplier right   */
  318.       bl = &work[1].man[0];                        /* Multiplicand left  */
  319.       br = &work[1].man[work[1].digits - 1];       /* Multiplicand right */
  320.       }
  321.    else {
  322.       work[1].digits = (work[1].digits + 2) / 3;   /* Round up to mod 3  */
  323.       work[1].digits *= 3;
  324.       al = &work[1].man[0];                        /* Multiplier left    */
  325.       ar = &work[1].man[work[1].digits - 1];       /* Multiplier right   */
  326.       bl = &work[0].man[0];                        /* Multiplicand left  */
  327.       br = &work[0].man[work[0].digits - 1];       /* Multiplicand right */
  328.       }
  329.  
  330.    work[2].digits = work[0].digits + work[1].digits;
  331.    cl = &work[2].man[0];                           /* Product left  */
  332.    cr = &work[2].man[work[2].digits - 1];          /* Product right */
  333.  
  334.  
  335.    /*
  336.     *    **************************************************
  337.     *    *                                                *
  338.     *    *  MULTIPLY: By Repeated 3 Group Multiplication  *
  339.     *    *                                                *
  340.     *    *     Where: a = multiplier in groups of 3       *
  341.     *    *            b = multiplicand                    *
  342.     *    *            c = product                         *
  343.     *    *                                                *
  344.     *    *     bbbbbbbbb     bbbbbbbbb     bbbbbbbbb      *
  345.     *    *                                                *
  346.     *    *           ---        ---        ---            *
  347.     *    *   x aaaaaaaaa   x aaaaaaaaa   x aaaaaaaaa      *
  348.     *    *   -----------   -----------   -----------      *
  349.     *    *     ccccccccc     ccccccccc     ccccccccc      *
  350.     *    *                ddddddddd     ddddddddd         *
  351.     *    *                           eeeeeeeee            *
  352.     *    *                                                *
  353.     *    **************************************************
  354.     */
  355.  
  356.    ax = ar;                                  /* Multiplier (ax) loop */
  357.    do {
  358.       cx = cr;
  359.       cr -= 3;                                  /* Back up each cycle */
  360.  
  361.                                                 /* Multiplier digits */
  362.       ad = *(ax - 2) * 100 + *(ax - 1) * 10 + *ax;
  363.  
  364.       bx = br;                                  /* Multiplicand (bx) loop */
  365.       do {
  366.          if ((*(cx--) += (*bx * ad)) >= 10000) {   /* Partial carry to      */
  367.             *cx += 1000;                           /*  prevent digit overflow */
  368.             *(cx + 1) -= 10000;
  369.             }
  370.          } while (--bx >= bl);                  /* Multiplicand loop end */
  371.  
  372.       if (KeyPressed() == ESCAPE) {             /* User abort */
  373.          Beep();
  374.          return(FALSE);
  375.          }
  376.  
  377.       } while ((ax -= 3) >= al);                /* Multiplier loop end */
  378.  
  379.  
  380.    cx = &work[2].man[work[2].digits];        /* Satisfy all carrys */
  381.    while (--cx > cl) {
  382.       if (*cx > 9) {
  383.          *(cx - 1) += *cx / 10;
  384.          *cx %= 10;
  385.          }
  386.       }
  387.  
  388.    return(Normalize(2));
  389.  
  390. }
  391.  
  392.  
  393.  
  394.  
  395. /*
  396.  *    **************************************************
  397.  *    *                                                *
  398.  *    *                Extended Divide                 *
  399.  *    *                                                *
  400.  *    *         work[2] = work[1] / work[0]            *
  401.  *    *                                                *
  402.  *    **************************************************
  403.  */
  404.  
  405. extern int ExtendedDivide(void)
  406.  
  407. {
  408.  
  409.    WORKDIGIT *ax, *bx, *cx;        /* Indexes               */
  410.    WORKDIGIT *al, *bl;             /* left pointers         */
  411.    WORKDIGIT *ar,      *cr;        /* right pointers        */
  412.    WORKDIGIT *bb, *bt;             /* borrow, temp pointers */
  413.    WORKDIGIT factor, val1, val2;   /* subtract factor       */
  414.    BOOLEAN underflow;              /* underflow flag        */
  415.    long exponent;                  /* long prevents oflow   */
  416.  
  417.  
  418.          /* Special Cases */
  419.  
  420.    if (! work[0].digits) {             /* Divide by zero: error */
  421.       DivideByZero();
  422.       return(FALSE);
  423.       }
  424.  
  425.    ClearWork(2);
  426.  
  427.    if (! work[1].digits) {             /* Divide into zero: result zero */
  428.       return(TRUE);
  429.       }
  430.  
  431.    exponent = work[1].exp - work[0].exp + 1L;
  432.  
  433.    if (exponent < MINEXP + 1L) {       /* If underflow, */
  434.       return(TRUE);                    /*  return zero  */
  435.       }
  436.  
  437.    if (exponent > (MAXEXP + 2L)) {
  438.       Overflow();
  439.       return(FALSE);
  440.       }
  441.  
  442.    work[2].exp = exponent;
  443.  
  444.    if (work[0].sign == work[1].sign)
  445.       work[2].sign = '+';
  446.    else
  447.       work[2].sign = '-';
  448.  
  449.    work[2].digits = compprec + 2;
  450.  
  451.    al = &work[0].man[0];                  /* Divisor left pointer   */
  452.    bl = &work[1].man[0];                  /* Dividend left pointer  */
  453.  
  454.    bb = &work[1].man[0];                  /* Leftmost limit of borrow */
  455.  
  456.    ar = &work[0].man[work[0].digits - 1]; /* Divisor right pointer  */
  457.    cr = &work[2].man[compprec + 1];       /* Quotient right pointer */
  458.  
  459.    cx = &work[2].man[0];                  /* Quotient index pointer */
  460.  
  461.    work[1].man[workprec - 1] = 1;   /* Insignificant digit simplifies scan */
  462.                                     /*  for leftmost digit of dividend     */
  463.  
  464.  
  465.    /*
  466.     *    **************************************************
  467.     *    *                                                *
  468.     *    *  DIVIDE:  By Repeated Weighted Subtraction     *
  469.     *    *                                                *
  470.     *    *     Where: a = divisor * factor                *
  471.     *    *            b = dividend                        *
  472.     *    *            f = factor                          *
  473.     *    *                                                *
  474.     *    *     bbbbbbbb      0bbbbbbb      00bbbbbb       *
  475.     *    *   - aaaaaa       - aaaaaa       - aaaaaa       *
  476.     *    *          |              |              |       *
  477.     *    *          v              v              v       *
  478.     *    *          f             ff            fff       *
  479.     *    *                                                *
  480.     *    **************************************************
  481.     */
  482.  
  483.  
  484.    do {                    /* Divide by repeated subtraction */
  485.       ax = al;
  486.       bx = bl;
  487.                               /* Calculate subtract factor */
  488.       val1 = *bx * 100 + *(bx + 1) * 10 + *(bx + 2);
  489.       if (bb < bx)
  490.          val1 += *bb * 1000;
  491.       val2 = *ax * 100 + *(ax + 1) * 10 + *(ax + 2);
  492.       factor = val1 / val2;
  493.  
  494.       if (factor)
  495.          do {                                /* Subtraction Cycle */
  496.             underflow = FALSE;
  497.             do {
  498.                if ((*bx -= (*ax * factor)) < 0) {
  499.                   bt = bx;
  500.  
  501.                   do {                             /* Begin borrow */
  502.                      if (bt > bb) {
  503.                         *(bt - 1) -= (int)((9 - *bt) / 10);
  504.                         *bt += ((int)((9 - *bt) / 10)) * 10 ;
  505.                         }
  506.                      else {                           /* Borrow failed */
  507.  
  508.                         while (ax >= al) {               /* Undo column subs */
  509.                            if ((*bx += (*ax * factor)) > 9) {
  510.                               bt = bx;                      /* Undo prev borrow */
  511.                               do {
  512.                                  *(bt - 1) += *bt / 10;
  513.                                  *bt %= 10;
  514.                                  } while (*(--bt) > 9);
  515.                               }                             /* Undo prev borrow end */
  516.                            ax--;
  517.                            bx--;
  518.                            }                             /* Undo column subs end */
  519.                         underflow = TRUE;
  520.                         break;
  521.                         }                             /* Borrow failed end */
  522.                      } while (*(--bt) < 0);
  523.                   }                                /* Borrow end */
  524.                ax++;
  525.                bx++;
  526.                } while ((! underflow) && (ax <= ar));
  527.  
  528.             factor -= underflow;                /* Subtracts if TRUE (== 1) */
  529.  
  530.             } while (underflow && factor);   /* Subtract cycle end */
  531.  
  532.       *cx = factor;           /* Result of subtraction */
  533.  
  534.       while (! *bb)           /* Find leftmost non zero is      */
  535.          bb++;                /*  always successful (see above) */
  536.  
  537.       if (bb > bl) {          /* If 1st non zero is right of current shift */
  538.          cx += (bb - bl);     /*  then shift right to it,                  */
  539.          bl = bb;
  540.          }
  541.       else {                  /*  else */
  542.          cx++;                /*  shift one column to the right */
  543.          bl++;
  544.          }
  545.  
  546.       if (KeyPressed() == ESCAPE) { /* User abort */
  547.          Beep();
  548.          return(FALSE);
  549.          }
  550.  
  551.       } while (cx <= cr);  /* Divide loop end */
  552.  
  553.  
  554.    return(Normalize(2));
  555.  
  556. }
  557.  
  558.  
  559.  
  560.  
  561. /*
  562.  *    **************************************************
  563.  *    *                                                *
  564.  *    *             Extended Square Root               *
  565.  *    *                                                *
  566.  *    *            work[2] = Sqrt(work[0])             *
  567.  *    *                                                *
  568.  *    *              (Uses 2 temp areas)               *
  569.  *    *                                                *
  570.  *    **************************************************
  571.  */
  572.  
  573. extern int ExtendedSquareRoot(void)
  574.  
  575. {
  576.  
  577.    long nl, apxl, lapxl;
  578.    int i, j, result;
  579.    int oldprec, cmpsize;
  580.    COMPTYPE *temp;
  581.    COMPTYPE *arg, *apx;
  582.  
  583.  
  584.          /* Special Cases */
  585.  
  586.    if (work[0].sign == '-') {                /* Negative: error */
  587.       NegativeArgument();
  588.       return(FALSE);
  589.       }
  590.  
  591.    if (! work[0].digits) {                   /* Zero: return zero */
  592.       ClearWork(2);
  593.       return(TRUE);
  594.       }
  595.  
  596.  
  597.          /* Set up temp registers */
  598.  
  599.    if ((temp = GETCOMPTEMP(2)) == NULL) {
  600.       MemoryError();
  601.       return(FALSE);
  602.       }
  603.    arg = temp;
  604.    apx = temp + 1;
  605.  
  606.    MoveWorkTemp(0, arg);                     /* Argument -> *arg */
  607.  
  608.  
  609.          /* Prepare first approximation */
  610.  
  611.    nl = 0L;                                  /* Extract 1st digits of N */
  612.    j = 8 + (int)(arg->exp & 0x01L);  /* 9 if odd, 8 if even */
  613.    for (i = 0; i < j; i++) {
  614.       nl = nl * 10L + (long)arg->man[i];
  615.       }
  616.  
  617.    if (j == 8L)                              /* Compute sqrt of digits */
  618.       apxl = 7071L;
  619.    else
  620.       apxl = 22360L;
  621.  
  622.    for (i = 1; i <= 6; i++) {
  623.       lapxl = apxl;
  624.       apxl = ((lapxl + (nl / lapxl)) >> 1);
  625.       }
  626.    SetWorkInteger(2, apxl);                  /* Put 1st appx in work[2] */
  627.  
  628.    if (arg->exp > 0L)                        /* Calc 1st approx exponent */
  629.       work[2].exp = (arg->exp + 1L) / 2L;
  630.    else
  631.       work[2].exp = arg->exp / 2L;
  632.    work[2].sign = '+';
  633.  
  634.  
  635.    /*
  636.     *    **************************************************
  637.     *    *                                                *
  638.     *    *  SQUARE ROOT:  Using Newton's Approximation    *
  639.     *    *                                                *
  640.     *    *     apx    = ( apx  + ( X / apx ) ) / 2        *
  641.     *    *        i+1        i            i               *
  642.     *    *                                                *
  643.     *    *                                                *
  644.     *    *  *arg = X                                      *
  645.     *    *                                                *
  646.     *    *  *apx = apx                                    *
  647.     *    *            i                                   *
  648.     *    *                                                *
  649.     *    **************************************************
  650.     */
  651.  
  652.    oldprec = normprec;                 /* Use progressive precision */
  653.    if (normprec > 8) {
  654.       normprec = 8;
  655.       compprec = COMPPREC;
  656.       workprec = WORKPREC;
  657.       }
  658.    cmpsize = compprec * sizeof(WORKDIGIT);
  659.  
  660.  
  661.    do {
  662.  
  663.       MoveWorkTemp(2, apx);               /* Save prev appx */
  664.  
  665.       MoveTempWork(arg, 1);               /* Divide arg */
  666.  
  667.       MoveWorkWork(2, 0);                 /*  by prev appx */
  668.  
  669.       if (! ExtendedDivide()) {
  670.          result = FALSE;
  671.          break;
  672.          }
  673.  
  674.  
  675.       MoveWorkWork(2, 1);                 /* Add result to prev appx */
  676.                                           /*  still in work[0]       */
  677.       if (! ExtendedAdd()) {
  678.          result = FALSE;
  679.          break;
  680.          }
  681.  
  682.  
  683.       MoveWorkWork(2, 1);                 /* Multiply result */
  684.  
  685.       ClearWork(2);                       /*  by .5 */
  686.       work[2].exp    = work[1].exp;
  687.       work[2].sign   = work[1].sign;
  688.       work[2].digits = work[1].digits + 1;
  689.       for (i = work[1].digits - 1; i >= 0; i--) {
  690.          if ((work[2].man[i+1] += work[1].man[i] * 5) > 9) {
  691.             work[2].man[i] = work[2].man[i+1] / 10;
  692.             work[2].man[i+1] %= 10;
  693.             }
  694.          }
  695.       if ( (! work[2].man[0])
  696.             ||
  697.            (work[2].digits > compprec) ) {
  698.          Normalize(2);
  699.          }
  700.  
  701.                                           /* Test for completion */
  702.       if (compprec >= arg->digits) {
  703.          if (! memcmp(work[2].man, apx->man, cmpsize) ) {
  704.             result = TRUE;
  705.             break;
  706.             }
  707.          }
  708.  
  709.       if (normprec < oldprec) {           /* Increase precision */
  710.          normprec *= 2;
  711.          if (normprec > oldprec)
  712.             normprec = oldprec;
  713.          compprec = COMPPREC;
  714.          workprec = WORKPREC;
  715.          cmpsize = compprec * sizeof(WORKDIGIT);
  716.          }
  717.  
  718.       } while (TRUE);               /* Loop until precision gained */
  719.  
  720.  
  721.    if (normprec < oldprec) {        /* Reset precision if needed */
  722.       normprec = oldprec;
  723.       compprec = COMPPREC;
  724.       workprec = WORKPREC;
  725.       }
  726.  
  727.    free(temp);
  728.  
  729.    return(result);
  730.  
  731. }
  732.  
  733.  
  734.  
  735.  
  736. /*
  737.  *    **************************************************
  738.  *    *                                                *
  739.  *    *                Extended Power                  *
  740.  *    *                                                *
  741.  *    *          work[2] = work[1] ^ work[0]           *
  742.  *    *                                                *
  743.  *    *              (Uses 1 temp area)                *
  744.  *    *                                                *
  745.  *    **************************************************
  746.  */
  747.  
  748. extern int ExtendedPower(void)
  749.  
  750. {
  751.  
  752.    int i;
  753.    unsigned long pow;
  754.    COMPTYPE *temp, *x, *y;
  755.    BOOLEAN recip = FALSE;
  756.  
  757.  
  758.          /* Special Cases */
  759.  
  760.    if (! work[1].digits) {                /* Y zero: return 0 */
  761.       ClearWork(2);
  762.       return(TRUE);
  763.       }
  764.  
  765.    if (! work[0].digits) {                /* X zero: return 1 */
  766.       SetWorkInteger(0, 1L);
  767.       return(TRUE);
  768.       }
  769.  
  770.    if (work[0].exp > MAXEXDIGITS) {          /* Overflow */
  771.       Overflow();
  772.       return(FALSE);
  773.       }
  774.  
  775.  
  776.          /* Loop if X 'small' integer */
  777.  
  778.    if ( (work[0].exp >= work[0].digits)   /* If X integer */
  779.          &&
  780.         (work[0].exp <= MAXEXDIGITS) ) {  /*  and 'small' */
  781.  
  782.  
  783.             /* Set up temp registers */
  784.  
  785.       if ((temp = GETCOMPTEMP(1)) == NULL) {
  786.          MemoryError();
  787.          return(FALSE);
  788.          }
  789.       y = temp;
  790.       MoveWorkTemp(1, y);                 /* Save Y */
  791.  
  792.       if (work[0].sign == '-')            /* If X neg, use reciprocal */
  793.          recip = TRUE;
  794.  
  795.  
  796.       /*
  797.        *    **************************************************
  798.        *    *                                                *
  799.        *    *  POWER: Y^X Using Square-Multiply algorithm:   *
  800.        *    *                                                *
  801.        *    *                                                *
  802.        *    *     Shift X left 1 bit                         *
  803.        *    *                                                *
  804.        *    *     Set rightmost bit of X to 1 (flag bit)     *
  805.        *    *                                                *
  806.        *    *     Shift X left until high order bit is 1     *
  807.        *    *                                                *
  808.        *    *     Shift X left again                         *
  809.        *    *                                                *
  810.        *    *     Set TEMP = Y                               *
  811.        *    *                                                *
  812.        *    *     While X <> (high order bit 1, rest 0)      *
  813.        *    *                                                *
  814.        *    *        Square TEMP                             *
  815.        *    *                                                *
  816.        *    *        If high order bit of X = 1              *
  817.        *    *                                                *
  818.        *    *           Multiply Y times TEMP                *
  819.        *    *                                                *
  820.        *    *        Shift X left 1 bit                      *
  821.        *    *                                                *
  822.        *    *     Result in TEMP is Y^X                      *
  823.        *    *                                                *
  824.        *    *                                                *
  825.        *    *     pow = X                                    *
  826.        *    *                                                *
  827.        *    **************************************************
  828.        */
  829.  
  830.       pow = 0;                            /* X -> pow */
  831.       for (i = 0; i < (int)work[0].exp; i++) {
  832.          pow = pow * 10 + work[0].man[i];
  833.          }
  834.       pow = (pow << 1) + 1;               /* Add flag bit & left justify */
  835.       while (pow < (unsigned long)0x80000000L) {
  836.          pow <<= 1;
  837.          }
  838.       pow <<= 1;
  839.  
  840.       MoveWorkWork(1, 2);                 /* Y -> work[2] (TEMP) */
  841.  
  842.  
  843.       while (pow != (unsigned long)0x80000000L) {
  844.  
  845.          MoveWorkWork(2, 0);              /* Square TEMP */
  846.          MoveWorkWork(2, 1);
  847.          if (! ExtendedMultiply()) {
  848.             free(temp);
  849.             return(FALSE);
  850.             }
  851.  
  852.          if (! work[2].digits) {          /* Test for underflow */
  853.             free(temp);
  854.             if (recip)
  855.                return(FALSE);             /* Divide by zero */
  856.             else
  857.                return(TRUE);
  858.             }
  859.  
  860.          if (pow & (unsigned long)0x80000000L) {
  861.  
  862.             MoveTempWork(y, 0);              /* Multiply by Y */
  863.             MoveWorkWork(2, 1);
  864.             if (! ExtendedMultiply()) {
  865.                free(temp);
  866.                return(FALSE);
  867.                }
  868.             if (! work[2].digits) {          /* Test for underflow */
  869.                free(temp);
  870.                if (recip)
  871.                   return(FALSE);             /* Divide by zero */
  872.                else
  873.                   return(TRUE);
  874.                }
  875.             }
  876.  
  877.          pow <<= 1;                       /* Shift pow left 1 bit */
  878.  
  879.          }
  880.  
  881.       free(temp);
  882.  
  883.       if (recip) {                        /* X neg: use reciprocal */
  884.          MoveWorkWork(2, 0);
  885.          SetWorkInteger(1, 1L);
  886.          if (! ExtendedDivide()) {
  887.             return(FALSE);
  888.             }
  889.          }
  890.  
  891.       return(TRUE);
  892.  
  893.       }  /* if */
  894.  
  895.  
  896.    /*
  897.     *    **************************************************
  898.     *    *                                                *
  899.     *    *  POWER: Using logarithms:                      *
  900.     *    *                                                *
  901.     *    *     Y^X = e^(ln(Y^X))                          *
  902.     *    *                                                *
  903.     *    *         = e^(X * ln(Y))                        *
  904.     *    *                                                *
  905.     *    *                                                *
  906.     *    *     *x = X                                     *
  907.     *    *                                                *
  908.     *    **************************************************
  909.     */
  910.  
  911.  
  912.          /* Set up temp registers */
  913.  
  914.    if ((temp = GETCOMPTEMP(1)) == NULL) {
  915.       MemoryError();
  916.       return(FALSE);
  917.       }
  918.    x = temp;
  919.  
  920.  
  921.          /* Save X */
  922.  
  923.    MoveWorkTemp(0, x);
  924.  
  925.  
  926.          /* Compute ln(Y) */
  927.  
  928.    MoveWorkWork(1, 0);              /* Get Y */
  929.    if (! ExtendedLn()) {
  930.       free(temp);
  931.       return(FALSE);
  932.       }
  933.  
  934.          /* Compute X * ln(Y) */
  935.  
  936.    MoveWorkWork(2, 1);              /* Get ln(Y) */
  937.    MoveTempWork(x, 0);              /* Get X */
  938.    if (! ExtendedMultiply()) {
  939.       free(temp);
  940.       return(FALSE);
  941.       }
  942.  
  943.    free(temp);
  944.  
  945.          /* Compute exp(X * ln(Y) */
  946.  
  947.    MoveWorkWork(2, 0);              /* Get X * ln(Y) */
  948.    if (! ExtendedExpE()) {
  949.       return(FALSE);
  950.       }
  951.  
  952.    return(TRUE);
  953.  
  954. }
  955.  
  956.  
  957.  
  958.  
  959. /*
  960.  *    **************************************************
  961.  *    *                                                *
  962.  *    *                  Sine/Cosine                   *
  963.  *    *                                                *
  964.  *    *             work[2] = Sin(work[0])             *
  965.  *    *                 if scflag = 0                  *
  966.  *    *                                                *
  967.  *    *             work[2] = Cos(work[0])             *
  968.  *    *                 if scflag = 1                  *
  969.  *    *                                                *
  970.  *    *              (Uses 2 temp areas)               *
  971.  *    *                                                *
  972.  *    **************************************************
  973.  */
  974.  
  975. extern int ExtendedSinCos(int scflag)
  976.  
  977. {
  978.  
  979.    int i;
  980.    long n, count = 0;
  981.    int sign = '+';
  982.    COMPTYPE *temp;
  983.    COMPTYPE *x, *xsq, *sum;
  984.    BOOLEAN add = FALSE;
  985.  
  986.  
  987.          /* Special Cases */
  988.  
  989.    if (! work[0].digits) {                /* X zero */
  990.       if (scflag)
  991.          SetWorkInteger(2, 1L);           /*  cos(0) = 1 */
  992.       else
  993.          ClearWork(2);                    /*  sin(0) = 0 */
  994.       return(TRUE);
  995.       }
  996.  
  997.    if (work[0].exp > MAXEXDIGITS) {       /* Overflow */
  998.       Overflow();
  999.       return(FALSE);
  1000.       }
  1001.  
  1002.  
  1003.          /* Set up temp registers */
  1004.  
  1005.    if ((temp = GETCOMPTEMP(2)) == NULL) {
  1006.       MemoryError();
  1007.       return(FALSE);
  1008.       }
  1009.    x = temp;                              /* x and xsq never used together */
  1010.    xsq  = temp;
  1011.    sum  = temp + 1;
  1012.  
  1013.  
  1014.    /*
  1015.     *    **************************************************
  1016.     *    *                                                *
  1017.     *    *  Scale: -pi/4 < X < pi/4  like this:           *
  1018.     *    *                                                *
  1019.     *    *                                                *
  1020.     *    *     X = X - (pi/2 * int(X / pi/2))             *
  1021.     *    *                                                *
  1022.     *    *     if (|X| > pi/4)                            *
  1023.     *    *                                                *
  1024.     *    *        X = X - (pi/2 * sign(X))                *
  1025.     *    *                                                *
  1026.     *    *                                                *
  1027.     *    *     *x = X                                     *
  1028.     *    *                                                *
  1029.     *    **************************************************
  1030.     */
  1031.  
  1032.  
  1033.    if (work[0].exp > 1) {              /* Check magnitude of X */
  1034.  
  1035.       MoveWorkTemp(0, x);              /* Save X */
  1036.  
  1037.       MoveWorkWork(0, 1);              /* Calculate int(X/(pi/2)) */
  1038.       ExtendedRecallHalfPi(0);
  1039.       if (! ExtendedDivide()) {
  1040.          free(temp);
  1041.          return(FALSE);
  1042.          }
  1043.       MoveWorkWork(2, 0);
  1044.       if (! ExtendedIntegerPart()) {
  1045.          free(temp);
  1046.          return(FALSE);
  1047.          }
  1048.       count = 0;                       /* int(X/(pi/2)) -> count */
  1049.       for (i = 0; i < (int)work[0].exp; i++) {
  1050.          count = count * 10 + work[0].man[i];
  1051.          }
  1052.  
  1053.       ExtendedRecallHalfPi(1);         /* Calculate pi/2 * int(X/(pi/2)) */
  1054.       if (! ExtendedMultiply()) {
  1055.          free(temp);
  1056.          return(FALSE);
  1057.          }
  1058.  
  1059.       MoveWorkWork(2, 0);              /* Calculate X - ((pi/2) * int(X/(pi/2))) */
  1060.       MoveTempWork(x, 1);
  1061.       if (! ExtendedSubtract()) {
  1062.          free(temp);
  1063.          return(FALSE);
  1064.          }
  1065.  
  1066.       MoveWorkWork(2, 0);              /* -pi/2 < X < pi/2 */
  1067.  
  1068.       }  /* end pi/2 scale */
  1069.  
  1070.  
  1071.    i = work[0].man[0] * 1000 +         /* Get 1st 4 digits of X */
  1072.        work[0].man[1] * 100 +
  1073.        work[0].man[2] * 10 +
  1074.        work[0].man[3];
  1075.  
  1076.    while ( (work[0].exp > 0)           /* Check magnitude of X */
  1077.            ||
  1078.            ( (work[0].exp == 0)
  1079.               &&
  1080.              (i > 7853) ) ) {
  1081.  
  1082.       MoveWorkWork(0, 1);              /* Calculate X - (pi/2 * sign X) */
  1083.       ExtendedRecallHalfPi(0);
  1084.       work[0].sign = work[1].sign;
  1085.       if (! ExtendedSubtract()) {
  1086.          free(temp);
  1087.          return(FALSE);
  1088.          }
  1089.  
  1090.       count++;
  1091.       MoveWorkWork(2, 0);              /* -pi/4 < X < pi/4 */
  1092.  
  1093.       i = work[0].man[0] * 1000 +      /* Get 1st 4 digits of X */
  1094.           work[0].man[1] * 100 +
  1095.           work[0].man[2] * 10 +
  1096.           work[0].man[3];
  1097.  
  1098.       }  /* end pi/4 scale */
  1099.  
  1100.  
  1101.  
  1102.    /*
  1103.     *    **************************************************
  1104.     *    *                                                *
  1105.     *    *  SIN(X): Using the power series:               *
  1106.     *    *                                                *
  1107.     *    *                                                *
  1108.     *    *                  X^3   X^5   X^7       X^n     *
  1109.     *    *     sin(X) = X - --- + --- - --- + ... ---     *
  1110.     *    *                   3!    5!    7!        n!     *
  1111.     *    *                                                *
  1112.     *    *  COS(X): Using the power series:               *
  1113.     *    *                                                *
  1114.     *    *                                                *
  1115.     *    *                  X^2   X^4   X^6       X^n     *
  1116.     *    *     cos(X) = 1 - --- + --- - --- + ... ---     *
  1117.     *    *                   2!    4!    6!        n!     *
  1118.     *    *                                                *
  1119.     *    *                                                *
  1120.     *    *     *xsq = X^2                                 *
  1121.     *    *                                                *
  1122.     *    *     *sum = Sum of series                       *
  1123.     *    *                                                *
  1124.     *    *     n    = Term number                         *
  1125.     *    *                                                *
  1126.     *    **************************************************
  1127.     */
  1128.  
  1129.  
  1130.  
  1131.    MoveWorkWork(0, 1);              /* Calculate X^2 */
  1132.    if (! ExtendedMultiply()) {
  1133.       free(temp);
  1134.       return(FALSE);
  1135.       }
  1136.    MoveWorkTemp(2, xsq);
  1137.  
  1138.    count += scflag;
  1139.  
  1140.    if (count & 2L) {                /* If count & 2: flip sign */
  1141.       sign = FlipSign(sign);
  1142.       }
  1143.  
  1144.    if (count & 1L) {                /* If count & 1: use cos series */
  1145.       SetWorkInteger(0, 1L);
  1146.       n = 0L;                       /* Even terms    */
  1147.       }
  1148.    else {
  1149.       n = 1L;                       /* Odd terms     */
  1150.       }
  1151.    MoveWorkTemp(0, sum);            /* X -> sum */
  1152.  
  1153.  
  1154.    while (TRUE) {                /* Until term < compprec */
  1155.  
  1156.       n += 2;                       /* n = term number */
  1157.  
  1158.          /* term *= X */
  1159.  
  1160.          /* term still in work[0] */
  1161.  
  1162.       MoveTempWork(xsq, 1);         /* Get xsq */
  1163.       if (! ExtendedMultiply()) {
  1164.          free(temp);
  1165.          return(FALSE);
  1166.          }
  1167.  
  1168.  
  1169.          /* term /= n! */
  1170.  
  1171.       SetWorkInteger(0, (n*(n-1))); /* Get n! */
  1172.       MoveWorkWork(2, 1);           /* Get term */
  1173.       if (! ExtendedDivide()) {
  1174.          free(temp);
  1175.          return(FALSE);
  1176.          }
  1177.  
  1178.  
  1179.          /* Add to/subtract from sum */
  1180.  
  1181.       MoveWorkWork(2, 0);           /* new term */
  1182.       MoveTempWork(sum, 1);         /* Get sum */
  1183.       if (add) {
  1184.          if (! ExtendedAdd()) {
  1185.             free(temp);
  1186.             return(FALSE);
  1187.             }
  1188.          add = FALSE;
  1189.          }
  1190.       else {
  1191.          if (! ExtendedSubtract()) {
  1192.             free(temp);
  1193.             return(FALSE);
  1194.             }
  1195.          add = TRUE;
  1196.          }
  1197.  
  1198.  
  1199.          /* Test for completion */
  1200.  
  1201.       if ( (labs(work[2].exp - work[0].exp) >= compprec)
  1202.             ||
  1203.            (! work[0].digits) ) {
  1204.          break;
  1205.          }
  1206.  
  1207.  
  1208.          /* Save new sum */
  1209.  
  1210.       MoveWorkTemp(2, sum);
  1211.  
  1212.       }  /* while */
  1213.  
  1214.    free(temp);
  1215.  
  1216.    if (sign == '-') {
  1217.       work[2].sign = FlipSign(work[2].sign);
  1218.       }
  1219.  
  1220.    return(TRUE);
  1221.  
  1222. }
  1223.  
  1224.  
  1225.  
  1226.  
  1227. /*
  1228.  *    **************************************************
  1229.  *    *                                                *
  1230.  *    *                   Tangent                      *
  1231.  *    *                                                *
  1232.  *    *             work[2] = Tan(work[0])             *
  1233.  *    *                                                *
  1234.  *    **************************************************
  1235.  */
  1236.  
  1237. extern int ExtendedTan(void)
  1238.  
  1239. {
  1240.  
  1241.    COMPTYPE *temp;
  1242.  
  1243.  
  1244.          /* Special Cases */
  1245.  
  1246.    if (! work[0].digits) {                /* X zero: return zero */
  1247.       ClearWork(2);
  1248.       return(TRUE);
  1249.       }
  1250.  
  1251.  
  1252.    /*
  1253.     *    **************************************************
  1254.     *    *                                                *
  1255.     *    *  TAN(X) Using these relations:                 *
  1256.     *    *                                                *
  1257.     *    *     tan(X) = sin(X) / cos(X)                   *
  1258.     *    *                                                *
  1259.     *    *                                                *
  1260.     *    *     *temp = x or sin(X)                        *
  1261.     *    *                                                *
  1262.     *    **************************************************
  1263.     */
  1264.  
  1265.  
  1266.          /* Set up temp register */
  1267.  
  1268.    if ((temp = GETCOMPTEMP(1)) == NULL) {
  1269.       MemoryError();
  1270.       return(FALSE);
  1271.       }
  1272.    temp = temp;
  1273.  
  1274.    MoveWorkTemp(0, temp);        /* Save X */
  1275.  
  1276.    if (! ExtendedSinCos(0)) {    /* Compute sin(X) */
  1277.       free(temp);
  1278.       return(FALSE);
  1279.       }
  1280.  
  1281.    MoveTempWork(temp, 0);        /* Get X */
  1282.    MoveWorkTemp(2, temp);        /* Save sin(X) */
  1283.  
  1284.    if (! ExtendedSinCos(1)) {    /* Compute cos(X) */
  1285.       free(temp);
  1286.       return(FALSE);
  1287.       }
  1288.  
  1289.    free(temp);
  1290.  
  1291.    MoveWorkWork(2, 0);           /* Compute sin(X) / cos(X) */
  1292.    MoveTempWork(temp, 1);
  1293.    if (! ExtendedDivide()) {
  1294.       return(FALSE);
  1295.       }
  1296.  
  1297.    return(TRUE);
  1298.  
  1299. }
  1300.  
  1301.  
  1302.  
  1303.  
  1304. /*
  1305.  *    **************************************************
  1306.  *    *                                                *
  1307.  *    *                Arc Sine/Cosine                 *
  1308.  *    *                                                *
  1309.  *    *            work[2] = ArcSin(work[0])           *
  1310.  *    *                 if scflag = 0                  *
  1311.  *    *                                                *
  1312.  *    *            work[2] = ArcCos(work[0])           *
  1313.  *    *                 if scflag = 1                  *
  1314.  *    *                                                *
  1315.  *    *               (uses 1 temp area)               *
  1316.  *    *                                                *
  1317.  *    **************************************************
  1318.  */
  1319.  
  1320. extern int ExtendedArcSinCos(int scflag)
  1321.  
  1322. {
  1323.  
  1324.    int sign = '+';
  1325.    COMPTYPE *temp;
  1326.  
  1327.  
  1328.  
  1329.          /* Special cases */
  1330.  
  1331.    if (work[0].sign == '-') {          /* Convert to positive */
  1332.       work[0].sign = '+';
  1333.       sign = '-';
  1334.       }
  1335.  
  1336.    if ( (work[0].exp > 1)                /* Arg > 1: return error */
  1337.          ||
  1338.         ( (work[0].exp == 1)
  1339.            &&
  1340.           ( (work[0].man[0] > 1)
  1341.              ||
  1342.             (work[0].digits > 1) ) ) ) {
  1343.       ArgumentInvalid();
  1344.       return(FALSE);
  1345.       }
  1346.  
  1347.  
  1348.  
  1349. /*
  1350.  *    **************************************************
  1351.  *    *                                                *
  1352.  *    *  Calculate arcsin, and convert if arccos       *
  1353.  *    *                                                *
  1354.  *    **************************************************
  1355.  */
  1356.  
  1357.    if (! work[0].digits) {             /* Zero: result 0 */
  1358.       ClearWork(2);
  1359.       }
  1360.    else if ( (work[0].exp == 1)        /* One: result = pi/2 */
  1361.               &&
  1362.              (work[0].man[0] == 1)
  1363.               &&
  1364.              (work[0].digits == 1) ) {
  1365.       ExtendedRecallHalfPi(2);
  1366.       }
  1367.    else {                              /* Calculate angle */
  1368.       if ((temp = GETCOMPTEMP(1)) == NULL) {
  1369.          MemoryError();
  1370.          return(FALSE);
  1371.          }
  1372.       MoveWorkTemp(0, temp);           /* Save 'sine' */
  1373.  
  1374.       MoveWorkWork(0, 1);              /* Calculate 'cosine' */
  1375.       if (! ExtendedMultiply()) {
  1376.          free(temp);
  1377.          return(FALSE);
  1378.          }
  1379.  
  1380.       MoveWorkWork(2, 0);
  1381.       SetWorkInteger(1, 1L);
  1382.       ExtendedSubtract();
  1383.  
  1384.       MoveWorkWork(2, 0);
  1385.       if (! ExtendedSquareRoot()) {    /* cos = sqrt(1 - sin^2) */
  1386.          free(temp);
  1387.          return(FALSE);
  1388.          }
  1389.  
  1390.       MoveTempWork(temp, 0);           /* 'sine' */
  1391.       MoveWorkWork(2, 1);              /* 'cosine' */
  1392.  
  1393.       free(temp);
  1394.  
  1395.       if (! ResolveAngle()) {
  1396.          return(FALSE);
  1397.          }
  1398.       }
  1399.  
  1400.    work[2].sign = sign;                /* Correct sign */
  1401.  
  1402.    if (scflag) {                       /* Convert to arccos */
  1403.       MoveWorkWork(2, 0);
  1404.       ExtendedRecallHalfPi(1);
  1405.       ExtendedSubtract();
  1406.       }
  1407.  
  1408.    return(TRUE);
  1409.  
  1410. }
  1411.  
  1412.  
  1413.  
  1414.  
  1415. /*
  1416.  *    **************************************************
  1417.  *    *                                                *
  1418.  *    *                  ArcTangent                    *
  1419.  *    *                                                *
  1420.  *    *           work[2] = ArcTan(work[0])            *
  1421.  *    *                                                *
  1422.  *    **************************************************
  1423.  */
  1424.  
  1425. extern int ExtendedArcTan(void)
  1426.  
  1427. {
  1428.  
  1429.    int sign = '+', recip = 0;
  1430.    COMPTYPE *temp;
  1431.  
  1432.  
  1433. /*
  1434.  *    **************************************************
  1435.  *    *                                                *
  1436.  *    *  Convert tan to sin and arcsin, then resolve   *
  1437.  *    *                                                *
  1438.  *    **************************************************
  1439.  */
  1440.  
  1441.  
  1442.          /* Special cases */
  1443.  
  1444.    if (work[0].sign == '-') {          /* Convert to positive */
  1445.       work[0].sign = '+';
  1446.       sign = '-';
  1447.       }
  1448.                            
  1449.    if ( (work[0].exp > 1)              /* Arg > 1: use reciprocal: */
  1450.          ||                            /* atan(1/x) = pi/2 - atan(x) */
  1451.         ( (work[0].exp == 1)
  1452.            &&
  1453.           ( (work[0].man[0] > 1)
  1454.              ||
  1455.             (work[0].digits > 1) ) ) ) {
  1456.       recip = 1;
  1457.       if (! ExtendedReciprocal())
  1458.          return(FALSE);
  1459.       MoveWorkWork(2, 0);
  1460.       }
  1461.  
  1462.  
  1463.          /* Get temp area */
  1464.  
  1465.    if ((temp = GETCOMPTEMP(1)) == NULL) {
  1466.       MemoryError();
  1467.       return(FALSE);
  1468.       }
  1469.    MoveWorkTemp(0, temp);              /* Save x */
  1470.  
  1471.  
  1472.          /* Calculate cos = 1 / sqrt(1 + tan^2) */
  1473.  
  1474.    MoveWorkWork(0, 1);
  1475.    if (! ExtendedMultiply()) {         /* tan^2 */
  1476.       free(temp);
  1477.       return(FALSE);
  1478.       }
  1479.    MoveWorkWork(2, 0);
  1480.    SetWorkInteger(1, 1L);
  1481.    ExtendedAdd();                      /* 1 + tan^2 */
  1482.    MoveWorkWork(2, 0);
  1483.    if (! ExtendedSquareRoot()) {       /* sqrt(1 + tan^2) */
  1484.       free(temp);
  1485.       return(FALSE);
  1486.       }
  1487.    MoveWorkWork(2, 0);
  1488.    SetWorkInteger(1, 1L);
  1489.    if (! ExtendedDivide()) {           /* cos = 1 / sqrt(1 + tan^2) */
  1490.       free(temp);
  1491.       return(FALSE);
  1492.       }
  1493.  
  1494.  
  1495.          /* Calculate sin = cos * tan */
  1496.  
  1497.    MoveWorkWork(2, 1);
  1498.    MoveTempWork(temp, 0);
  1499.    free(temp);
  1500.    if (! ExtendedMultiply())
  1501.       return(FALSE);                   /* cos still in work[1] */
  1502.    MoveWorkWork(2, 0);                 /* sin -> work[0] */
  1503.  
  1504.    if (! ResolveAngle())               /* Resolve angle */
  1505.       return(FALSE);
  1506.  
  1507.    if (recip) {                        /* Resolve reciprocal: */
  1508.       MoveWorkWork(2, 0);              /* atan(1/x) = pi/2 - atan(x) */
  1509.       ExtendedRecallHalfPi(1);
  1510.       ExtendedSubtract();
  1511.       }
  1512.  
  1513.    work[2].sign = sign;                /* Correct sign */
  1514.  
  1515.    return(TRUE);
  1516.  
  1517. }
  1518.  
  1519.  
  1520.  
  1521.  
  1522. /*
  1523.  *    **************************************************
  1524.  *    *                                                *
  1525.  *    *                Resolve Angle                   *
  1526.  *    *                                                *
  1527.  *    *  work[2] = Angle(sin = work[0], cos = work[1]) *
  1528.  *    *                                                *
  1529.  *    *         with angle in first quadrant           *
  1530.  *    *                                                *
  1531.  *    *              (uses 4 temp areas)               *
  1532.  *    *                                                *
  1533.  *    **************************************************
  1534.  */
  1535.  
  1536. static int ResolveAngle(void)
  1537.  
  1538. {
  1539.  
  1540.    long dif = 0, n, order;
  1541.    COMPTYPE *temp;
  1542.    COMPTYPE *sin, *cos, *tsin, *tcos;     /* This group of variables is  */
  1543.    COMPTYPE *x, *y, *sum;                 /*  never used with this group */
  1544.  
  1545.  
  1546.          /* Set up temp registers */
  1547.  
  1548.    if ((temp = GETCOMPTEMP(4)) == NULL) {
  1549.       MemoryError();
  1550.       return(FALSE);
  1551.       }
  1552.    sin  = temp;         /* This group used for scaling */
  1553.    cos  = temp + 1;
  1554.    tsin = temp + 2;
  1555.    tcos = temp + 3;
  1556.  
  1557.    x    = temp;         /* This group used in Euler series */
  1558.    y    = temp + 1;
  1559.    sum  = temp + 2;
  1560.  
  1561.  
  1562.          /* Save off sin and cos */
  1563.  
  1564.    MoveWorkTemp(0, sin);
  1565.    MoveWorkTemp(1, cos);
  1566.  
  1567.  
  1568.    /*
  1569.     *    **************************************************
  1570.     *    *                                                *
  1571.     *    *  Resolve angle A for which we have             *
  1572.     *    *                                                *
  1573.     *    *     sin(A) and cos(A):                         *
  1574.     *    *                                                *
  1575.     *    *                                                *
  1576.     *    *  Scale the arguments sin and cos like this:    *
  1577.     *    *                                                *
  1578.     *    *                                                *
  1579.     *    *  Choose angle A such that                      *
  1580.     *    *                                                *
  1581.     *    *     sin(A) = work[0]  (sin)                    *
  1582.     *    *  and                                           *
  1583.     *    *     cos(A) = work[1]  (cos)                    *
  1584.     *    *                                                *
  1585.     *    *                                                *
  1586.     *    *  Then use the relations:                       *
  1587.     *    *                                                *
  1588.     *    *     sin(A - B) = sin(A)cos(B) - cos(A)sin(B)   *
  1589.     *    *  and                                           *
  1590.     *    *     cos(A - B) = cos(A)cos(B) + sin(A)sin(B)   *
  1591.     *    *                                                *
  1592.     *    *                                                *
  1593.     *    *  Using convenient B (.1, .01 and .001)         *
  1594.     *    *                                                *
  1595.     *    *     and repeatedly adding each B to DIF,       *
  1596.     *    *                                                *
  1597.     *    *     the angle (A - DIF) is reduced until       *
  1598.     *    *                                                *
  1599.     *    *     0 <= sin(A - DIF) <= sin(.001)             *
  1600.     *    *                                                *
  1601.     *    *                                                *
  1602.     *    *  Then the arctangent of sin(A-DIF)/cos(A-DIF)  *
  1603.     *    *                                                *
  1604.     *    *     is calculated and added to DIF to obtain   *
  1605.     *    *                                                *
  1606.     *    *     the resolved angle.                        *
  1607.     *    *                                                *
  1608.     *    *                                                *
  1609.     *    *  dif = DIF * 1000                              *
  1610.     *    *                                                *
  1611.     *    **************************************************
  1612.     */
  1613.  
  1614.  
  1615.          /* Reduce with B = .1 */
  1616.  
  1617.    while (sin->exp > -1) {          /* sin(.1) appx .09 */
  1618.       dif += 100L;                  /* Add .1 * 1000    */
  1619.  
  1620.       MoveTempWork(sin, 0);         /* Calc sin(A)cos(B) */
  1621.       ExtendedRecallCosP1(1);
  1622.       if (! ExtendedMultiply()) {
  1623.          free(temp);
  1624.          return(FALSE);
  1625.          }
  1626.       MoveWorkTemp(2, tsin);        /* Save it */
  1627.  
  1628.       MoveTempWork(cos, 0);         /* Calc cos(A)sin(B) */
  1629.       ExtendedRecallSinP1(1);
  1630.       if (! ExtendedMultiply()) {
  1631.          free(temp);
  1632.          return(FALSE);
  1633.          }
  1634.  
  1635.       MoveWorkWork(2, 0);           /* Calc sin(A)cos(B) - cos(A)sin(B) */
  1636.       MoveTempWork(tsin, 1);
  1637.       ExtendedSubtract();
  1638.       MoveWorkTemp(2, tsin);        /* Save new sin */
  1639.  
  1640.  
  1641.       MoveTempWork(cos, 0);         /* Calc cos(A)cos(B) */
  1642.       ExtendedRecallCosP1(1);
  1643.       if (! ExtendedMultiply()) {
  1644.          free(temp);
  1645.          return(FALSE);
  1646.          }
  1647.       MoveWorkTemp(2, tcos);        /* Save it */
  1648.  
  1649.       MoveTempWork(sin, 0);         /* Calc sin(A)sin(B) */
  1650.       ExtendedRecallSinP1(1);
  1651.       if (! ExtendedMultiply()) {
  1652.          free(temp);
  1653.          return(FALSE);
  1654.          }
  1655.  
  1656.       MoveWorkWork(2, 0);           /* Calc cos(A)cos(B) + sin(A)sin(B) */
  1657.       MoveTempWork(tcos, 1);
  1658.       ExtendedAdd();
  1659.       MoveWorkTemp(2, cos);         /* Save new cos */
  1660.       MoveTempTemp(tsin, sin);      /*  and new sin */
  1661.       }
  1662.  
  1663.  
  1664.          /* Reduce with B = .01 */
  1665.  
  1666.    while (sin->exp > -2) {          /* sin(.1) appx .009 */
  1667.       dif += 10L;                   /* Add .01 * 1000    */
  1668.  
  1669.       MoveTempWork(sin, 0);         /* Calc sin(A)cos(B) */
  1670.       ExtendedRecallCosP01(1);
  1671.       if (! ExtendedMultiply()) {
  1672.          free(temp);
  1673.          return(FALSE);
  1674.          }
  1675.       MoveWorkTemp(2, tsin);        /* Save it */
  1676.  
  1677.       MoveTempWork(cos, 0);         /* Calc cos(A)sin(B) */
  1678.       ExtendedRecallSinP01(1);
  1679.       if (! ExtendedMultiply()) {
  1680.          free(temp);
  1681.          return(FALSE);
  1682.          }
  1683.  
  1684.       MoveWorkWork(2, 0);           /* Calc sin(A)cos(B) - cos(A)sin(B) */
  1685.       MoveTempWork(tsin, 1);
  1686.       ExtendedSubtract();
  1687.       MoveWorkTemp(2, tsin);        /* Save new sin */
  1688.  
  1689.  
  1690.       MoveTempWork(cos, 0);         /* Calc cos(A)cos(B) */
  1691.       ExtendedRecallCosP01(1);
  1692.       if (! ExtendedMultiply()) {
  1693.          free(temp);
  1694.          return(FALSE);
  1695.          }
  1696.       MoveWorkTemp(2, tcos);        /* Save it */
  1697.  
  1698.       MoveTempWork(sin, 0);         /* Calc sin(A)sin(B) */
  1699.       ExtendedRecallSinP01(1);
  1700.       if (! ExtendedMultiply()) {
  1701.          free(temp);
  1702.          return(FALSE);
  1703.          }
  1704.  
  1705.       MoveWorkWork(2, 0);           /* Calc cos(A)cos(B) + sin(A)sin(B) */
  1706.       MoveTempWork(tcos, 1);
  1707.       ExtendedAdd();
  1708.       MoveWorkTemp(2, cos);         /* Save new cos */
  1709.       MoveTempTemp(tsin, sin);      /*  and new sin */
  1710.       }
  1711.  
  1712.  
  1713.          /* Reduce with B = .001 */
  1714.  
  1715.    while (sin->exp > -3) {          /* sin(.1) appx .0009 */
  1716.       dif += 1L;                    /* Add .001 * 1000    */
  1717.  
  1718.       MoveTempWork(sin, 0);         /* Calc sin(A)cos(B) */
  1719.       ExtendedRecallCosP001(1);
  1720.       if (! ExtendedMultiply()) {
  1721.          free(temp);
  1722.          return(FALSE);
  1723.          }
  1724.       MoveWorkTemp(2, tsin);        /* Save it */
  1725.  
  1726.       MoveTempWork(cos, 0);         /* Calc cos(A)sin(B) */
  1727.       ExtendedRecallSinP001(1);
  1728.       if (! ExtendedMultiply()) {
  1729.          free(temp);
  1730.          return(FALSE);
  1731.          }
  1732.  
  1733.       MoveWorkWork(2, 0);           /* Calc sin(A)cos(B) - cos(A)sin(B) */
  1734.       MoveTempWork(tsin, 1);
  1735.       ExtendedSubtract();
  1736.       MoveWorkTemp(2, tsin);        /* Save new sin */
  1737.  
  1738.  
  1739.       MoveTempWork(cos, 0);         /* Calc cos(A)cos(B) */
  1740.       ExtendedRecallCosP001(1);
  1741.       if (! ExtendedMultiply()) {
  1742.          free(temp);
  1743.          return(FALSE);
  1744.          }
  1745.       MoveWorkTemp(2, tcos);        /* Save it */
  1746.  
  1747.       MoveTempWork(sin, 0);         /* Calc sin(A)sin(B) */
  1748.       ExtendedRecallSinP001(1);
  1749.       if (! ExtendedMultiply()) {
  1750.          free(temp);
  1751.          return(FALSE);
  1752.          }
  1753.  
  1754.       MoveWorkWork(2, 0);           /* Calc cos(A)cos(B) + sin(A)sin(B) */
  1755.       MoveTempWork(tcos, 1);
  1756.       ExtendedAdd();
  1757.       MoveWorkTemp(2, cos);         /* Save new cos */
  1758.       MoveTempTemp(tsin, sin);      /*  and new sin */
  1759.       }
  1760.  
  1761.  
  1762.          /* Calculate tan(A) = sin(A)/cos(A) */
  1763.  
  1764.       MoveTempWork(sin, 1);         /* With angle <= .01 */
  1765.       MoveTempWork(cos, 0);         /*  cos(A) > 0       */
  1766.       if (! ExtendedDivide()) {
  1767.          free(temp);
  1768.          return(FALSE);
  1769.          }
  1770.       MoveWorkWork(2, 0);           /* Put in Work[0] */
  1771.  
  1772.  
  1773.    /*
  1774.     *    *******************************************************
  1775.     *    *                                                     *
  1776.     *    *  ARCTAN(X): Using Euler's series:                   *
  1777.     *    *                                                     *
  1778.     *    *                                                     *
  1779.     *    *            y        2    2*4      2*4*6             *
  1780.     *    *  atan(x) = - * (1 + -y + ---y^2 + -----y^3 + ... )  *
  1781.     *    *            x        3    3*5      3*5*7             *
  1782.     *    *                                                     *
  1783.     *    *                                                     *
  1784.     *    *     where = y = x^2 / (1 + x^2)                     *
  1785.     *    *                                                     *
  1786.     *    *                                                     *
  1787.     *    *     *x    = Argument                                *
  1788.     *    *                                                     *
  1789.     *    *     *y    = x^2 / (1 + x^2)                         *
  1790.     *    *                                                     *
  1791.     *    *     n     = Term counter                            *
  1792.     *    *                                                     *
  1793.     *    *******************************************************
  1794.     */
  1795.  
  1796.  
  1797.    MoveWorkTemp(0, x);        /* Save x */
  1798.  
  1799.  
  1800.          /* Calculate y = x^2 / (1 + x^2) */
  1801.  
  1802.    MoveWorkWork(0, 1);        /* Calculate x^2 */
  1803.    if (! ExtendedMultiply()) {
  1804.       free(temp);
  1805.       return(FALSE);
  1806.       }
  1807.  
  1808.    MoveWorkWork(2, 1);        /* Calculate 1 + x^2 */
  1809.    SetWorkInteger(0, 1L);
  1810.    ExtendedAdd();             /* Doesn't change x^2 in work[1] */
  1811.  
  1812.    MoveWorkWork(2, 0);
  1813.    if (! ExtendedDivide()) {
  1814.       free(temp);
  1815.       return(FALSE);
  1816.       }
  1817.    MoveWorkTemp(2, y);        /* Store y */
  1818.  
  1819.  
  1820.    SetTempInteger(sum, 1L);   /* Sum */
  1821.  
  1822.    n = 0;                     /* N */
  1823.  
  1824.    SetWorkInteger(0, 1L);     /* Term in work[0] */
  1825.  
  1826.  
  1827.  
  1828.    while (TRUE) {             /* Until term < compprec */
  1829.  
  1830.       n += 2;                       /* n = term number */
  1831.  
  1832.          /* term *= n  (term still in work[0]) */
  1833.  
  1834.       SetWorkInteger(1, n);         /* Put n in work[1] */
  1835.       if (! ExtendedMultiply()) {
  1836.          free(temp);
  1837.          return(FALSE);
  1838.          }
  1839.  
  1840.  
  1841.          /* term *= y */
  1842.  
  1843.       MoveTempWork(y, 0);           /* Get y */
  1844.       MoveWorkWork(2, 1);           /* Get term */
  1845.       if (! ExtendedMultiply()) {
  1846.          free(temp);
  1847.          return(FALSE);
  1848.          }
  1849.  
  1850.  
  1851.          /* term /= (n + 1) */
  1852.  
  1853.       SetWorkInteger(0, (n + 1L));  /* Get n + 1 */
  1854.       MoveWorkWork(2, 1);           /* Get term */
  1855.       if (! ExtendedDivide()) {
  1856.          free(temp);
  1857.          return(FALSE);
  1858.          }
  1859.  
  1860.  
  1861.          /* sum += term */
  1862.  
  1863.       MoveWorkWork(2, 0);           /* Get term */
  1864.       MoveTempWork(sum, 1);         /* Get sum */
  1865.       if (! ExtendedAdd()) {
  1866.          free(temp);
  1867.          return(FALSE);
  1868.          }
  1869.  
  1870.  
  1871.          /* Test for completion */
  1872.  
  1873.       if ( (labs(work[2].exp - work[0].exp) >= compprec)
  1874.             ||
  1875.            (! work[0].digits) ) {
  1876.          break;
  1877.          }
  1878.  
  1879.  
  1880.          /* Save new sum */
  1881.  
  1882.       MoveWorkTemp(2, sum);
  1883.  
  1884.       }  /* while */
  1885.  
  1886.  
  1887.          /* Compute angle = sum * y / x */
  1888.  
  1889.    MoveWorkWork(2, 0);        /* Get sum */
  1890.    MoveTempWork(y, 1);
  1891.    if (! ExtendedMultiply()) {
  1892.       free(temp);
  1893.       return(FALSE);
  1894.       }
  1895.    MoveWorkWork(2, 1);
  1896.    MoveTempWork(x, 0);
  1897.    if (! ExtendedDivide()) {
  1898.       free(temp);
  1899.       return(FALSE);
  1900.       }
  1901.  
  1902.    free(temp);
  1903.  
  1904.          /* Resolve dif scaling */
  1905.  
  1906.    MoveWorkWork(2, 0);        /* Get angle */
  1907.  
  1908.    work[1].exp  = 1L;         /* Put dif in work[1] */
  1909.    work[1].sign = '+';
  1910.    work[1].digits = 0;
  1911.    order = 1000L;
  1912.    while (order > dif) {
  1913.       order /= 10;
  1914.       work[1].exp--;
  1915.       }
  1916.    while (order) {
  1917.       work[1].man[work[1].digits] = (WORKDIGIT) (dif / order);
  1918.       dif %= order;
  1919.       order /= 10L;
  1920.       work[1].digits++;
  1921.       }
  1922.    memset(&work[1].man[work[1].digits], 0, 
  1923.                ((workprec - work[1].digits) * sizeof(WORKDIGIT)));
  1924.  
  1925.    return(ExtendedAdd());
  1926.  
  1927. }
  1928.  
  1929.  
  1930.  
  1931.  
  1932. /*
  1933.  *    **************************************************
  1934.  *    *                                                *
  1935.  *    *           Extended Common Logarithm            *
  1936.  *    *                                                *
  1937.  *    *             work[2] = Log(work[0])             *
  1938.  *    *                                                *
  1939.  *    **************************************************
  1940.  */
  1941.  
  1942. extern int ExtendedLog(void)
  1943.  
  1944. {
  1945.  
  1946.  
  1947.    /*
  1948.     *    **************************************************
  1949.     *    *                                                *
  1950.     *    *  LOG  (X): Using this relation:                *
  1951.     *    *     10                                         *
  1952.     *    *                                                *
  1953.     *    *              lnX                               *
  1954.     *    *     log  X = ---                               *
  1955.     *    *        10    ln10                              *
  1956.     *    *                                                *
  1957.     *    **************************************************
  1958.     */
  1959.  
  1960.    if (! ExtendedLn()) {         /* Get lnX */
  1961.       return(FALSE);
  1962.       }
  1963.    MoveWorkWork(2, 1);
  1964.  
  1965.    ExtendedRecallLn10(0);        /* Get ln10 */
  1966.  
  1967.    if (! ExtendedDivide()) {
  1968.       return(FALSE);
  1969.       }
  1970.  
  1971.    return(TRUE);
  1972.  
  1973.  
  1974. }
  1975.  
  1976.  
  1977.  
  1978.  
  1979. /*
  1980.  *    **************************************************
  1981.  *    *                                                *
  1982.  *    *          Extended Common Antilogarithm         *
  1983.  *    *                                                *
  1984.  *    *              work[2] = 10^(work[0])            *
  1985.  *    *                                                *
  1986.  *    **************************************************
  1987.  */
  1988.  
  1989. extern int ExtendedExp10(void)
  1990.  
  1991. {
  1992.  
  1993.    int i;
  1994.    long pow = 0L;
  1995.  
  1996.  
  1997.          /* Special Cases */
  1998.  
  1999.    if (! work[0].digits) {                /* Arg zero: return one */
  2000.       SetWorkInteger(2, 1L);
  2001.       return(TRUE);
  2002.       }
  2003.  
  2004.    if (! work[0].exp > MAXEXDIGITS) {
  2005.       if (! work[0].sign == '-') {        /* Underflow: return zero */
  2006.          ClearWork(2);
  2007.          return(TRUE);
  2008.          }
  2009.       else {
  2010.          Overflow();                      /* Overflow */
  2011.          return(FALSE);
  2012.          }
  2013.       }
  2014.  
  2015.  
  2016.  
  2017.          /* If X integer: X -> exp */
  2018.  
  2019.    if (work[0].exp >= work[0].digits) {
  2020.  
  2021.       pow = 0;                            /* X -> pow */
  2022.       for (i = 0; i < (int)work[0].exp; i++) {
  2023.          pow = pow * 10 + work[0].man[i];
  2024.          }
  2025.       if (work[0].sign == '-') {
  2026.          pow = - pow;
  2027.          }
  2028.       pow++;
  2029.  
  2030.       work[2].exp = pow;
  2031.       work[2].sign = '+';
  2032.       work[2].digits = 1;
  2033.       work[2].man[0] = 1;
  2034.  
  2035.       return(TRUE);
  2036.  
  2037.       }  /* if */
  2038.  
  2039.  
  2040.    /*
  2041.     *    **************************************************
  2042.     *    *                                                *
  2043.     *    *  ALOG: 10^X using this:                        *
  2044.     *    *                                                *
  2045.     *    *     10^X = e^(X*ln10)                          *
  2046.     *    *                                                *
  2047.     *    **************************************************
  2048.     */
  2049.  
  2050.  
  2051.          /* Compute X * ln10 */
  2052.  
  2053.    ExtendedRecallLn10(1);           /* Get ln10 */
  2054.    if (! ExtendedMultiply()) {
  2055.       return(FALSE);
  2056.       }
  2057.  
  2058.          /* Compute e^(X * ln10) */
  2059.  
  2060.    MoveWorkWork(2, 0);              /* Get X * ln10 */
  2061.    if (! ExtendedExpE()) {
  2062.       return(FALSE);
  2063.       }
  2064.  
  2065.    return(TRUE);
  2066.  
  2067. }
  2068.  
  2069.  
  2070.  
  2071.  
  2072. /*
  2073.  *    **************************************************
  2074.  *    *                                                *
  2075.  *    *          Extended Natural Logarithm            *
  2076.  *    *                                                *
  2077.  *    *             work[2] = ln(work[0])              *
  2078.  *    *                                                *
  2079.  *    *              (Uses 3 temp areas)               *
  2080.  *    *                                                *
  2081.  *    **************************************************
  2082.  */
  2083.  
  2084. extern int ExtendedLn(void)
  2085.  
  2086. {
  2087.  
  2088.    int mulp9 = 0, mulp99 = 0, mulp999 = 0, mulp9999 = 0, mulp99999 = 0;
  2089.    int dig, bor;
  2090.    long den, pow10 = 0;
  2091.    WORKDIGIT *pl, *pr;
  2092.    COMPTYPE *temp;
  2093.    COMPTYPE *ln, *ysq, *num;
  2094.  
  2095.  
  2096.          /* Special Cases */
  2097.  
  2098.    if (! work[0].digits) {          /* Zero: error */
  2099.       ZeroArgument();
  2100.       return(FALSE);
  2101.       }
  2102.  
  2103.    if (work[0].sign == '-') {       /* Negative: error */
  2104.       NegativeArgument();
  2105.       return(FALSE);
  2106.       }
  2107.  
  2108.    if ( (work[0].digits == 1)       /* One: return zero */
  2109.          &&
  2110.         (work[0].exp == 1L)
  2111.          &&
  2112.         (work[0].man[0] == 1) ) {
  2113.       ClearWork(2);
  2114.       return(TRUE);
  2115.       }
  2116.  
  2117.  
  2118.    /*
  2119.     *    **************************************************
  2120.     *    *                                                *
  2121.     *    *  Scale: 1 <= X < 10  like this:                *
  2122.     *    *                                                *
  2123.     *    *     M = Mantissa, EX = Exponent                *
  2124.     *    *                                                *
  2125.     *    *     ln(X) = ln(M * 10^EX)                      *
  2126.     *    *                                                *
  2127.     *    *           = ln(M) + ln(10^EX)                  *
  2128.     *    *                                                *
  2129.     *    *           = ln(M) + EX * ln(10)                *
  2130.     *    *                                                *
  2131.     *    *                                                *
  2132.     *    *     pow10 = EX - 1  (so 1 <= M < 10)           *
  2133.     *    *                                                *
  2134.     *    **************************************************
  2135.     */
  2136.  
  2137.    if (work[0].exp != 1L) {
  2138.       pow10 = work[0].exp - 1;
  2139.       work[0].exp = 1;
  2140.       }
  2141.  
  2142.  
  2143.    /*
  2144.     *    **************************************************
  2145.     *    *                                                *
  2146.     *    *  Scale: 1 <= X <= 1.001 like this:             *
  2147.     *    *                                                *
  2148.     *    *     ln(X) = ln(X*F/F)                          *
  2149.     *    *                                                *
  2150.     *    *           = ln(X*F) - ln(F)                    *
  2151.     *    *                                                *
  2152.     *    *                                                *
  2153.     *    *     mulpF = nbr times X*F for X <= 1.00001     *
  2154.     *    *                                                *
  2155.     *    **************************************************
  2156.     */
  2157.  
  2158.  
  2159.  
  2160.    /*
  2161.     *    **************************************************
  2162.     *    *                                                *
  2163.     *    *  Calculate X * .9, .99, .999, etc. like this:  *
  2164.     *    *                                                *
  2165.     *    *     X * .9 = X - X/10                          *
  2166.     *    *                                                *
  2167.     *    *     X * .99 = X - X/100                        *
  2168.     *    *                                                *
  2169.     *    *     X * .999 = X - X/1000                      *
  2170.     *    *                                                *
  2171.     *    *     X * .9999 = X - X/10000                    *
  2172.     *    *                                                *
  2173.     *    *     X * .99999 = X - X/100000                  *
  2174.     *    *                                                *
  2175.     *    *                                                *
  2176.     *    *  The subtractions are done in place            *
  2177.     *    *                                                *
  2178.     *    **************************************************
  2179.     */
  2180.  
  2181.             /* F = .9 */
  2182.  
  2183.    while ( (work[0].man[0] > 1)        /* While X > 1.1 */
  2184.             ||
  2185.            (work[0].man[1] > 1) ) {
  2186.       pl = &work[0].man[0];
  2187.       pr = &work[0].man[work[0].digits];
  2188.       if (work[0].digits < compprec)
  2189.          work[0].digits += 1;
  2190.       dig = 0;
  2191.       bor = 0;
  2192.       while (--pr >= pl) {
  2193.          dig = *(pr + 1) - *pr - bor;
  2194.          if (dig < 0) {
  2195.             bor = 1;
  2196.             dig += 10;
  2197.             }
  2198.          else {
  2199.             bor = 0;
  2200.             }
  2201.          *(pr + 1) = dig;
  2202.          }
  2203.       work[0].man[0] -= bor;
  2204.       mulp9++;
  2205.       }
  2206.  
  2207.  
  2208.             /* F = .99 */
  2209.  
  2210.    while ( (work[0].man[1] > 0)        /* While X > 1.01 */
  2211.             ||
  2212.            (work[0].man[2] > 1) ) {
  2213.       pl = &work[0].man[0];
  2214.       pr = &work[0].man[work[0].digits];
  2215.       if (work[0].digits < compprec)
  2216.          work[0].digits += 2;
  2217.       dig = 0;
  2218.       bor = 0;
  2219.       while (--pr >= pl) {
  2220.          dig = *(pr + 2) - *pr - bor;
  2221.          if (dig < 0) {
  2222.             bor = 1;
  2223.             dig += 10;
  2224.             }
  2225.          else {
  2226.             bor = 0;
  2227.             }
  2228.          *(pr + 2) = dig;
  2229.          }
  2230.       work[0].man[1] -= bor;
  2231.       mulp99++;
  2232.       }
  2233.  
  2234.  
  2235.             /* F = .999 */
  2236.  
  2237.    while ( (work[0].man[2] > 0)        /* While X > 1.001 */
  2238.             ||
  2239.            (work[0].man[3] > 1) ) {
  2240.       pl = &work[0].man[0];
  2241.       pr = &work[0].man[work[0].digits];
  2242.       if (work[0].digits < compprec)
  2243.          work[0].digits += 3;
  2244.       dig = 0;
  2245.       bor = 0;
  2246.       while (--pr >= pl) {
  2247.          dig = *(pr + 3) - *pr - bor;
  2248.          if (dig < 0) {
  2249.             bor = 1;
  2250.             dig += 10;
  2251.             }
  2252.          else {
  2253.             bor = 0;
  2254.             }
  2255.          *(pr + 3) = dig;
  2256.          }
  2257.       work[0].man[2] -= bor;
  2258.       mulp999++;
  2259.       }
  2260.  
  2261.  
  2262.             /* F = .9999 */
  2263.  
  2264.    while ( (work[0].man[3] > 0)        /* While X > 1.0001 */
  2265.             ||
  2266.            (work[0].man[4] > 1) ) {
  2267.       pl = &work[0].man[0];
  2268.       pr = &work[0].man[work[0].digits];
  2269.       if (work[0].digits < compprec)
  2270.          work[0].digits += 4;
  2271.       dig = 0;
  2272.       bor = 0;
  2273.       while (--pr >= pl) {
  2274.          dig = *(pr + 4) - *pr - bor;
  2275.          if (dig < 0) {
  2276.             bor = 1;
  2277.             dig += 10;
  2278.             }
  2279.          else {
  2280.             bor = 0;
  2281.             }
  2282.          *(pr + 4) = dig;
  2283.          }
  2284.       work[0].man[3] -= bor;
  2285.       mulp9999++;
  2286.       }
  2287.  
  2288.  
  2289.             /* F = .99999 */
  2290.  
  2291.    while ( (work[0].man[4] > 0)        /* While X > 1.00001 */
  2292.             ||
  2293.            (work[0].man[5] > 1) ) {
  2294.       pl = &work[0].man[0];
  2295.       pr = &work[0].man[work[0].digits];
  2296.       if (work[0].digits < compprec)
  2297.          work[0].digits += 5;
  2298.       dig = 0;
  2299.       bor = 0;
  2300.       while (--pr >= pl) {
  2301.          dig = *(pr + 5) - *pr - bor;
  2302.          if (dig < 0) {
  2303.             bor = 1;
  2304.             dig += 10;
  2305.             }
  2306.          else {
  2307.             bor = 0;
  2308.             }
  2309.          *(pr + 5) = dig;
  2310.          }
  2311.       work[0].man[4] -= bor;
  2312.       mulp99999++;
  2313.       }
  2314.  
  2315.  
  2316.          /* Set up temp registers */
  2317.  
  2318.    if ((temp = GETCOMPTEMP(3)) == NULL) {
  2319.       MemoryError();
  2320.       return(FALSE);
  2321.       }
  2322.    ln  = temp;
  2323.    ysq = temp + 1;
  2324.    num = temp + 2;
  2325.  
  2326.  
  2327.  
  2328.    /*
  2329.     *    **************************************************
  2330.     *    *                                                *
  2331.     *    *  LOG (X): Using the power series:              *
  2332.     *    *     e                                          *
  2333.     *    *                                                *
  2334.     *    *                X-1          1+Y                *
  2335.     *    *  We choose Y = --- then X = --- so ln X =      *
  2336.     *    *                X+1          1-Y                *
  2337.     *    *                                                *
  2338.     *    *                                                *
  2339.     *    *     1+Y             Y^3   Y^5   Y^7            *
  2340.     *    *  ln --- = 2 * ( Y + --- + --- + --- + ... )    *
  2341.     *    *     1-Y              3     5     7             *
  2342.     *    *                                                *
  2343.     *    *                                                *
  2344.     *    *  *ln  = Sum of series                          *
  2345.     *    *                                                *
  2346.     *    *  *ysq = y^2                                    *
  2347.     *    *                                                *
  2348.     *    *  *num = term numerator                         *
  2349.     *    *                                                *
  2350.     *    **************************************************
  2351.     */
  2352.  
  2353.  
  2354.          /* Compute Y = (X - 1) / (X + 1) */
  2355.  
  2356.    MoveWorkWork(0, 1);              /* X -> work[1] */
  2357.    SetWorkInteger(0, 1L);           /* 1 -> work[0] */
  2358.    if (! ExtendedSubtract()) {      /* Compute (X - 1) */
  2359.       free(temp);
  2360.       return(FALSE);
  2361.       }
  2362.    MoveWorkTemp(2, ysq);            /* Save (X - 1) */
  2363.  
  2364.  
  2365.    if (! ExtendedAdd()) {           /* Compute (X + 1) */
  2366.       free(temp);
  2367.       return(FALSE);
  2368.       }
  2369.    MoveWorkWork(2, 0);              /* (X + 1) -> work[0] */
  2370.    MoveTempWork(ysq, 1);            /* (X - 1) -> work[1] */
  2371.    if (! ExtendedDivide()) {        /* Compute (X - 1) / (X + 1) */
  2372.       free(temp);
  2373.       return(FALSE);
  2374.       }
  2375.  
  2376.    MoveWorkTemp(2, num);            /* Y -> num */
  2377.    MoveWorkTemp(2, ln);             /* Y -> ln */
  2378.  
  2379.  
  2380.          /* Compute Y^2 */
  2381.  
  2382.    MoveWorkWork(2, 0);              /* Get Y */
  2383.    MoveWorkWork(2, 1);
  2384.    if (! ExtendedMultiply()) {
  2385.       free(temp);                   /* Compute Y^2 */
  2386.       return(FALSE);
  2387.       }
  2388.    MoveWorkTemp(2, ysq);            /* Save Y^2 */
  2389.  
  2390.  
  2391.    den = 1L;                        /* den = 1 */
  2392.  
  2393.  
  2394.    /*
  2395.     *    **************************************************
  2396.     *    *                                                *
  2397.     *    *   LN (X): Where X = (1 + Y) / (1 - Y):         *
  2398.     *    *     e                                          *
  2399.     *    *                                                *
  2400.     *    *      1+Y             Y^3   Y^5   Y^7           *
  2401.     *    *   ln --- = 2 * ( Y + --- + --- + --- + ... )   *
  2402.     *    *      1-Y              3     5     7            *
  2403.     *    *                                                *
  2404.     *    **************************************************
  2405.     */
  2406.  
  2407.    while (TRUE) {             /* Until term < compprec */
  2408.  
  2409.  
  2410.             /* num *= Y^2 */
  2411.  
  2412.       MoveTempWork(num, 0);
  2413.       MoveTempWork(ysq, 1);
  2414.       if (! ExtendedMultiply()) {
  2415.          free(temp);
  2416.          return(FALSE);
  2417.          }
  2418.       MoveWorkTemp(2, num);            /* Save new num */
  2419.  
  2420.  
  2421.             /* den += 2 */
  2422.  
  2423.       den += 2;
  2424.  
  2425.  
  2426.             /* num / den */
  2427.  
  2428.       SetWorkInteger(0, den);          /* Get den */
  2429.       MoveWorkWork(2, 1);              /* Get num */
  2430.       if (! ExtendedDivide()) {
  2431.          free(temp);
  2432.          return(FALSE);
  2433.          }
  2434.  
  2435.             /* ln + num / den */
  2436.  
  2437.       MoveWorkWork(2, 0);
  2438.       MoveTempWork(ln, 1);
  2439.       if (! ExtendedAdd()) {
  2440.          free(temp);
  2441.          return(FALSE);
  2442.          }
  2443.  
  2444.             /* Test for completion */
  2445.  
  2446.       if ( (labs(work[2].exp - work[0].exp) >= compprec)
  2447.             ||
  2448.            (! work[2].digits) ) {
  2449.          break;
  2450.          }
  2451.  
  2452.       MoveWorkTemp(2, ln);             /* Save new ln */
  2453.  
  2454.       }  /* while */
  2455.  
  2456.  
  2457.             /* Calculate 2 * ( ... ) */
  2458.  
  2459.    MoveWorkWork(2, 0);
  2460.    MoveWorkWork(2, 0);
  2461.    if (! ExtendedAdd()) {
  2462.       free(temp);
  2463.       return(FALSE);
  2464.       }
  2465.  
  2466.  
  2467.             /* Resolve .99999 scaling */
  2468.  
  2469.    if (mulp99999) {
  2470.       if (mulp99999 == 1) {
  2471.          ExtendedRecallLnP99999(0);       /* ln(.99999) -> work[0] */
  2472.          MoveWorkWork(2, 1);
  2473.          }
  2474.       else {
  2475.          MoveWorkTemp(2, ln);             /* Save ln */
  2476.          SetWorkInteger(0, (long)mulp99999);
  2477.          ExtendedRecallLnP99999(1);       /* ln(.99999) -> work[1] */
  2478.          if (! ExtendedMultiply()) {
  2479.             free(temp);
  2480.             return(FALSE);
  2481.             }
  2482.          MoveWorkWork(2, 0);
  2483.          MoveTempWork(ln, 1);
  2484.          }
  2485.       if (! ExtendedSubtract()) {
  2486.          free(temp);
  2487.          return(FALSE);
  2488.          }
  2489.       }
  2490.  
  2491.  
  2492.             /* Resolve .9999 scaling */
  2493.  
  2494.    if (mulp9999) {
  2495.       if (mulp9999 == 1) {
  2496.          ExtendedRecallLnP9999(0);        /* ln(.9999) -> work[0] */
  2497.          MoveWorkWork(2, 1);
  2498.          }
  2499.       else {
  2500.          MoveWorkTemp(2, ln);             /* Save ln */
  2501.          SetWorkInteger(0, (long)mulp9999);
  2502.          ExtendedRecallLnP9999(1);        /* ln(.9999) -> work[1] */
  2503.          if (! ExtendedMultiply()) {
  2504.             free(temp);
  2505.             return(FALSE);
  2506.             }
  2507.          MoveWorkWork(2, 0);
  2508.          MoveTempWork(ln, 1);
  2509.          }
  2510.       if (! ExtendedSubtract()) {
  2511.          free(temp);
  2512.          return(FALSE);
  2513.          }
  2514.       }
  2515.  
  2516.  
  2517.             /* Resolve .999 scaling */
  2518.  
  2519.    if (mulp999) {
  2520.       if (mulp999 == 1) {
  2521.          ExtendedRecallLnP999(0);         /* ln(.999) -> work[0] */
  2522.          MoveWorkWork(2, 1);
  2523.          }
  2524.       else {
  2525.          MoveWorkTemp(2, ln);             /* Save ln */
  2526.          SetWorkInteger(0, (long)mulp999);
  2527.          ExtendedRecallLnP999(1);         /* ln(.999) -> work[1] */
  2528.          if (! ExtendedMultiply()) {
  2529.             free(temp);
  2530.             return(FALSE);
  2531.             }
  2532.          MoveWorkWork(2, 0);
  2533.          MoveTempWork(ln, 1);
  2534.          }
  2535.       if (! ExtendedSubtract()) {
  2536.          free(temp);
  2537.          return(FALSE);
  2538.          }
  2539.       }
  2540.  
  2541.  
  2542.             /* Resolve .99 scaling */
  2543.  
  2544.    if (mulp99) {
  2545.       if (mulp99 == 1) {
  2546.          ExtendedRecallLnP99(0);          /* ln(.99) -> work[0] */
  2547.          MoveWorkWork(2, 1);
  2548.          }
  2549.       else {
  2550.          MoveWorkTemp(2, ln);             /* Save ln */
  2551.          SetWorkInteger(0, (long)mulp99);
  2552.          ExtendedRecallLnP99(1);          /* ln(.99) -> work[1] */
  2553.          if (! ExtendedMultiply()) {
  2554.             free(temp);
  2555.             return(FALSE);
  2556.             }
  2557.          MoveWorkWork(2, 0);
  2558.          MoveTempWork(ln, 1);
  2559.          }
  2560.       if (! ExtendedSubtract()) {
  2561.          free(temp);
  2562.          return(FALSE);
  2563.          }
  2564.       }
  2565.  
  2566.  
  2567.             /* Resolve .9 scaling */
  2568.  
  2569.    if (mulp9) {
  2570.       if (mulp9 == 1) {
  2571.          ExtendedRecallLnP9(0);           /* ln(.9) -> work[0] */
  2572.          MoveWorkWork(2, 1);
  2573.          }
  2574.       else {
  2575.          MoveWorkTemp(2, ln);             /* Save ln */
  2576.          SetWorkInteger(0, (long)mulp9);
  2577.          ExtendedRecallLnP9(1);           /* ln(.9) -> work[1] */
  2578.          if (! ExtendedMultiply()) {
  2579.             free(temp);
  2580.             return(FALSE);
  2581.             }
  2582.          MoveWorkWork(2, 0);
  2583.          MoveTempWork(ln, 1);
  2584.          }
  2585.       if (! ExtendedSubtract()) {
  2586.          free(temp);
  2587.          return(FALSE);
  2588.          }
  2589.       }
  2590.  
  2591.  
  2592.             /* Resolve pow10 scaling */
  2593.  
  2594.    if (pow10) {
  2595.       MoveWorkTemp(2, ln);             /* Save ln */
  2596.       ExtendedRecallLn10(0);           /* ln(10) -> work[0] */
  2597.       SetWorkInteger(1, pow10);        /* pow10 -> work[1] */
  2598.       if (! ExtendedMultiply()) {
  2599.          free(temp);
  2600.          return(FALSE);
  2601.          }
  2602.       MoveTempWork(ln, 0);             /* ln -> work[0] */
  2603.       MoveWorkWork(2, 1);              /* (ln(10) * pow10) -> work[1] */
  2604.       if (! ExtendedAdd()) {
  2605.          free(temp);
  2606.          return(FALSE);
  2607.          }
  2608.       }
  2609.  
  2610.    free(temp);
  2611.  
  2612.    return(TRUE);
  2613.  
  2614.  
  2615. }
  2616.  
  2617.  
  2618.  
  2619.  
  2620. /*
  2621.  *    **************************************************
  2622.  *    *                                                *
  2623.  *    *         Extended Natural Antilogarithm         *
  2624.  *    *                                                *
  2625.  *    *              work[2] = e^(work[0])             *
  2626.  *    *                                                *
  2627.  *    *               (Uses 2 temp areas)              *
  2628.  *    *                                                *
  2629.  *    **************************************************
  2630.  */
  2631.  
  2632. extern int ExtendedExpE(void)
  2633.  
  2634. {
  2635.  
  2636.    int i, j;
  2637.    long n, pow2 = 1, powe = 0;
  2638.    COMPTYPE *temp;
  2639.    COMPTYPE *arg, *exp;
  2640.    BOOLEAN recip = FALSE;
  2641.  
  2642.  
  2643.          /* Special Cases */
  2644.  
  2645.    if (! work[0].digits) {                   /* X zero: return 1 */
  2646.       SetWorkInteger(2, 1L);
  2647.       return(TRUE);
  2648.       }
  2649.  
  2650.    if (work[0].exp > MAXEXDIGITS) {          /* Overflow */
  2651.       Overflow();
  2652.       return(FALSE);
  2653.       }
  2654.  
  2655.  
  2656.    if ( (work[0].exp >= work[0].digits)      /* If X integer */
  2657.          &&
  2658.         (work[0].exp <= MAXEXDIGITS) ) {     /*  and 'small' */
  2659.       ExtendedRecallE(1);
  2660.       return(ExtendedPower());               /*  compute directly */
  2661.       }
  2662.  
  2663.  
  2664.          /* Process neg arg as pos & take recip at exit */
  2665.  
  2666.    if (work[0].sign == '-') {
  2667.       recip = TRUE;
  2668.       work[0].sign = '+';
  2669.       }
  2670.  
  2671.  
  2672.    /*
  2673.     *    **************************************************
  2674.     *    *                                                *
  2675.     *    *  Scale X like this:                            *
  2676.     *    *                                                *
  2677.     *    *     e^X = e^(iii.ddd)                          *
  2678.     *    *                                                *
  2679.     *    *         = e^(iii + .ddd)                       *
  2680.     *    *                                                *
  2681.     *    *         = e^(iii) * e^(.ddd)                   *
  2682.     *    *                                                *
  2683.     *    *                                                *
  2684.     *    *     powe = iii                                 *
  2685.     *    *                                                *
  2686.     *    *     *arg = .ddd                                *
  2687.     *    *                                                *
  2688.     *    **************************************************
  2689.     */
  2690.  
  2691.    if (work[0].exp > 0) {                 /* Test for iii part */
  2692.  
  2693.       j = (int)work[0].exp;               /* Extract iii -> powe */
  2694.       for (i = 0; i < j; i++) {
  2695.          powe = powe * 10 + work[0].man[i];
  2696.          work[0].man[i] = 0;
  2697.          }
  2698.  
  2699.       Normalize(0);                       /* Normalize .ddd */
  2700.       }
  2701.  
  2702.    /*
  2703.     *    **************************************************
  2704.     *    *                                                *
  2705.     *    *  Scale X like this:                            *
  2706.     *    *                                                *
  2707.     *    *     e^X = e^(X/n * n)                          *
  2708.     *    *                                                *
  2709.     *    *         = (e^(X/n))^n                          *
  2710.     *    *                                                *
  2711.     *    *                                                *
  2712.     *    *     pow2 = n (powers of 2)                     *
  2713.     *    *                                                *
  2714.     *    **************************************************
  2715.     */
  2716.  
  2717.    while (work[0].exp > -7L) {            /* Scale to X < .0000001 */
  2718.  
  2719.       ClearWork(2);                       /* Multiply X by .5 */
  2720.       work[2].exp    = work[0].exp;
  2721.       work[2].sign   = work[0].sign;
  2722.       work[2].digits = work[0].digits + 1;
  2723.       for (i = work[0].digits - 1; i >= 0; i--) {
  2724.          if ((work[2].man[i+1] += work[0].man[i] * 5) > 9) {
  2725.             work[2].man[i] = work[2].man[i+1] / 10;
  2726.             work[2].man[i+1] %= 10;
  2727.             }
  2728.          }
  2729.       if ( (! work[2].man[0])
  2730.             ||
  2731.            (work[2].digits > compprec) ) {
  2732.          Normalize(2);
  2733.          }
  2734.  
  2735.       MoveWorkWork(2, 0);                 /* Save new arg */
  2736.       pow2 *= 2;
  2737.       }
  2738.  
  2739.  
  2740.          /* Set up temp registers */
  2741.  
  2742.    if ((temp = GETCOMPTEMP(2)) == NULL) {
  2743.       MemoryError();
  2744.       return(FALSE);
  2745.       }
  2746.    arg  = temp;
  2747.    exp  = temp + 1;
  2748.  
  2749.    MoveWorkTemp(0, arg);                  /* Save arg */
  2750.  
  2751.  
  2752.          /* Calculate exp to first 2 terms */
  2753.  
  2754.    SetWorkInteger(1, 1L);                 /* Compute sum of terms 0 & 1: */
  2755.    if (! ExtendedAdd()) {                 /*  1 + X   (work[0] == X) */
  2756.       free(temp);
  2757.       return(FALSE);
  2758.       }
  2759.    MoveWorkTemp(2, exp);                  /*  and -> exp */
  2760.  
  2761.    n = 1L;                                /* Terms 0 & 1 done */
  2762.  
  2763.  
  2764.    /*
  2765.     *    **************************************************
  2766.     *    *                                                *
  2767.     *    *  EXPONENT (e^X): Using the power series:       *
  2768.     *    *                                                *
  2769.     *    *                                                *
  2770.     *    *                   X^2   X^3   X^4              *
  2771.     *    *     e^X = 1 + X + --- + --- + --- + ...        *
  2772.     *    *                    2!    3!    4!              *
  2773.     *    *                                                *
  2774.     *    *                                                *
  2775.     *    *     *arg = X                                   *
  2776.     *    *                                                *
  2777.     *    *     *exp = Sum of series                       *
  2778.     *    *                                                *
  2779.     *    *     n    = Term number                         *
  2780.     *    *                                                *
  2781.     *    **************************************************
  2782.     */
  2783.  
  2784.  
  2785.    while (TRUE) {             /* Until term < compprec */
  2786.  
  2787.       n++;                          /* n = term number */
  2788.  
  2789.          /* term *= X */
  2790.  
  2791.          /* term still in work[0] */
  2792.  
  2793.       MoveTempWork(arg, 1);         /* Get arg */
  2794.       if (! ExtendedMultiply()) {
  2795.          free(temp);
  2796.          return(FALSE);
  2797.          }
  2798.  
  2799.  
  2800.          /* term /= n */
  2801.  
  2802.       SetWorkInteger(0, n);         /* Get n */
  2803.       MoveWorkWork(2, 1);           /* Get term */
  2804.       if (! ExtendedDivide()) {
  2805.          free(temp);
  2806.          return(FALSE);
  2807.          }
  2808.  
  2809.  
  2810.          /* Add to exp */
  2811.  
  2812.       MoveWorkWork(2, 0);           /* new term */
  2813.       MoveTempWork(exp, 1);         /* Get exp */
  2814.       if (! ExtendedAdd()) {
  2815.          free(temp);
  2816.          return(FALSE);
  2817.          }
  2818.  
  2819.  
  2820.          /* Test for completion */
  2821.  
  2822.       if ( (labs(work[2].exp - work[0].exp) >= compprec)
  2823.             ||
  2824.            (! work[0].digits) ) {
  2825.          break;
  2826.          }
  2827.  
  2828.  
  2829.          /* Save new exp */
  2830.  
  2831.       MoveWorkTemp(2, exp);
  2832.  
  2833.       }  /* while */
  2834.  
  2835.  
  2836.  
  2837.          /* Resolve pow2 scaling */
  2838.  
  2839.    if (pow2 > 1) {
  2840.       SetWorkInteger(0, pow2);
  2841.       MoveWorkWork(2, 1);
  2842.       if (! ExtendedPower()) {
  2843.          free(temp);
  2844.          return(FALSE);
  2845.          }
  2846.       }
  2847.  
  2848.  
  2849.          /* Resolve powe scaling */
  2850.  
  2851.    if (powe) {
  2852.       MoveWorkTemp(2, arg);            /* Save e^(.ddd) part */
  2853.  
  2854.       SetWorkInteger(0, powe);         /* Calculate e^(iii) */
  2855.       ExtendedRecallE(1);
  2856.       if (! ExtendedPower()) {
  2857.          free(temp);
  2858.          return(FALSE);
  2859.          }
  2860.  
  2861.       MoveWorkWork(2, 0);              /* Calculate e^(iii) * e^(.ddd) */
  2862.       MoveTempWork(arg, 1);
  2863.       if (! ExtendedMultiply()) {
  2864.          free(temp);
  2865.          return(FALSE);
  2866.          }
  2867.       }
  2868.  
  2869.    free(temp);
  2870.  
  2871.          /* Resolve reciprocal */
  2872.  
  2873.    if (recip) {
  2874.       MoveWorkWork(2, 0);
  2875.       SetWorkInteger(1, 1L);
  2876.       if (! ExtendedDivide()) {
  2877.          return(FALSE);
  2878.          }
  2879.       }
  2880.  
  2881.    return(TRUE);
  2882.  
  2883. }
  2884.  
  2885.  
  2886.  
  2887.  
  2888. /*
  2889.  *    **************************************************
  2890.  *    *                                                *
  2891.  *    *             Extended Reciprocal                *
  2892.  *    *                                                *
  2893.  *    *            work[2] = 1 / work[0]               *
  2894.  *    *                                                *
  2895.  *    **************************************************
  2896.  */
  2897.  
  2898. extern int ExtendedReciprocal(void)
  2899.  
  2900. {
  2901.  
  2902.    SetWorkInteger(1, 1L);        /* 1 -> work[1] */
  2903.  
  2904.    return(ExtendedDivide());
  2905.  
  2906. }
  2907.  
  2908.  
  2909.  
  2910.  
  2911. /*
  2912.  *    **************************************************
  2913.  *    *                                                *
  2914.  *    *              Extended Factorial                *
  2915.  *    *                                                *
  2916.  *    *              work[2] = work[0]!                *
  2917.  *    *                                                *
  2918.  *    **************************************************
  2919.  */
  2920.  
  2921. extern int ExtendedFactorial(void)
  2922.  
  2923. {
  2924.  
  2925.    int i;
  2926.    long x, n;
  2927.  
  2928.  
  2929.          /* Special Cases */
  2930.  
  2931.    if (! work[0].digits) {                /* Test for zero (0! == 1) */
  2932.       SetWorkInteger(2, 1L);
  2933.       return(TRUE);
  2934.       }
  2935.  
  2936.    if (work[0].sign == '-') {             /* Test for negative: error */
  2937.       NegativeArgument();
  2938.       return(FALSE);
  2939.       }
  2940.  
  2941.    if (work[0].exp < work[0].digits) {    /* Test for non integer: error */
  2942.       ArgumentNotInteger();
  2943.       return(FALSE);
  2944.       }
  2945.  
  2946.    if ( (work[0].exp == 1)                /* if X = 1 or 2: X! = X */
  2947.          &&
  2948.         (work[0].man[0] <= 2) ) {
  2949.       MoveWorkWork(0, 2);
  2950.       return(TRUE);
  2951.       }
  2952.  
  2953.    if (work[0].exp > 6) {                 /* Rough test for overflow */
  2954.       Overflow();
  2955.       return(FALSE);
  2956.       }
  2957.  
  2958.  
  2959.    /*
  2960.     *    **************************************************
  2961.     *    *                                                *
  2962.     *    *  FACTORIAL:  X! = 2 * 3 * ... * (X-1) * X      *
  2963.     *    *                                                *
  2964.     *    **************************************************
  2965.     */
  2966.  
  2967.    x = 0;                              /* X -> x */
  2968.    for (i = 0; i < (int)work[0].exp; i++) {
  2969.       x = x * 10 + work[0].man[i];
  2970.       }
  2971.  
  2972.    SetWorkInteger(2, 6L);              /* work[2] = N! (3!) */
  2973.  
  2974.    for (n = 4L; n <= x; n++) {         /* N = 4, ..., X */
  2975.       SetWorkInteger(0, n);            /* work[0] = N */
  2976.       MoveWorkWork(2, 1);              /* work[1] = (N - 1)! */
  2977.       if (! ExtendedMultiply()) {      /* work[2] = N * (N - 1) = N! */
  2978.          return(FALSE);
  2979.          }
  2980.       }
  2981.  
  2982.    return(TRUE);
  2983.  
  2984. }
  2985.