home *** CD-ROM | disk | FTP | other *** search
/ Computerworld 1996 March / Computerworld_1996-03_cd.bin / idg_cd3 / grafika / fraktaly / frasr192 / bignum.c < prev    next >
C/C++ Source or Header  |  1995-03-04  |  40KB  |  1,245 lines
  1. /* bignum.c - C routines for bignumbers */
  2.  
  3. /*
  4. Wesley Loewer's Big Numbers.        (C) 1994, Wesley B. Loewer
  5.  
  6. The bignumber format is simply a signed integer of variable length.  The
  7. bytes are stored in reverse order (least significant byte first, most
  8. significant byte last).  The sign bit is the highest bit of the most
  9. significant byte.  Negatives are stored in 2's complement form.  The byte
  10. length of the bignumbers must be a multiple of 4 for 386+ operations, and
  11. a multiple of 2 for 8086/286 and non 80x86 machines.
  12.  
  13. Some of the arithmetic operations coded here may alter some of the
  14. operands used.  Therefore, take note of the SIDE-EFFECTS listed with each
  15. procedure.  If the value of an operand needs to be retained, just use
  16. copy_bn() first.  This was done for speed sake to avoid unnecessary
  17. copying.  If space is at such a premium that copying it would be
  18. difficult, some of the operations only change the sign of the value.  In
  19. this case, the original could be optained by calling neg_a_bn().
  20.  
  21. Most of the bignumber routines operate on true integers.  Some of the
  22. procedures, however, are designed specifically for fixed decimal point
  23. operations.  The results and how the results are interpreted depend on
  24. where the implied decimal point is located.  The routines that depend on
  25. where the decimal is located are:  strtobn(), bntostr(), bntoint(), inttobn(),
  26. bntofloat(), floattobn(), inv_bn(), div_bn().  The number of bytes
  27. designated for the integer part may be 1, 2, or 4.
  28.  
  29.  
  30. BIGNUMBER FORMAT:
  31. The following is a discription of the bignumber format and associated
  32. variables.  The number is stored in reverse order (Least Significant Byte,
  33. LSB, stored first in memory, Most Significant Byte, MSB, stored last).
  34. Each '_' below represents a block of memory used for arithmetic (1 block =
  35. 4 bytes on 386+, 2 bytes on 286-).  All lengths are given in bytes.
  36.  
  37. LSB                                MSB
  38.   _  _  _  _  _  _  _  _  _  _  _  _
  39. n <----------- bnlength ----------->
  40.                     intlength  ---> <---
  41.  
  42.   bnlength  = the length in bytes of the bignumber
  43.   intlength = the number of bytes used to represent the integer part of
  44.               the bignumber.  Possible values are 1, 2, or 4.  This
  45.               determines the largest number that can be represented by
  46.               the bignumber.
  47.                 intlength = 1, max value = 127.99...
  48.                 intlength = 2, max value = 32,767.99...
  49.                 intlength = 4, max value = 2,147,483,647.99...
  50.  
  51.  
  52. FULL DOUBLE PRECISION MULTIPLICATION:
  53. ( full_mult_bn(), full_square_bn() )
  54.  
  55. The product of two bignumbers, n1 and n2, will be a result, r, which is
  56. a double wide bignumber.  The integer part will also be twice as wide,
  57. thereby eliminating the possiblity of overflowing the number.
  58.  
  59. LSB                                                                    MSB
  60.   _  _  _  _  _  _  _  _  _  _  _  _  _  _  _  _  _  _  _  _  _  _  _  _
  61. r <--------------------------- 2*bnlength ----------------------------->
  62.                                                      2*intlength  --->  <---
  63.  
  64. If this double wide bignumber, r, needs to be converted to a normal,
  65. single width bignumber, this is easily done with pointer arithmetic.  The
  66. converted value starts at r+shiftfactor (where shiftfactor =
  67. bnlength-intlength) and continues for bnlength bytes.  The lower order
  68. bytes and the upper integer part of the double wide number can then be
  69. ignored.
  70.  
  71. LSB                                                                    MSB
  72.   _  _  _  _  _  _  _  _  _  _  _  _  _  _  _  _  _  _  _  _  _  _  _  _
  73. r <--------------------------- 2*bnlength ----------------------------->
  74.                                                      2*intlength  --->  <---
  75.                                  LSB                                  MSB
  76.                r+shiftfactor   <----------  bnlength  ------------>
  77.                                                        intlength  ---> <---
  78.  
  79.  
  80. PARTIAL PRECISION MULTIPLICATION:
  81. ( mult_bn(), square_bn() )
  82.  
  83. In most cases, full double precision multiplication is not necessary.  The
  84. lower order bytes are usually thrown away anyway.  The non-"full"
  85. multiplication routines only calculate rlength bytes in the result.  The
  86. value of rlength must be in the range: 2*bnlength <= rlength < bnlength.
  87. The amount by which rlength exceeds bnlength accounts for the extra bytes
  88. that must be multiplied so that the first bnlength bytes are correct.
  89. These extra bytes are refered to in the code as the "padding," that is:
  90. rlength=bnlength+padding.
  91.  
  92. All three of the values, bnlength, rlength, and therefore padding, must be
  93. multiples of the size of memory blocks being used for arithmetic (2 on
  94. 8086/286 and 4 on 386+).  Typically, the padding is 2*blocksize.  In the
  95. case where bnlength=blocksize, padding can only be blocksize to keep
  96. rlength from being too big.
  97.  
  98. The product of two bignumbers, n1 and n2, will then be a result, r, which
  99. is of length rlength.  The integer part will be twice as wide, thereby
  100. eliminating the possiblity of overflowing the number.
  101.  
  102.              LSB                                      MSB
  103.                _  _  _  _  _  _  _  _  _  _  _  _  _  _
  104.             r  <---- rlength = bnlength+padding ------>
  105.                                     2*intlength  --->  <---
  106.  
  107. If r needs to be converted to a normal, single width bignumber, this is
  108. easily done with pointer arithmetic.  The converted value starts at
  109. r+shiftfactor (where shiftfactor = padding-intlength) and continues for
  110. bnlength bytes.  The lower order bytes and the upper integer part of the
  111. double wide number can then be ignored.
  112.  
  113.              LSB                                      MSB
  114.                _  _  _  _  _  _  _  _  _  _  _  _  _  _
  115.             r  <---- rlength = bnlength+padding ------>
  116.                                     2*intlength  --->  <---
  117.                    LSB                                  MSB
  118.       r+shiftfactor  <----------  bnlength  --------->
  119.                                        intlength ---> <---
  120. */
  121.  
  122. /************************************************************************/
  123. /* There are three parts to the bignum library:                         */
  124. /*                                                                      */
  125. /* 1) bignum.c - initialization, general routines, routines that would  */
  126. /*    not be speeded up much with assembler.                            */
  127. /*                                                                      */
  128. /* 2) bignuma.asm - hand coded assembler routines.                      */
  129. /*                                                                      */
  130. /* 3) bignumc.c - portable C versions of routines in bignuma.asm        */
  131. /*                                                                      */
  132. /************************************************************************/
  133.  
  134. #include <stdio.h>
  135. #include <stdlib.h>
  136. #include <string.h>
  137. #include <malloc.h>
  138. #include "prototyp.h"
  139. #include "bignum.h"
  140. #ifdef sign
  141. #undef sign
  142. #endif
  143.  
  144. /* globals */
  145. extern int bnstep, bnlength, intlength, rlength, padding, shiftfactor, decimals;
  146. extern int bflength, rbflength, bfshiftfactor, bfdecimals;
  147.  
  148. /* used internally by bignum.c routines */
  149. extern bn_t bntmp1, bntmp2, bntmp3, bntmp4, bntmp5, bntmp6; /* rlength  */
  150. extern bn_t bntmpcpy1, bntmpcpy2;                           /* bnlength */
  151.  
  152. /* used by other routines */
  153. extern  bn_t bn_pi;
  154. extern  bn_t bnxmin, bnxmax, bnymin, bnymax, bnx3rd, bny3rd;        /* bnlength */
  155. extern  bn_t bnxdel, bnydel, bnxdel2, bnydel2, bnclosenuff;         /* bnlength */
  156. extern  bn_t bntmpsqrx, bntmpsqry, bntmp;                           /* rlength  */
  157. extern  _BNCMPLX bnold, bnnew, bnparm, bnsaved;                     /* bnlength */
  158.  
  159. #define LOG10_256 2.4082399653118
  160. #define LOG_256   5.5451774444795
  161.  
  162. /*************************************************************************
  163. * The original bignumber code was written specifically for a Little Endian
  164. * system (80x86).  The following is not particularly efficient, but was
  165. * simple to incorporate.  If speed with a Big Endian machine is critical,
  166. * the bignumber format could be reversed.
  167. **************************************************************************/
  168. #ifdef BIG_ENDIAN /* Big Endian: &variable == &(most significant byte), moterola style */
  169. U32 access32(BYTE far *addr)
  170.     {
  171.     return addr[0] | ((U32)addr[1]<<8) | ((U32)addr[2]<<16) | ((U32)addr[3]<<24);
  172.     }
  173.  
  174. U16 access16(BYTE far *addr)
  175.     {
  176.     return (U16)addr[0] | ((U16)addr[1]<<8);
  177.     }
  178.  
  179. S16 accessS16(S16 far *addr)
  180.     {
  181.     return (S16)((BYTE *)addr)[0] | ((S16)((BYTE *)addr)[1]<<8);
  182.     }
  183.  
  184. U32 set32(BYTE far *addr, U32 val)
  185.     {
  186.     addr[0] = val&0xff;
  187.     addr[1] = (val>>8)&0xff;
  188.     addr[2] = (val>>16)&0xff;
  189.     addr[3] = (val>>24)&0xff;
  190.     return val;
  191.     }
  192.  
  193. U16 set16(BYTE far *addr, U16 val)
  194.     {
  195.     addr[0] = val&0xff;
  196.     addr[1] = (val>>8)&0xff;
  197.     return val;
  198.     }
  199.  
  200. S16 setS16(S16 far *addr, S16 val)
  201.     {
  202.     ((BYTE *)addr)[0] = val&0xff;
  203.     ((BYTE *)addr)[1] = (val>>8)&0xff;
  204.     return val;
  205.     }
  206. #endif
  207.  
  208. #if (_MSC_VER >= 700)
  209. #pragma code_seg ("bignum1_text")     /* place following in an overlay */
  210. #endif
  211.  
  212. /************************************************************************/
  213. /* convert_bn  -- convert bignum numbers from old to new lengths        */
  214. int convert_bn(bn_t new, bn_t old, int newbnlength, int newintlength,
  215.                                    int oldbnlength, int oldintlength)
  216.    {
  217.    int savebnlength, saveintlength;
  218.  
  219.    /* save lengths so not dependent on external environment */
  220.    saveintlength = intlength;
  221.    savebnlength  = bnlength;
  222.  
  223.    intlength     = newintlength;
  224.    bnlength      = newbnlength;
  225.    clear_bn(new);
  226.  
  227.    if(newbnlength - newintlength > oldbnlength - oldintlength)
  228.       {
  229.  
  230.       /* This will keep the integer part from overflowing past the array. */
  231.       bnlength = oldbnlength - oldintlength + min(oldintlength, newintlength);
  232.  
  233.       _fmemcpy(new+newbnlength-newintlength-oldbnlength+oldintlength,
  234.                old,bnlength);
  235.       }
  236.    else
  237.       {
  238.       bnlength = newbnlength - newintlength + min(oldintlength, newintlength);
  239.       _fmemcpy(new,old+oldbnlength-oldintlength-newbnlength+newintlength,
  240.                bnlength);
  241.       }
  242.    intlength = saveintlength;
  243.    bnlength  = savebnlength;
  244.    return(0);
  245.    }
  246.  
  247. /********************************************************************/
  248. /* bn_hexdump() - for debugging, dumps to stdout                    */
  249.  
  250. void bn_hexdump(bn_t r)
  251.     {
  252.     int i;
  253.  
  254.     for (i=0; i<bnlength; i++)
  255.         printf("%02X ", r[i]);
  256.     printf("\n");
  257.     return;
  258.     }
  259.  
  260. /**********************************************************************/
  261. /* strtobn() - converts a string into a bignumer                       */
  262. /*   r - pointer to a bignumber                                       */
  263. /*   s - string in the floating point format [-][digits].[digits]     */
  264. /*   note: the string may not be empty or have extra space and may    */
  265. /*   not use scientific notation (2.3e4).                             */
  266.  
  267. bn_t strtobn(bn_t r, char *s)
  268.     {
  269.     unsigned l;
  270.     bn_t onesbyte;
  271.     int signflag=0;
  272.     long longval;
  273.  
  274.     clear_bn(r);
  275.     onesbyte = r + bnlength - intlength;
  276.  
  277.     if (s[0] == '+')    /* for + sign */
  278.         {
  279.         s++;
  280.         }
  281.     else if (s[0] == '-')    /* for neg sign */
  282.         {
  283.         signflag = 1;
  284.         s++;
  285.         }
  286.  
  287.     if (strchr(s, '.') != NULL) /* is there a decimal point? */
  288.         {
  289.         l = strlen(s) - 1;      /* start with the last digit */
  290.         while (s[l] >= '0' && s[l] <= '9') /* while a digit */
  291.             {
  292.             *onesbyte = (BYTE)(s[l--] - '0');
  293.             div_a_bn_int(r, 10);
  294.             }
  295.  
  296.         if (s[l] == '.')
  297.             {
  298.             longval = atol(s);
  299.             switch (intlength)
  300.                 { /* only 1, 2, or 4 are allowed */
  301.                 case 1:
  302.                     *onesbyte = (BYTE)longval;
  303.                     break;
  304.                 case 2:
  305.                     set16(onesbyte, (U16)longval);
  306.                     break;
  307.                 case 4:
  308.                     set32(onesbyte, longval);
  309.                     break;
  310.                 }
  311.             }
  312.         }
  313.     else
  314.         {
  315.         longval = atol(s);
  316.         switch (intlength)
  317.             { /* only 1, 2, or 4 are allowed */
  318.             case 1:
  319.                 *onesbyte = (BYTE)longval;
  320.                 break;
  321.             case 2:
  322.                 set16(onesbyte, (U16)longval);
  323.                 break;
  324.             case 4:
  325.                 set32(onesbyte, longval);
  326.                 break;
  327.             }
  328.         }
  329.  
  330.  
  331.     if (signflag)
  332.         neg_a_bn(r);
  333.  
  334.     return r;
  335.     }
  336.  
  337. /********************************************************************/
  338. /* strlen_needed() - returns string length needed to hold bignumber */
  339.  
  340. int strlen_needed()
  341.    {
  342.    int length;
  343.  
  344.    /* first space for integer part */
  345.    switch(intlength)
  346.       {
  347.       case 1:
  348.          length = 3;  /* max 127 */
  349.          break;
  350.       case 2:
  351.          length = 5;  /* max 32767 */
  352.          break;
  353.       case 4:
  354.          length = 10; /* max 2147483647 */
  355.          break;
  356.       }
  357.     length += decimals;  /* decimal part */
  358.     length += 2;         /* decimal point and sign */
  359.     length += 4;         /* null and a little extra for safety */
  360.     return(length);
  361.     }
  362.          
  363. /********************************************************************/
  364. /* bntostr() - converts a bignumber into a string                    */
  365. /*   s - string, must be large enough to hold the number.           */
  366. /*   r - bignumber                                                  */
  367. /*   will covert to a floating point notation                       */
  368. /*   SIDE-EFFECT: the bignumber, r, is destroyed.                   */
  369. /*                Copy it first if necessary.                       */
  370.  
  371. char *unsafe_bntostr(char *s, int dec, bn_t r)
  372.     {
  373.     int l=0, d;
  374.     bn_t onesbyte;
  375.     long longval;
  376.  
  377.     if (dec == 0)
  378.         dec = decimals;
  379.     onesbyte = r + bnlength - intlength;
  380.  
  381.     if (is_bn_neg(r))
  382.         {
  383.         neg_a_bn(r);
  384.         *(s++) = '-';
  385.         }
  386.     switch (intlength)
  387.         { /* only 1, 2, or 4 are allowed */
  388.         case 1:
  389.             longval = *onesbyte;
  390.             break;
  391.         case 2:
  392.             longval = access16(onesbyte);
  393.             break;
  394.         case 4:
  395.             longval = access32(onesbyte);
  396.             break;
  397.         }
  398.     ltoa(longval, s, 10);
  399.     l = strlen(s);
  400.     s[l++] = '.';
  401.     for (d=0; d < dec; d++)
  402.         {
  403.         *onesbyte = 0;  /* clear out highest byte */
  404.         mult_a_bn_int(r, 10);
  405.         if (is_bn_zero(r))
  406.             break;
  407.         s[l++] = (BYTE)(*onesbyte + '0');
  408.         }
  409.     s[l] = '\0'; /* don't forget nul char */
  410.  
  411.     return s;
  412.     }
  413.  
  414. #if (_MSC_VER >= 700)
  415. #pragma code_seg ( )         /* back to normal segment */
  416. #endif
  417.  
  418. /*********************************************************************/
  419. /*  b = l                                                            */
  420. /*  Converts a long to a bignumber                                   */
  421. bn_t inttobn(bn_t r, long longval)
  422.     {
  423.     bn_t onesbyte;
  424.  
  425.     clear_bn(r);
  426.     onesbyte = r + bnlength - intlength;
  427.     switch (intlength)
  428.         { /* only 1, 2, or 4 are allowed */
  429.         case 1:
  430.             *onesbyte = (BYTE)longval;
  431.             break;
  432.         case 2:
  433.             set16(onesbyte, (U16)longval);
  434.             break;
  435.         case 4:
  436.             set32(onesbyte, longval);
  437.             break;
  438.         }
  439.     return r;
  440.     }
  441.  
  442. /*********************************************************************/
  443. /*  l = floor(b), floor rounds down                                  */
  444. /*  Converts the integer part a bignumber to a long                  */
  445. long bntoint(bn_t n)
  446.     {
  447.     bn_t onesbyte;
  448.     long longval;
  449.  
  450.     onesbyte = n + bnlength - intlength;
  451.     switch (intlength)
  452.         { /* only 1, 2, or 4 are allowed */
  453.         case 1:
  454.             longval = *onesbyte;
  455.             break;
  456.         case 2:
  457.             longval = access16(onesbyte);
  458.             break;
  459.         case 4:
  460.             longval = access32(onesbyte);
  461.             break;
  462.         }
  463.     return longval;
  464.     }
  465.  
  466. /*********************************************************************/
  467. /*  b = f                                                            */
  468. /*  Converts a double to a bignumber                                 */
  469. bn_t floattobn(bn_t r, LDBL f)
  470.     {
  471. #ifndef USE_BIGNUM_C_CODE
  472.     floattobf(bftmp1, f);
  473.     return bftobn(r, bftmp1);
  474.     
  475. #else
  476.     bn_t onesbyte;
  477.     int i;
  478.     int signflag=0;
  479.  
  480.     clear_bn(r);
  481.     onesbyte = r + bnlength - intlength;
  482.  
  483.     if (f < 0)
  484.         {
  485.         signflag = 1;
  486.         f = -f;
  487.         }
  488.  
  489.     switch (intlength)
  490.         { /* only 1, 2, or 4 are allowed */
  491.         case 1:
  492.             *onesbyte = (BYTE)f;
  493.             break;
  494.         case 2:
  495.             set16(onesbyte, (U16)f);
  496.             break;
  497.         case 4:
  498.             set32(onesbyte, (U32)f);
  499.             break;
  500.         }
  501.  
  502.     f -= (long)f; /* keep only the decimal part */
  503.     for (i = bnlength - intlength - 1; i >= 0 && f != 0.0; i--)
  504.         {
  505.         f *= 256;
  506.         r[i] = (BYTE)f;  /* keep use the integer part */
  507.         f -= (BYTE)f; /* now throw away the integer part */
  508.         }
  509.  
  510.     if (signflag)
  511.         neg_a_bn(r);
  512.  
  513.     return r;
  514. #endif /* USE_BIGNUM_C_CODE */
  515.     }
  516.  
  517. /********************************************************************/
  518. /* sign(r)                                                          */
  519. int sign_bn(bn_t n)
  520.     {
  521.     return is_bn_neg(n) ? -1 : is_bn_not_zero(n) ? 1 : 0;
  522.     }
  523.  
  524. /********************************************************************/
  525. /* r = |n|                                                          */
  526. bn_t abs_bn(bn_t r, bn_t n)
  527.     {
  528.     copy_bn(r,n);    
  529.     if (is_bn_neg(r))
  530.         neg_a_bn(r);
  531.     return r;
  532.     }
  533.  
  534. /********************************************************************/
  535. /* r = |r|                                                          */
  536. bn_t abs_a_bn(bn_t r)
  537.     {
  538.     if (is_bn_neg(r))
  539.         neg_a_bn(r);
  540.     return r;
  541.     }
  542.  
  543. /********************************************************************/
  544. /* r = 1/n                                                          */
  545. /* uses bntmp1 - bntmp3 - global temp bignumbers                    */
  546. /*  SIDE-EFFECTS:                                                   */
  547. /*      n ends up as |n|    Make copy first if necessary.           */
  548. bn_t unsafe_inv_bn(bn_t r, bn_t n)
  549.     {
  550.     int signflag=0, i;
  551.     long maxval;
  552.     LDBL f;
  553.  
  554.     /* use Newton's recursive method for zeroing in on 1/n : r=r(2-rn) */
  555.  
  556.     if (is_bn_neg(n))
  557.         { /* will be a lot easier to deal with just positives */
  558.         signflag = 1;
  559.         neg_a_bn(n);
  560.         }
  561.  
  562.     f = bntofloat(n);
  563.     if (f == 0) /* division by zero */
  564.         {
  565.         max_bn(r);
  566.         return r;
  567.         }
  568.     f = 1/f; /* approximate inverse */
  569.     maxval = (1L << ((intlength<<3)-1)) - 1;
  570.     if (f > maxval) /* check for overflow */
  571.         {
  572.         max_bn(r);
  573.         return r;
  574.         }
  575.     else if (f <= -maxval)
  576.         {
  577.         max_bn(r);
  578.         neg_a_bn(r);
  579.         return r;
  580.         }
  581.  
  582.     floattobn(r, f); /* start with approximate inverse */
  583.     clear_bn(bntmp2); /* will be used as 1.0 and 2.0 */
  584.  
  585.     for (i=0; i<25; i++) /* safety net, this shouldn't ever be needed */
  586.         {
  587.         unsafe_mult_bn(bntmp1, r, n); /* bntmp1=rn */
  588.         bntmp2[bnlength-intlength] = 1; /* bntmp2 = 1.0 */
  589.         if ( !cmp_bn(bntmp2, bntmp1+shiftfactor) ) /* if not different */
  590.             break;  /* they must be the same */
  591.         bntmp2[bnlength-intlength] = 2; /* bntmp2 = 2.0 */
  592.         sub_bn(bntmp3, bntmp2, bntmp1+shiftfactor); /* bntmp3=2-rn */
  593.         unsafe_mult_bn(bntmp1, r, bntmp3); /* bntmp1=r(2-rn) */
  594.         copy_bn(r, bntmp1+shiftfactor); /* r = bntmp1 */
  595.         }
  596.  
  597.     if (signflag)
  598.         {
  599.         neg_a_bn(r);
  600.         }
  601.     return r;
  602.     }
  603.  
  604. /********************************************************************/
  605. /* r = n1/n2                                                        */
  606. /*      r - result of length bnlength                               */
  607. /* uses bntmp1 - bntmp3 - global temp bignumbers                    */
  608. /*  SIDE-EFFECTS:                                                   */
  609. /*      n1, n2 end up as |n1|, |n2|                                 */
  610. /*      Make copies first if necessary.                             */
  611. bn_t unsafe_div_bn(bn_t r, bn_t n1, bn_t n2)
  612.     {
  613.     int scale1, scale2, scale, sign=0, i;
  614.     long maxval;
  615.     LDBL a, b, f;
  616.  
  617.     /* first, check for valid data */
  618.     a = bntofloat(n1);
  619.     if (a == 0) /* division into zero */
  620.         {
  621.         clear_bn(r); /* return 0 */
  622.         return r;
  623.         }
  624.     b = bntofloat(n2);
  625.     if (b == 0) /* division by zero */
  626.         {
  627.         max_bn(r);
  628.         return r;
  629.         }
  630.     f = a/b; /* approximate quotient */
  631.     maxval = (1L << ((intlength<<3)-1)) - 1;
  632.     if (f > maxval) /* check for overflow */
  633.         {
  634.         max_bn(r);
  635.         return r;
  636.         }
  637.     else if (f <= -maxval)
  638.         {
  639.         max_bn(r);
  640.         neg_a_bn(r);
  641.         return r;
  642.         }
  643.     /* appears to be ok, do division */
  644.  
  645.     if (is_bn_neg(n1))
  646.         {
  647.         neg_a_bn(n1);
  648.         sign = !sign;
  649.         }
  650.     if (is_bn_neg(n2))
  651.         {
  652.         neg_a_bn(n2);
  653.         sign = !sign;
  654.         }
  655.  
  656.     /* scale n1 and n2 so: |n| >= 1/256 */
  657.     /* scale = (int)(log(1/fabs(a))/LOG_256) = LOG_256(1/|a|) */
  658.     i = bnlength-1;
  659.     while (i >= 0 && n1[i] == 0)
  660.         i--;
  661.     scale1 = bnlength - i - 2;
  662.     if (scale1 < 0)
  663.        scale1 = 0;
  664.     i = bnlength-1;
  665.     while (i >= 0 && n2[i] == 0)
  666.         i--;
  667.     scale2 = bnlength - i - 2;
  668.     if (scale2 < 0)
  669.        scale2 = 0;
  670.  
  671.     /* important!, use memmove(), not memcpy() */
  672.     _fmemmove(n1+scale1, n1, bnlength-scale1); /* shift bytes over */
  673.     _fmemset(n1, 0, scale1);  /* zero out the rest */
  674.     _fmemmove(n2+scale2, n2, bnlength-scale2); /* shift bytes over */
  675.     _fmemset(n2, 0, scale2);  /* zero out the rest */
  676.  
  677.     unsafe_inv_bn(r, n2);
  678.     unsafe_mult_bn(bntmp1, n1, r);
  679.     copy_bn(r, bntmp1+shiftfactor); /* r = bntmp1 */
  680.  
  681.     if (scale1 != scale2)
  682.         {
  683.         /* Rescale r back to what it should be.  Overflow has already been checked */
  684.         if (scale1 > scale2) /* answer is too big, adjust it */
  685.             {
  686.             scale = scale1-scale2;
  687.             _fmemmove(r, r+scale, bnlength-scale); /* shift bytes over */
  688.             _fmemset(r+bnlength-scale, 0, scale);  /* zero out the rest */
  689.             }
  690.         else if (scale1 < scale2) /* answer is too small, adjust it */
  691.             {
  692.             scale = scale2-scale1;
  693.             _fmemmove(r+scale, r, bnlength-scale); /* shift bytes over */
  694.             _fmemset(r, 0, scale);                 /* zero out the rest */
  695.             }
  696.         /* else scale1 == scale2 */
  697.  
  698.         }
  699.  
  700.     if (sign)
  701.         neg_a_bn(r);
  702.  
  703.     return r;
  704.     }
  705.  
  706. /********************************************************************/
  707. /* sqrt(r)                                                          */
  708. /* uses bntmp1 - bntmp4 - global temp bignumbers                    */
  709. /*  SIDE-EFFECTS:                                                   */
  710. /*      n ends up as |n|                                            */
  711. bn_t sqrt_bn(bn_t r, bn_t n)
  712.     {
  713.     int i, comp, almost_match=0;
  714.     LDBL f;
  715.  
  716. /* use Newton's recursive method for zeroing in on sqrt(n): r=.5(r+n/r) */
  717.  
  718.     if (is_bn_neg(n))
  719.         { /* sqrt of a neg, return 0 */
  720.         clear_bn(r);
  721.         return r;
  722.         }
  723.  
  724.     f = bntofloat(n);
  725.     if (f == 0) /* division by zero will occur */
  726.         {
  727.         clear_bn(r); /* sqrt(0) = 0 */
  728.         return r;
  729.         }
  730.     f = sqrtl(f); /* approximate square root */
  731.     /* no need to check overflow */
  732.  
  733.     floattobn(r, f); /* start with approximate sqrt */
  734.     copy_bn(bntmp4, r);
  735.  
  736.     for (i=0; i<25; i++) /* safety net, this shouldn't ever be needed */
  737.         {
  738.         unsafe_div_bn(bntmp4, n, r);
  739.         add_a_bn(r, bntmp4);
  740.         half_a_bn(r);
  741.         if ((comp=abs(cmp_bn(r, bntmp4))) <= 2 ) /* if match or almost match */
  742.             {
  743.             if (comp == 0)  /* perfect match */
  744.                 break;
  745.             if (almost_match == 2) /* did they almost match before? */
  746.                 { /* yes, this must be the third time */
  747.                 /* average last 4 bytes */
  748.                 *(U32 far *)r = (*(U32 far *)r + *(U32 far *)bntmp4) / 2;
  749.                 break;
  750.                 }
  751.             else /* this is the first or second time they almost matched */
  752.                 almost_match++;
  753.             }
  754.         }
  755.  
  756.     return r;
  757.     }
  758.  
  759. /********************************************************************/
  760. /* exp(r)                                                           */
  761. /* uses bntmp1, bntmp2, bntmp3 - global temp bignumbers             */
  762. bn_t exp_bn(bn_t r, bn_t n)
  763.     {
  764.     U16 fact=1;
  765.  
  766.     if (is_bn_zero(n))
  767.         {
  768.         inttobn(r, 1);
  769.         return r;
  770.         }
  771.  
  772. /* use Taylor Series (very slow convergence) */
  773.     inttobn(r, 1); /* start with r=1.0 */
  774.     copy_bn(bntmp2, r);
  775.     for (;;)
  776.         {
  777.         /* copy n, if n is negative, mult_bn() alters n */
  778.         unsafe_mult_bn(bntmp3, bntmp2, copy_bn(bntmp1, n));
  779.         copy_bn(bntmp2, bntmp3+shiftfactor);
  780.         div_a_bn_int(bntmp2, fact);
  781.         if (!is_bn_not_zero(bntmp2))
  782.             break; /* too small to register */
  783.         add_a_bn(r, bntmp2);
  784.         fact++;
  785.         }
  786.     return r;
  787.     }
  788.  
  789. /********************************************************************/
  790. /* ln(r)                                                            */
  791. /* uses bntmp1 - bntmp6 - global temp bignumbers                    */
  792. /*  SIDE-EFFECTS:                                                   */
  793. /*      n ends up as |n|                                            */
  794. bn_t unsafe_ln_bn(bn_t r, bn_t n)
  795.     {
  796.     int i, comp, almost_match=0;
  797.     long maxval;
  798.     LDBL f;
  799.  
  800. /* use Newton's recursive method for zeroing in on ln(n): r=r+n*exp(-r)-1 */
  801.  
  802.     if (is_bn_neg(n) || is_bn_zero(n))
  803.         { /* error, return largest neg value */
  804.         max_bn(r);
  805.         neg_a_bn(r);
  806.         return r;
  807.         }
  808.  
  809.     f = bntofloat(n);
  810.     f = logl(f); /* approximate ln(x) */
  811.     maxval = (1L << ((intlength<<3)-1)) - 1;
  812.     if (f > maxval) /* check for overflow */
  813.         {
  814.         max_bn(r);
  815.         return r;
  816.         }
  817.     else if (f <= -maxval)
  818.         {
  819.         max_bn(r);
  820.         neg_a_bn(r);
  821.         return r;
  822.         }
  823.     /* appears to be ok, do ln */
  824.  
  825.     floattobn(r, f); /* start with approximate ln */
  826.     neg_a_bn(r); /* -r */
  827.     copy_bn(bntmp5, r); /* -r */
  828.     inttobn(bntmp4, 1);
  829.  
  830.     for (i=0; i<25; i++) /* safety net, this shouldn't ever be needed */
  831.         {
  832.         exp_bn(bntmp6, r);     /* exp(-r) */
  833.         unsafe_mult_bn(bntmp2, bntmp6, n);  /* n*exp(-r) */
  834.         sub_a_bn(bntmp2+shiftfactor, bntmp4);   /* n*exp(-r) - 1 */
  835.         sub_a_bn(r, bntmp2+shiftfactor);        /* -r - (n*exp(-r) - 1) */
  836.  
  837.         if ((comp=abs(cmp_bn(r, bntmp5))) <= 2 ) /* if match or almost match */
  838.             {
  839.             if (comp == 0)  /* perfect match */
  840.                 break;
  841.             if (almost_match == 2) /* did they almost match twice before? */
  842.                 { /* yes, this must be the third time */
  843.                 /* average last 4 bytes */
  844.                 *(U32 far *)r = (*(U32 far *)r + *(U32 far *)bntmp5) / 2;
  845.                 break;
  846.                 }
  847.             else /* this is the first or second time they almost matched */
  848.                 almost_match++;
  849.             }
  850.         copy_bn(bntmp5, r); /* -r */
  851.         }
  852.     neg_a_bn(r); /* -(-r) */
  853.     return r;
  854.     }
  855.  
  856. /********************************************************************/
  857. /* sincos_bn(r)                                                     */
  858. /* uses bntmp1 - bntmp2 - global temp bignumbers                    */
  859. /*  SIDE-EFFECTS:                                                   */
  860. /*      n ends up as |n| mod (pi/4)                                 */
  861. bn_t unsafe_sincos_bn(bn_t s, bn_t c, bn_t n)
  862.     {
  863.     U16 fact=2;
  864.     int k=0, i, halves;
  865.     int signcos=0, signsin=0, switch_sincos=0;
  866.  
  867.     /* assure range 0 <= x < pi/4 */
  868.  
  869.     if (is_bn_zero(n))
  870.         {
  871.         clear_bn(s);    /* sin(0) = 0 */
  872.         inttobn(c, 1);  /* cos(0) = 1 */
  873.         return s;
  874.         }
  875.  
  876.     if (is_bn_neg(n))
  877.         {
  878.         signsin = !signsin; /* sin(-x) = -sin(x), odd; cos(-x) = cos(x), even */
  879.         neg_a_bn(n);
  880.         }
  881.     /* n >= 0 */
  882.  
  883.     double_bn(bntmp1, bn_pi); /* 2*pi */
  884.     /* this could be done with remainders, but it would probably be slower */
  885.     while (cmp_bn(n, bntmp1) >= 0) /* while n >= 2*pi */
  886.         sub_a_bn(n, bntmp1);
  887.     /* 0 <= n < 2*pi */
  888.  
  889.     copy_bn(bntmp1, bn_pi); /* pi */
  890.     if (cmp_bn(n, bntmp1) >= 0) /* if n >= pi */
  891.         {
  892.         sub_a_bn(n, bntmp1);
  893.         signsin = !signsin;
  894.         signcos = !signcos;
  895.         }
  896.     /* 0 <= n < pi */
  897.  
  898.     half_bn(bntmp1, bn_pi); /* pi/2 */
  899.     if (cmp_bn(n, bntmp1) > 0) /* if n > pi/2 */
  900.         {
  901.         sub_bn(n, bn_pi, n);   /* pi - n */
  902.         signcos = !signcos;
  903.         }
  904.     /* 0 <= n < pi/2 */
  905.  
  906.     half_bn(bntmp1, bn_pi); /* pi/2 */
  907.     half_a_bn(bntmp1);      /* pi/4 */
  908.     if (cmp_bn(n, bntmp1) > 0) /* if n > pi/4 */
  909.         {
  910.         half_bn(bntmp1, bn_pi); /* pi/2 */
  911.         sub_bn(n, bntmp1, n);  /* pi/2 - n */
  912.         switch_sincos = !switch_sincos;
  913.         }
  914.     /* 0 <= n < pi/4 */
  915.  
  916.     /* this looks redundant, but n could now be zero when it wasn't before */
  917.     if (is_bn_zero(n))
  918.         {
  919.         clear_bn(s);    /* sin(0) = 0 */
  920.         inttobn(c, 1);  /* cos(0) = 1 */
  921.         return s;
  922.         }
  923.  
  924.  
  925. /* at this point, the double angle trig identities could be used as many  */
  926. /* times as desired to reduce the range to pi/8, pi/16, etc...  Each time */
  927. /* the range is cut in half, the number of iterations required is reduced */
  928. /* by "quite a bit."  It's just a matter of testing to see what gives the */
  929. /* optimal results.                                                       */
  930.    /* halves = bnlength / 10; */ /* this is experimental */
  931.    halves = 1;
  932.    for (i = 0; i < halves; i++)
  933.        half_a_bn(n);
  934.  
  935. /* use Taylor Series (very slow convergence) */
  936.     copy_bn(s, n); /* start with s=n */
  937.     inttobn(c, 1); /* start with c=1 */
  938.     copy_bn(bntmp1, n); /* the current x^n/n! */
  939.  
  940.     for (;;)
  941.         {
  942.         /* even terms for cosine */
  943.         unsafe_mult_bn(bntmp2, bntmp1, n);
  944.         copy_bn(bntmp1, bntmp2+shiftfactor);
  945.         div_a_bn_int(bntmp1, fact++);
  946.         if (!is_bn_not_zero(bntmp1))
  947.             break; /* too small to register */
  948.         if (k) /* alternate between adding and subtracting */
  949.             add_a_bn(c, bntmp1);
  950.         else
  951.             sub_a_bn(c, bntmp1);
  952.  
  953.         /* odd terms for sine */
  954.         unsafe_mult_bn(bntmp2, bntmp1, n);
  955.         copy_bn(bntmp1, bntmp2+shiftfactor);
  956.         div_a_bn_int(bntmp1, fact++);
  957.         if (!is_bn_not_zero(bntmp1))
  958.             break; /* too small to register */
  959.         if (k) /* alternate between adding and subtracting */
  960.             add_a_bn(s, bntmp1);
  961.         else
  962.             sub_a_bn(s, bntmp1);
  963.         k = !k; /* toggle */
  964.         }
  965.  
  966.      /* now need to undo what was done by cutting angles in half */
  967.      inttobn(bntmp1, 1);
  968.      for (i = 0; i < halves; i++)
  969.          {
  970.          unsafe_mult_bn(bntmp2, s, c); /* no need for safe mult */
  971.          double_bn(s, bntmp2+shiftfactor); /* sin(2x) = 2*sin(x)*cos(x) */
  972.          unsafe_square_bn(bntmp2,c);
  973.          double_a_bn(bntmp2+shiftfactor);
  974.          sub_bn(c, bntmp2+shiftfactor, bntmp1); /* cos(2x) = 2*cos(x)*cos(x) - 1 */
  975.          }
  976.  
  977.     if (switch_sincos)
  978.          {
  979.          copy_bn(bntmp1, s);
  980.          copy_bn(s, c);
  981.          copy_bn(c, bntmp1);
  982.          }
  983.      if (signsin)
  984.          neg_a_bn(s);
  985.      if (signcos)
  986.          neg_a_bn(c);
  987.      return s; /* return sine I guess */
  988.     }
  989.  
  990. /********************************************************************/
  991. /* atan(r)                                                          */
  992. /* uses bntmp1 - bntmp5 - global temp bignumbers                    */
  993. /*  SIDE-EFFECTS:                                                   */
  994. /*      n ends up as |n|                                            */
  995. bn_t unsafe_atan_bn(bn_t r, bn_t n)
  996.     {
  997.     int i, comp, almost_match=0, signflag=0;
  998.     LDBL f;
  999.  
  1000. /* use Newton's recursive method for zeroing in on atan(n): r=r-cos(r)(sin(r)-n*cos(r)) */
  1001.  
  1002.     if (is_bn_neg(n))
  1003.         {
  1004.         signflag = 1;
  1005.         neg_a_bn(n);
  1006.         }
  1007.  
  1008.     f = bntofloat(n);
  1009.     f = atanl(f); /* approximate arctangent */
  1010.     /* no need to check overflow */
  1011.  
  1012.     floattobn(r, f); /* start with approximate atan */
  1013.     copy_bn(bntmp3, r);
  1014.  
  1015.     for (i=0; i<25; i++) /* safety net, this shouldn't ever be needed */
  1016.         {
  1017.         unsafe_sincos_bn(bntmp4, bntmp5, bntmp3);   /* sin(r), cos(r) */
  1018.         copy_bn(bntmp3, r); /* restore bntmp3 from sincos_bn() */
  1019.         copy_bn(bntmp1, bntmp5);
  1020.         unsafe_mult_bn(bntmp2, n, bntmp1);     /* n*cos(r) */
  1021.         sub_a_bn(bntmp4, bntmp2+shiftfactor); /* sin(r) - n*cos(r) */
  1022.         unsafe_mult_bn(bntmp1, bntmp5, bntmp4); /* cos(r) * (sin(r) - n*cos(r)) */
  1023.         sub_a_bn(r, bntmp1+shiftfactor); /* r - cos(r) * (sin(r) - n*cos(r)) */
  1024.  
  1025.         if ((comp=abs(cmp_bn(r, bntmp3))) <= 2 ) /* if match or almost match */
  1026.             {
  1027.             if (comp == 0)  /* perfect match */
  1028.                 break;
  1029.             if (almost_match == 2) /* did they almost match before? */
  1030.                 { /* yes, this must be the third time */
  1031.                 /* average last 4 bytes */
  1032.                 *(U32 far *)r = (*(U32 far *)r + *(U32 far *)bntmp3) / 2;
  1033.                 break;
  1034.                 }
  1035.             else /* this is the first or second time they almost matched */
  1036.                 almost_match++;
  1037.             }
  1038.         copy_bn(bntmp3, r); /* make a copy for later comparison */
  1039.         }
  1040.  
  1041.     if (signflag)
  1042.         neg_a_bn(r);
  1043.     return r;
  1044.     }
  1045.  
  1046. /********************************************************************/
  1047. /* atan2(r,ny,nx)                                                     */
  1048. /* uses bntmp1 - bntmp6 - global temp bigfloats                     */
  1049. bn_t unsafe_atan2_bn(bn_t r, bn_t ny, bn_t nx)
  1050.    {
  1051.    int signx, signy;
  1052.  
  1053.    signx = sign_bn(nx);
  1054.    signy = sign_bn(ny);
  1055.  
  1056.    if (signy == 0)
  1057.       {
  1058.       if (signx < 0)
  1059.          copy_bn(r, bn_pi); /* negative x axis, 180 deg */
  1060.       else    /* signx >= 0    positive x axis, 0 */
  1061.          clear_bn(r);
  1062.       return(r);
  1063.       }
  1064.    if (signx == 0)
  1065.       {
  1066.       copy_bn(r, bn_pi); /* y axis */
  1067.       half_a_bn(r);      /* +90 deg */
  1068.       if (signy < 0)
  1069.          neg_a_bn(r);    /* -90 deg */
  1070.       return(r);
  1071.       }
  1072.  
  1073.    if (signy < 0)
  1074.       neg_a_bn(ny);
  1075.    if (signx < 0)
  1076.       neg_a_bn(nx);
  1077.    unsafe_div_bn(bntmp6,ny,nx);
  1078.    unsafe_atan_bn(r, bntmp6);
  1079.    if (signx < 0)
  1080.       sub_bn(r,bn_pi,r);
  1081.    if(signy < 0)
  1082.       neg_a_bn(r);
  1083.    return(r);
  1084.    }
  1085.  
  1086.  
  1087. /**********************************************************************/
  1088. /* The rest of the functions are "safe" versions of the routines that */
  1089. /* have side effects which alter the parameters                       */
  1090. /**********************************************************************/
  1091.  
  1092. /**********************************************************************/
  1093. bn_t full_mult_bn(bn_t r, bn_t n1, bn_t n2)
  1094.     {
  1095.     int sign1, sign2;
  1096.  
  1097.     sign1 = is_bn_neg(n1);
  1098.     sign2 = is_bn_neg(n2);
  1099.     unsafe_full_mult_bn(r, n1, n2);
  1100.     if (sign1)
  1101.         neg_a_bn(n1);
  1102.     if (sign2)
  1103.         neg_a_bn(n2);
  1104.     return r;
  1105.     }
  1106.  
  1107. /**********************************************************************/
  1108. bn_t mult_bn(bn_t r, bn_t n1, bn_t n2)
  1109.     {
  1110.     int sign1, sign2;
  1111.  
  1112.     /* TW ENDFIX */   
  1113.     sign1 = is_bn_neg(n1);
  1114.     sign2 = is_bn_neg(n2);
  1115.     unsafe_mult_bn(r, n1, n2);
  1116.     if (sign1)
  1117.         neg_a_bn(n1);
  1118.     if (sign2)
  1119.         neg_a_bn(n2);
  1120.     return r;
  1121.     }
  1122.  
  1123. /**********************************************************************/
  1124. bn_t full_square_bn(bn_t r, bn_t n)
  1125.     {
  1126.     int sign;
  1127.  
  1128.     sign = is_bn_neg(n);
  1129.     unsafe_full_square_bn(r, n);
  1130.     if (sign)
  1131.         neg_a_bn(n);
  1132.     return r;
  1133.     }
  1134.  
  1135. /**********************************************************************/
  1136. bn_t square_bn(bn_t r, bn_t n)
  1137.     {
  1138.     int sign;
  1139.  
  1140.     sign = is_bn_neg(n);
  1141.     unsafe_square_bn(r, n);
  1142.     if (sign)
  1143.         neg_a_bn(n);
  1144.     return r;
  1145.     }
  1146.  
  1147. /**********************************************************************/
  1148. bn_t div_bn_int(bn_t r, bn_t n, U16 u)
  1149.     {
  1150.     int sign;
  1151.  
  1152.     sign = is_bn_neg(n);
  1153.     unsafe_div_bn_int(r, n, u);
  1154.     if (sign)
  1155.         neg_a_bn(n);
  1156.     return r;
  1157.     }
  1158.  
  1159. #if (_MSC_VER >= 700)
  1160. #pragma code_seg ("bignum1_text")     /* place following in an overlay */
  1161. #endif
  1162.  
  1163. /**********************************************************************/
  1164. char *bntostr(char *s, int dec, bn_t r)
  1165.     {
  1166.     return unsafe_bntostr(s, dec, copy_bn(bntmpcpy2, r));
  1167.     }
  1168.  
  1169. #if (_MSC_VER >= 700)
  1170. #pragma code_seg ( )        /* back to normal segment */
  1171. #endif
  1172.  
  1173. /**********************************************************************/
  1174. bn_t inv_bn(bn_t r, bn_t n)
  1175.     {
  1176.     int sign;
  1177.  
  1178.     sign = is_bn_neg(n);
  1179.     unsafe_inv_bn(r, n);
  1180.     if (sign)
  1181.         neg_a_bn(n);
  1182.     return r;
  1183.     }
  1184.  
  1185. /**********************************************************************/
  1186. bn_t div_bn(bn_t r, bn_t n1, bn_t n2)
  1187.     {
  1188.     copy_bn(bntmpcpy1, n1);
  1189.     copy_bn(bntmpcpy2, n2);
  1190.     return unsafe_div_bn(r, bntmpcpy1, bntmpcpy2);
  1191.     }
  1192.  
  1193. /**********************************************************************/
  1194. bn_t ln_bn(bn_t r, bn_t n)
  1195.     {
  1196. #if 0
  1197.     int sign;
  1198.  
  1199.     sign = is_bn_neg(n);
  1200.     unsafe_ln_bn(r, n);
  1201.     if (sign)
  1202.         neg_a_bn(n);
  1203. #endif
  1204.     copy_bn(bntmpcpy1, n); /* allows r and n to overlap memory */
  1205.     unsafe_ln_bn(r, bntmpcpy1);
  1206.     return r;
  1207.     }
  1208.  
  1209. /**********************************************************************/
  1210. bn_t sincos_bn(bn_t s, bn_t c, bn_t n)
  1211.     {
  1212.     return unsafe_sincos_bn(s, c, copy_bn(bntmpcpy1, n));
  1213.     }
  1214.  
  1215. /**********************************************************************/
  1216. bn_t atan_bn(bn_t r, bn_t n)
  1217.     {
  1218.     int sign;
  1219.  
  1220.     sign = is_bn_neg(n);
  1221.     unsafe_atan_bn(r, n);
  1222.     if (sign)
  1223.         neg_a_bn(n);
  1224.     return r;
  1225.     }
  1226.  
  1227. /**********************************************************************/
  1228. bn_t atan2_bn(bn_t r, bn_t ny, bn_t nx)
  1229.    {
  1230.     copy_bn(bntmpcpy1, ny);
  1231.     copy_bn(bntmpcpy2, nx);
  1232.     unsafe_atan2_bn(r, bntmpcpy1, bntmpcpy2);
  1233.     return r;
  1234.     }
  1235.  
  1236.  
  1237. /**********************************************************************/
  1238. /* Tim's miscellaneous stuff follows                                  */
  1239.  
  1240. /**********************************************************************/
  1241. int is_bn_zero(bn_t n)
  1242.     {
  1243.     return !is_bn_not_zero(n);
  1244.     }
  1245.