home *** CD-ROM | disk | FTP | other *** search
/ Computerworld 1996 March / Computerworld_1996-03_cd.bin / idg_cd3 / grafika / fraktaly / frasr192 / bigflt.c < prev    next >
C/C++ Source or Header  |  1995-03-04  |  50KB  |  1,646 lines
  1. /* bigflt.c - C routines for big floating point numbers */
  2.  
  3. /*
  4. Wesley Loewer's Big Numbers.        (C) 1994, Wesley B. Loewer
  5.  
  6. Big Floating Point Number Format:
  7.  
  8. A big floating point number consists of a signed big integer of length
  9. bflength and of intlength 2 (see bignum.c for details) followed by a
  10. signed 16 bit exponent.  The value of the big floating point number is:
  11.     value = mantissa * 256^exponent
  12.     where the absolute value of the mantissa is in the range 1<=m<256.
  13.  
  14. Notice that this differs from the IEEE format where
  15.     value = mantissa * 2^exponent
  16.     where the absolute value of the mantissa is in the range 1<=m<2.
  17. */
  18.  
  19. #include <stdio.h>
  20. #include <stdlib.h>
  21. #include <string.h>
  22. #include <malloc.h>
  23. #include "prototyp.h"
  24.  
  25. #ifndef _BIGNUM_H
  26. #include "bignum.h"
  27. #endif
  28.  
  29. #ifndef _BIGFLT_H
  30. #include "bigflt.h"
  31. #endif
  32.  
  33. #include "biginit.h"
  34.  
  35. #define LOG10_256 2.4082399653118
  36. #define LOG_256   5.5451774444795
  37.  
  38. #if (_MSC_VER >= 700)
  39. #pragma code_seg ("bigflt1_text")     /* place following in an overlay */
  40. #endif
  41.  
  42. /********************************************************************/
  43. /* bf_hexdump() - for debugging, dumps to stdout                    */
  44.  
  45. void bf_hexdump(bf_t r)
  46.     {
  47.     int i;
  48.  
  49.     for (i=0; i<bflength; i++)
  50.         printf("%02X ", *(r+i));
  51.     printf(" e %04hX ", (S16)access16(r+bflength));
  52.     printf("\n");
  53.     return;
  54.     }
  55.  
  56. /**********************************************************************/
  57. /* strtobf() - converts a string into a bigfloat                       */
  58. /*   r - pointer to a bigfloat                                        */
  59. /*   s - string in the floating point format [+-][dig].[dig]e[+-][dig]*/
  60. /*   note: the string may not be empty or have extra space.           */
  61. /*         It may use scientific notation.                            */
  62. /* USES: bftmp1                                                       */
  63.  
  64. bf_t strtobf(bf_t r, char *s)
  65.     {
  66.     BYTE onesbyte;
  67.     int signflag=0;
  68.     char *l, *d, *e; /* pointer to s, ".", "[eE]" */
  69.     int powerten=0, keeplooping;
  70.  
  71.     clear_bf(r);
  72.  
  73.     if (s[0] == '+')    /* for + sign */
  74.         {
  75.         s++;
  76.         }
  77.     else if (s[0] == '-')    /* for neg sign */
  78.         {
  79.         signflag = 1;
  80.         s++;
  81.         }
  82.  
  83.     d = strchr(s, '.');
  84.     e = strchr(s, 'e');
  85.     if (e == NULL)
  86.         e = strchr(s, 'E');
  87.     if (e != NULL)
  88.         {
  89.         powerten = atoi(e+1);    /* read in the e (x10^) part */
  90.         l = e - 1; /* just before e */
  91.         }
  92.     else
  93.         l = s + strlen(s) - 1;  /* last digit */
  94.  
  95.     if (d != NULL) /* is there a decimal point? */
  96.         {
  97.         while (*l >= '0' && *l <= '9') /* while a digit */
  98.             {
  99.             onesbyte = (BYTE)(*(l--) - '0');
  100.             inttobf(bftmp1,onesbyte);
  101.             unsafe_add_a_bf(r, bftmp1);
  102.             div_a_bf_int(r, 10);
  103.             }
  104.  
  105.         if (*(l--) == '.') /* the digit was found */
  106.             {
  107.             keeplooping = *l >= '0' && *l <= '9' && l>=s;
  108.             while (keeplooping) /* while a digit */
  109.                 {
  110.                 onesbyte = (BYTE)(*(l--) - '0');
  111.                 inttobf(bftmp1,onesbyte);
  112.                 unsafe_add_a_bf(r, bftmp1);
  113.                 keeplooping = *l >= '0' && *l <= '9' && l>=s;
  114.                 if (keeplooping)
  115.                     {
  116.                     div_a_bf_int(r, 10);
  117.                     powerten++;    /* increase the power of ten */
  118.                     }
  119.                 }
  120.             }
  121.         }
  122.     else
  123.         {
  124.         keeplooping = *l >= '0' && *l <= '9' && l>=s;
  125.         while (keeplooping) /* while a digit */
  126.             {
  127.             onesbyte = (BYTE)(*(l--) - '0');
  128.             inttobf(bftmp1,onesbyte);
  129.             unsafe_add_a_bf(r, bftmp1);
  130.             keeplooping = *l >= '0' && *l <= '9' && l>=s;
  131.             if (keeplooping)
  132.                 {
  133.                 div_a_bf_int(r, 10);
  134.                 powerten++;    /* increase the power of ten */
  135.                 }
  136.             }
  137.         }
  138.  
  139.     if (powerten > 0)
  140.         {
  141.         for (; powerten>0; powerten--)
  142.             {
  143.             mult_a_bf_int(r, 10);
  144.             }
  145.         }
  146.     else if (powerten < 0)
  147.         {
  148.         for (; powerten<0; powerten++)
  149.             {
  150.             div_a_bf_int(r, 10);
  151.             }
  152.         }
  153.     if (signflag)
  154.         neg_a_bf(r);
  155.  
  156.     return r;
  157.     }
  158.  
  159. /********************************************************************/
  160. /* strlen_needed() - returns string length needed to hold bigfloat */
  161.  
  162. int strlen_needed_bf()
  163.    {
  164.    int length;
  165.  
  166.    /* first space for integer part */
  167.     length = 1;
  168.     length += decimals;  /* decimal part */
  169.     length += 2;         /* decimal point and sign */
  170.     length += 2;         /* e and sign */
  171.     length += 4;         /* exponent */
  172.     length += 4;         /* null and a little extra for safety */
  173.     return(length);
  174.     }
  175.  
  176. /********************************************************************/
  177. /* bftostr() - converts a bigfloat into a scientific notation string */
  178. /*   s - string, must be large enough to hold the number.           */
  179. /* dec - decimal places, 0 for max                                  */
  180. /*   r - bigfloat                                                   */
  181. /*   will convert to a floating point notation                      */
  182. /*   SIDE-EFFECT: the bigfloat, r, is destroyed.                    */
  183. /*                Copy it first if necessary.                       */
  184. /* USES: bftmp1 - bftmp2                                            */
  185. /********************************************************************/
  186.  
  187. char *unsafe_bftostr(char *s, int dec, bf_t r)
  188.     {
  189.     LDBL value;
  190.     int power;
  191.     int saved;
  192.     bf10_t bf10tmp;
  193.  
  194.     value = bftofloat(r);
  195.     if (value == 0.0)
  196.         {
  197.         strcpy(s, "0.0");
  198.         return s;
  199.         }
  200.  
  201.     saved = save_stack();
  202.     bf10tmp = alloc_stack(dec+4); /* probably ought to be allocated at init */
  203.     copy_bf(bftmp1, r);
  204.     unsafe_bftobf10(bf10tmp, dec, bftmp1);
  205.     power = (S16)access16(bf10tmp+dec+2); /* where the exponent is stored */
  206.     if (power > -4 && power < 6) /* tinker with this */
  207.         bf10tostr_f(s, dec, bf10tmp);
  208.     else
  209.         bf10tostr_e(s, dec, bf10tmp);
  210.     restore_stack(saved);
  211.     return s;
  212.     }
  213.  
  214.  
  215. /********************************************************************/
  216. /* the e version puts it in scientific notation, (like printf's %e) */
  217. char *unsafe_bftostr_e(char *s, int dec, bf_t r)
  218.     {
  219.     LDBL value;
  220.     int saved;
  221.     bf10_t bf10tmp;
  222.  
  223.     value = bftofloat(r);
  224.     if (value == 0.0)
  225.         {
  226.         strcpy(s, "0.0");
  227.         return s;
  228.         }
  229.  
  230.     saved = save_stack();
  231.     bf10tmp = alloc_stack(dec+4); /* probably ought to be allocated at init */
  232.     copy_bf(bftmp1, r);
  233.     unsafe_bftobf10(bf10tmp, dec, bftmp1);
  234.     bf10tostr_e(s, dec, bf10tmp);
  235.     restore_stack(saved);
  236.     return s;
  237.     }
  238.     
  239. /********************************************************************/
  240. /* the f version puts it in decimal notation, (like printf's %f) */
  241. char *unsafe_bftostr_f(char *s, int dec, bf_t r)
  242.     {
  243.     LDBL value;
  244.     int saved;
  245.     bf10_t bf10tmp;
  246.  
  247.     value = bftofloat(r);
  248.     if (value == 0.0)
  249.         {
  250.         strcpy(s, "0.0");
  251.         return s;
  252.         }
  253.  
  254.     saved = save_stack();
  255.     bf10tmp = alloc_stack(dec+4); /* probably ought to be allocated at init */
  256.     copy_bf(bftmp1, r);
  257.     unsafe_bftobf10(bf10tmp, dec, bftmp1);
  258.     bf10tostr_f(s, dec, bf10tmp);
  259.     restore_stack(saved);
  260.     return s;
  261.     }
  262.  
  263. #if (_MSC_VER >= 700)
  264. #pragma code_seg ( )       /* back to normal segment */
  265. #endif
  266.  
  267. /*********************************************************************/
  268. /*  bn = floor(bf)                                                   */
  269. /*  Converts a bigfloat to a bignumber (integer)                     */
  270. /*  bflength must be at least bnlength+2                             */
  271. bn_t bftobn(bn_t n, bf_t f)
  272.     {
  273.     int fexp;
  274.     int movebytes;
  275.     BYTE hibyte;
  276.  
  277.     fexp = (S16)access16(f+bflength);
  278.     if (fexp >= intlength)
  279.         { /* if it's too big, use max value */
  280.         max_bn(n);
  281.         if (is_bf_neg(f))
  282.             neg_a_bn(n);
  283.         return n;
  284.         }
  285.  
  286.     if (-fexp > bnlength - intlength) /* too small, return zero */
  287.         {
  288.         clear_bn(n);
  289.         return n;
  290.         }
  291.  
  292.     /* already checked for over/underflow, this should be ok */
  293.     movebytes = bnlength - intlength + fexp + 1;
  294.     _fmemcpy(n, f+bflength-movebytes-1, movebytes);
  295.     hibyte = *(f+bflength-1);
  296.     _fmemset(n+movebytes, hibyte, bnlength-movebytes); /* sign extends */
  297.     return n;
  298.     }
  299.  
  300. /*********************************************************************/
  301. /*  bf = bn                                                          */
  302. /*  Converts a bignumber (integer) to a bigfloat                     */
  303. /*  bflength must be at least bnlength+2                             */
  304. bf_t bntobf(bf_t f, bn_t n)
  305.     {
  306.     _fmemcpy(f+bflength-bnlength-1, n, bnlength);
  307.     _fmemset(f, 0, bflength - bnlength - 1);
  308.     *(f+bflength-1) = (BYTE)(is_bn_neg(n) ? 0xFF : 0x00); /* sign extend */
  309.     set16(f+bflength, (S16)(intlength - 1)); /* exp */
  310.     norm_bf(f);
  311.     return f;
  312.     }
  313.  
  314. /*********************************************************************/
  315. /*  b = l                                                            */
  316. /*  Converts a long to a bigfloat                                   */
  317. bf_t inttobf(bf_t r, long longval)
  318.     {
  319.     clear_bf(r);
  320.     set32(r+bflength-4, (S32)longval);
  321.     set16(r+bflength, (S16)2);
  322.     norm_bf(r);
  323.     return r;
  324.     }
  325.  
  326. /*********************************************************************/
  327. /*  l = floor(b), floor rounds down                                  */
  328. /*  Converts a bigfloat to a long                                    */
  329. /*  note: a bf value of 2.999... will be return a value of 2, not 3  */
  330. long bftoint(bf_t f)
  331.     {
  332.     int fexp;
  333.     long longval;
  334.  
  335.     fexp = (S16)access16(f+bflength);
  336.     if (fexp > 2)
  337.         {
  338.         longval = 0x7FFFFFFFL;
  339.         if (is_bf_neg(f))
  340.             longval = -longval;
  341.         return longval;
  342.         }
  343.     longval = *(S32 far *)(f+bflength-4);
  344.     longval >>= 8*(2-fexp);
  345.     return longval;
  346.     }
  347.  
  348. /********************************************************************/
  349. /* sign(r)                                                          */
  350. int sign_bf(bf_t n)
  351.     {
  352.     return is_bf_neg(n) ? -1 : is_bf_not_zero(n) ? 1 : 0;
  353.     }
  354.  
  355. /********************************************************************/
  356. /* r = |n|                                                          */
  357. bf_t abs_bf(bf_t r, bf_t n)
  358.     {
  359.     copy_bf(r,n); 
  360.     if (is_bf_neg(r))
  361.        {
  362.        neg_a_bf(r);
  363.        }
  364.     return r;
  365.     }
  366.  
  367. /********************************************************************/
  368. /* r = |r|                                                          */
  369. bf_t abs_a_bf(bf_t r)
  370.     {
  371.     if (is_bf_neg(r))
  372.         neg_a_bf(r);
  373.     return r;
  374.     }
  375.  
  376. /********************************************************************/
  377. /* r = 1/n                                                          */
  378. /* uses bftmp1 - bftmp2 - global temp bigfloats                     */
  379. /*  SIDE-EFFECTS:                                                   */
  380. /*      n ends up as |n|/256^exp    Make copy first if necessary.   */
  381. bf_t unsafe_inv_bf(bf_t r, bf_t n)
  382.     {
  383.     int signflag=0, i;
  384.     int fexp, rexp;
  385.     LDBL f;
  386.     bf_t orig_r, orig_n, orig_bftmp1;
  387.     int  orig_bflength,
  388.          orig_bnlength,
  389.          orig_padding,
  390.          orig_rlength,
  391.          orig_shiftfactor,
  392.          orig_rbflength,
  393.          orig_bfshiftfactor;
  394.  
  395.     /* use Newton's recursive method for zeroing in on 1/n : r=r(2-rn) */
  396.  
  397.     if (is_bf_neg(n))
  398.         { /* will be a lot easier to deal with just positives */
  399.         signflag = 1;
  400.         neg_a_bf(n);
  401.         }
  402.  
  403.     fexp = (S16)access16(n+bflength);
  404.     set16(n+bflength, (S16)0); /* put within LDBL range */
  405.  
  406.     f = bftofloat(n);
  407.     if (f == 0) /* division by zero */
  408.         {
  409.         max_bf(r);
  410.         return r;
  411.         }
  412.     f = 1/f; /* approximate inverse */
  413.  
  414.     /* With Newton's Method, there is no need to calculate all the digits */
  415.     /* every time.  The precision approximately doubles each iteration.   */
  416.     /* Save original values. */
  417.     orig_bflength      = bflength;
  418.     orig_bnlength      = bnlength;
  419.     orig_padding       = padding;
  420.     orig_rlength       = rlength;
  421.     orig_shiftfactor   = shiftfactor;
  422.     orig_rbflength     = rbflength;
  423.     orig_bfshiftfactor = bfshiftfactor;
  424.     orig_r             = r;
  425.     orig_n             = n;
  426.     orig_bftmp1        = bftmp1;
  427.  
  428.     /* calculate new starting values */
  429.     bnlength = intlength + (int)(LDBL_DIG/LOG10_256) + 1; /* round up */
  430.     if (bnlength > orig_bnlength)
  431.         bnlength = orig_bnlength;
  432.     calc_lengths();
  433.  
  434.     /* adjust pointers */
  435.     r = orig_r + orig_bflength - bflength;
  436.     n = orig_n + orig_bflength - bflength;
  437.     bftmp1 = orig_bftmp1 + orig_bflength - bflength;
  438.  
  439.     floattobf(r, f); /* start with approximate inverse */
  440.  
  441.     for (i=0; i<25; i++) /* safety net, this shouldn't ever be needed */
  442.         {
  443.         /* adjust lengths */
  444.         bnlength <<= 1; /* double precision */
  445.         if (bnlength > orig_bnlength)
  446.            bnlength = orig_bnlength;
  447.         calc_lengths();
  448.         r = orig_r + orig_bflength - bflength;
  449.         n = orig_n + orig_bflength - bflength;
  450.         bftmp1 = orig_bftmp1 + orig_bflength - bflength;
  451.  
  452.         unsafe_mult_bf(bftmp1, r, n); /* bftmp1=rn */
  453.         inttobf(bftmp2, 1); /* will be used as 1.0 */
  454.         if (cmp_bf(bftmp1, bftmp2) == 0)
  455.             break;
  456.  
  457.         inttobf(bftmp2, 2); /* will be used as 2.0 */
  458.         unsafe_sub_a_bf(bftmp2, bftmp1); /* bftmp2=2-rn */
  459.         unsafe_mult_bf(bftmp1, r, bftmp2); /* bftmp1=r(2-rn) */
  460.         copy_bf(r, bftmp1); /* r = bftmp1 */
  461.         }
  462.  
  463.     /* restore original values */
  464.     bflength      = orig_bflength;
  465.     bnlength      = orig_bnlength;
  466.     padding       = orig_padding;
  467.     rlength       = orig_rlength;
  468.     shiftfactor   = orig_shiftfactor;
  469.     rbflength     = orig_rbflength;
  470.     bfshiftfactor = orig_bfshiftfactor;
  471.     r             = orig_r;
  472.     n             = orig_n;
  473.     bftmp1        = orig_bftmp1;
  474.  
  475.     if (signflag)
  476.         {
  477.         neg_a_bf(r);
  478.         }
  479.     rexp = (S16)access16(r+bflength);
  480.     rexp -= fexp;
  481.     set16(r+bflength, (S16)rexp); /* adjust result exponent */
  482.     return r;
  483.     }
  484.  
  485. /********************************************************************/
  486. /* r = n1/n2                                                        */
  487. /*      r - result of length bflength                               */
  488. /* uses bftmp1 - bftmp2 - global temp bigfloats                     */
  489. /*  SIDE-EFFECTS:                                                   */
  490. /*      n1, n2 end up as |n1|/256^x, |n2|/256^x                     */
  491. /*      Make copies first if necessary.                             */
  492. bf_t unsafe_div_bf(bf_t r, bf_t n1, bf_t n2)
  493.     {
  494.     int aexp, bexp, rexp;
  495.     LDBL a, b;
  496.  
  497.     /* first, check for valid data */
  498.  
  499.     aexp = (S16)access16(n1+bflength);
  500.     set16(n1+bflength, (S16)0); /* put within LDBL range */
  501.  
  502.     a = bftofloat(n1);
  503.     if (a == 0) /* division into zero */
  504.         {
  505.         clear_bf(r); /* return 0 */
  506.         return r;
  507.         }
  508.  
  509.     bexp = (S16)access16(n2+bflength);
  510.     set16(n2+bflength, (S16)0); /* put within LDBL range */
  511.  
  512.     b = bftofloat(n2);
  513.     if (b == 0) /* division by zero */
  514.         {
  515.         max_bf(r);
  516.         return r;
  517.         }
  518.  
  519.     unsafe_inv_bf(r, n2);
  520.     unsafe_mult_bf(bftmp1, n1, r);
  521.     copy_bf(r, bftmp1); /* r = bftmp1 */
  522.  
  523.     rexp = (S16)access16(r+bflength);
  524.     rexp += aexp - bexp;
  525.     set16(r+bflength, (S16)rexp); /* adjust result exponent */
  526.  
  527.     return r;
  528.     }
  529.  
  530. /********************************************************************/
  531. /* sqrt(r)                                                          */
  532. /* uses bftmp1 - bftmp3 - global temp bigfloats                    */
  533. /*  SIDE-EFFECTS:                                                   */
  534. /*      n ends up as |n|                                            */
  535. bf_t unsafe_sqrt_bf(bf_t r, bf_t n)
  536.     {
  537.     int i, comp, almost_match=0;
  538.     LDBL f;
  539.     bf_t orig_r, orig_n;
  540.     int  orig_bflength,
  541.          orig_bnlength,
  542.          orig_padding,
  543.          orig_rlength,
  544.          orig_shiftfactor,
  545.          orig_rbflength,
  546.          orig_bfshiftfactor;
  547.  
  548. /* use Newton's recursive method for zeroing in on sqrt(n): r=.5(r+n/r) */
  549.  
  550.     if (is_bf_neg(n))
  551.         { /* sqrt of a neg, return 0 */
  552.         clear_bf(r);
  553.         return r;
  554.         }
  555.  
  556.     f = bftofloat(n);
  557.     if (f == 0) /* division by zero will occur */
  558.         {
  559.         clear_bf(r); /* sqrt(0) = 0 */
  560.         return r;
  561.         }
  562.     f = sqrtl(f); /* approximate square root */
  563.     /* no need to check overflow */
  564.  
  565.     /* With Newton's Method, there is no need to calculate all the digits */
  566.     /* every time.  The precision approximately doubles each iteration.   */
  567.     /* Save original values. */
  568.     orig_bflength      = bflength;
  569.     orig_bnlength      = bnlength;
  570.     orig_padding       = padding;
  571.     orig_rlength       = rlength;
  572.     orig_shiftfactor   = shiftfactor;
  573.     orig_rbflength     = rbflength;
  574.     orig_bfshiftfactor = bfshiftfactor;
  575.     orig_r             = r;
  576.     orig_n             = n;
  577.  
  578.     /* calculate new starting values */
  579.     bnlength = intlength + (int)(LDBL_DIG/LOG10_256) + 1; /* round up */
  580.     if (bnlength > orig_bnlength)
  581.         bnlength = orig_bnlength;
  582.     calc_lengths();
  583.  
  584.     /* adjust pointers */
  585.     r = orig_r + orig_bflength - bflength;
  586.     n = orig_n + orig_bflength - bflength;
  587.  
  588.     floattobf(r, f); /* start with approximate sqrt */
  589.  
  590.     for (i=0; i<25; i++) /* safety net, this shouldn't ever be needed */
  591.         {
  592.         /* adjust lengths */
  593.         bnlength <<= 1; /* double precision */
  594.         if (bnlength > orig_bnlength)
  595.            bnlength = orig_bnlength;
  596.         calc_lengths();
  597.         r = orig_r + orig_bflength - bflength;
  598.         n = orig_n + orig_bflength - bflength;
  599.  
  600.         unsafe_div_bf(bftmp3, n, r);
  601.         unsafe_add_a_bf(r, bftmp3);
  602.         half_a_bf(r);
  603.         if ((comp=abs(cmp_bf(r, bftmp3))) <= 2 ) /* if match or almost match */
  604.             {
  605.             if (comp == 0)  /* perfect match */
  606.                 break;
  607.             if (almost_match == 2) /* did they almost match before? */
  608.                 { /* yes, this must be the third time */
  609.                 break;
  610.                 }
  611.             else /* this is the first or second time they almost matched */
  612.                 almost_match++;
  613.             }
  614.         }
  615.  
  616.     /* restore original values */
  617.     bflength      = orig_bflength;
  618.     bnlength      = orig_bnlength;
  619.     padding       = orig_padding;
  620.     rlength       = orig_rlength;
  621.     shiftfactor   = orig_shiftfactor;
  622.     rbflength     = orig_rbflength;
  623.     bfshiftfactor = orig_bfshiftfactor;
  624.     r             = orig_r;
  625.     n             = orig_n;
  626.  
  627.     return r;
  628.     }
  629.  
  630. /********************************************************************/
  631. /* exp(r)                                                           */
  632. /* uses bftmp1, bftmp2, bftmp3 - global temp bigfloats             */
  633. bf_t exp_bf(bf_t r, bf_t n)
  634.     {
  635.     U16 fact=1;
  636.     S16 far * testexp, far * rexp;
  637.  
  638.     testexp = (S16 far *)(bftmp2+bflength);
  639.     rexp = (S16 far *)(r+bflength);
  640.  
  641.     if (is_bf_zero(n))
  642.         {
  643.         inttobf(r, 1);
  644.         return r;
  645.         }
  646.  
  647. /* use Taylor Series (very slow convergence) */
  648.     inttobf(r, 1); /* start with r=1.0 */
  649.     copy_bf(bftmp2, r);
  650.     for (;;)
  651.         {
  652.         copy_bf(bftmp1, n);
  653.         unsafe_mult_bf(bftmp3, bftmp2, bftmp1);
  654.         unsafe_div_bf_int(bftmp2, bftmp3, fact);
  655.         if (*testexp < *rexp-(bflength-2))
  656.             break; /* too small to register */
  657.         unsafe_add_a_bf(r, bftmp2);
  658.         fact++;
  659.         }
  660.  
  661.     return r;
  662.     }
  663.  
  664. /********************************************************************/
  665. /* ln(r)                                                            */
  666. /* uses bftmp1 - bftmp6 - global temp bigfloats                     */
  667. /*  SIDE-EFFECTS:                                                   */
  668. /*      n ends up as |n|                                            */
  669. bf_t unsafe_ln_bf(bf_t r, bf_t n)
  670.     {
  671.     int i, comp, almost_match=0;
  672.     LDBL f;
  673.     bf_t orig_r, orig_n, orig_bftmp5;
  674.     int  orig_bflength,
  675.          orig_bnlength,
  676.          orig_padding,
  677.          orig_rlength,
  678.          orig_shiftfactor,
  679.          orig_rbflength,
  680.          orig_bfshiftfactor;
  681.  
  682. /* use Newton's recursive method for zeroing in on ln(n): r=r+n*exp(-r)-1 */
  683.  
  684.     if (is_bf_neg(n) || is_bf_zero(n))
  685.         { /* error, return largest neg value */
  686.         max_bf(r);
  687.         neg_a_bf(r);
  688.         return r;
  689.         }
  690.  
  691.     f = bftofloat(n);
  692.     f = logl(f); /* approximate ln(x) */
  693.     /* no need to check overflow */
  694.     /* appears to be ok, do ln */
  695.  
  696.     /* With Newton's Method, there is no need to calculate all the digits */
  697.     /* every time.  The precision approximately doubles each iteration.   */
  698.     /* Save original values. */
  699.     orig_bflength      = bflength;
  700.     orig_bnlength      = bnlength;
  701.     orig_padding       = padding;
  702.     orig_rlength       = rlength;
  703.     orig_shiftfactor   = shiftfactor;
  704.     orig_rbflength     = rbflength;
  705.     orig_bfshiftfactor = bfshiftfactor;
  706.     orig_r             = r;
  707.     orig_n             = n;
  708.     orig_bftmp5        = bftmp5;
  709.  
  710.     /* calculate new starting values */
  711.     bnlength = intlength + (int)(LDBL_DIG/LOG10_256) + 1; /* round up */
  712.     if (bnlength > orig_bnlength)
  713.         bnlength = orig_bnlength;
  714.     calc_lengths();
  715.  
  716.     /* adjust pointers */
  717.     r = orig_r + orig_bflength - bflength;
  718.     n = orig_n + orig_bflength - bflength;
  719.     bftmp5 = orig_bftmp5 + orig_bflength - bflength;
  720.  
  721.     floattobf(r, f); /* start with approximate ln */
  722.     neg_a_bf(r); /* -r */
  723.     copy_bf(bftmp5, r); /* -r */
  724.  
  725.     for (i=0; i<25; i++) /* safety net, this shouldn't ever be needed */
  726.         {
  727.         /* adjust lengths */
  728.         bnlength <<= 1; /* double precision */
  729.         if (bnlength > orig_bnlength)
  730.            bnlength = orig_bnlength;
  731.         calc_lengths();
  732.         r = orig_r + orig_bflength - bflength;
  733.         n = orig_n + orig_bflength - bflength;
  734.         bftmp5 = orig_bftmp5 + orig_bflength - bflength;
  735.  
  736.         exp_bf(bftmp6, r);     /* exp(-r) */
  737.         unsafe_mult_bf(bftmp2, bftmp6, n);  /* n*exp(-r) */
  738.         inttobf(bftmp4, 1);
  739.         unsafe_sub_a_bf(bftmp2, bftmp4);   /* n*exp(-r) - 1 */
  740.         unsafe_sub_a_bf(r, bftmp2);        /* -r - (n*exp(-r) - 1) */
  741.         if ((comp=abs(cmp_bf(r, bftmp5))) <= 2 ) /* if match or almost match */
  742.             {
  743.             if (comp == 0)  /* perfect match */
  744.                 break;
  745.             if (almost_match == 2) /* did they almost match twice before? */
  746.                 { /* yes, this must be the third time */
  747.                 break;
  748.                 }
  749.             else /* this is the first or second time they almost matched */
  750.                 almost_match++;
  751.             }
  752.         copy_bf(bftmp5, r); /* -r */
  753.         }
  754.  
  755.     /* restore original values */
  756.     bflength      = orig_bflength;
  757.     bnlength      = orig_bnlength;
  758.     padding       = orig_padding;
  759.     rlength       = orig_rlength;
  760.     shiftfactor   = orig_shiftfactor;
  761.     rbflength     = orig_rbflength;
  762.     bfshiftfactor = orig_bfshiftfactor;
  763.     r             = orig_r;
  764.     n             = orig_n;
  765.     bftmp5        = orig_bftmp5;
  766.  
  767.     neg_a_bf(r); /* -(-r) */
  768.     return r;
  769.     }
  770.  
  771. /********************************************************************/
  772. /* sincos_bf(r)                                                     */
  773. /* uses bftmp1 - bftmp2 - global temp bigfloats                     */
  774. /*  SIDE-EFFECTS:                                                   */
  775. /*      n ends up as |n| mod (pi/4)                                 */
  776. bf_t unsafe_sincos_bf(bf_t s, bf_t c, bf_t n)
  777.     {
  778.     U16 fact=2;
  779.     int k=0, i, halves;
  780.     int signcos=0, signsin=0, switch_sincos=0;
  781.     int sin_done=0, cos_done=0;
  782.     S16 far * testexp, far * cexp, far * sexp;
  783.  
  784.     testexp = (S16 far *)(bftmp1+bflength);
  785.     cexp = (S16 far *)(c+bflength);
  786.     sexp = (S16 far *)(s+bflength);
  787.  
  788. #if !CALCULATING_BIG_PI
  789.     /* assure range 0 <= x < pi/4 */
  790.  
  791.     if (is_bf_zero(n))
  792.         {
  793.         clear_bf(s);    /* sin(0) = 0 */
  794.         inttobf(c, 1);  /* cos(0) = 1 */
  795.         return s;
  796.         }
  797.  
  798.     if (is_bf_neg(n))
  799.         {
  800.         signsin = !signsin; /* sin(-x) = -sin(x), odd; cos(-x) = cos(x), even */
  801.         neg_a_bf(n);
  802.         }
  803.     /* n >= 0 */
  804.  
  805.     double_bf(bftmp1, bf_pi); /* 2*pi */
  806.     /* this could be done with remainders, but it would probably be slower */
  807.     while (cmp_bf(n, bftmp1) >= 0) /* while n >= 2*pi */
  808.         {
  809.         copy_bf(bftmp2, bftmp1);
  810.         unsafe_sub_a_bf(n, bftmp2);
  811.         }
  812.     /* 0 <= n < 2*pi */
  813.  
  814.     copy_bf(bftmp1, bf_pi); /* pi */
  815.     if (cmp_bf(n, bftmp1) >= 0) /* if n >= pi */
  816.         {
  817.         unsafe_sub_a_bf(n, bftmp1);
  818.         signsin = !signsin;
  819.         signcos = !signcos;
  820.         }
  821.     /* 0 <= n < pi */
  822.  
  823.     half_bf(bftmp1, bf_pi); /* pi/2 */
  824.     if (cmp_bf(n, bftmp1) > 0) /* if n > pi/2 */
  825.         {
  826.         copy_bf(bftmp2, bf_pi);
  827.         unsafe_sub_bf(n, bftmp2, n);
  828.         signcos = !signcos;
  829.         }
  830.     /* 0 <= n < pi/2 */
  831.  
  832.     half_bf(bftmp1, bf_pi); /* pi/2 */
  833.     half_a_bf(bftmp1);      /* pi/4 */
  834.     if (cmp_bf(n, bftmp1) > 0) /* if n > pi/4 */
  835.         {
  836.         copy_bf(bftmp2, n);
  837.         half_bf(bftmp1, bf_pi); /* pi/2 */
  838.         unsafe_sub_bf(n, bftmp1, bftmp2);  /* pi/2 - n */
  839.         switch_sincos = !switch_sincos;
  840.         }
  841.     /* 0 <= n < pi/4 */
  842.  
  843.     /* this looks redundant, but n could now be zero when it wasn't before */
  844.     if (is_bf_zero(n))
  845.         {
  846.         clear_bf(s);    /* sin(0) = 0 */
  847.         inttobf(c, 1);  /* cos(0) = 1 */
  848.         return s;
  849.         }
  850.  
  851.  
  852. /* at this point, the double angle trig identities could be used as many  */
  853. /* times as desired to reduce the range to pi/8, pi/16, etc...  Each time */
  854. /* the range is cut in half, the number of iterations required is reduced */
  855. /* by "quite a bit."  It's just a matter of testing to see what gives the */
  856. /* optimal results.                                                       */
  857.     /* halves = bflength / 10; */ /* this is experimental */
  858.     halves = 1;
  859.     for (i = 0; i < halves; i++)
  860.         half_a_bf(n);
  861. #endif
  862.  
  863. /* use Taylor Series (very slow convergence) */
  864.     copy_bf(s, n); /* start with s=n */
  865.     inttobf(c, 1); /* start with c=1 */
  866.     copy_bf(bftmp1, n); /* the current x^n/n! */
  867.     do
  868.         {
  869.         /* even terms for cosine */
  870.         copy_bf(bftmp2, bftmp1);
  871.         unsafe_mult_bf(bftmp1, bftmp2, n);
  872.         div_a_bf_int(bftmp1, fact++);
  873.         if (!cos_done)
  874.             {
  875.             cos_done = (*testexp < *cexp-(bflength-2)); /* too small to register */
  876.             if (!cos_done)
  877.                 {
  878.                 if (k) /* alternate between adding and subtracting */
  879.                     unsafe_add_a_bf(c, bftmp1);
  880.                 else
  881.                     unsafe_sub_a_bf(c, bftmp1);
  882.                 }
  883.             }
  884.  
  885.         /* odd terms for sine */
  886.         copy_bf(bftmp2, bftmp1);
  887.         unsafe_mult_bf(bftmp1, bftmp2, n);
  888.         div_a_bf_int(bftmp1, fact++);
  889.         if (!sin_done)
  890.             {
  891.             sin_done = (*testexp < *sexp-(bflength-2)); /* too small to register */
  892.             if (!sin_done)
  893.                 {
  894.                 if (k) /* alternate between adding and subtracting */
  895.                     unsafe_add_a_bf(s, bftmp1);
  896.                 else
  897.                     unsafe_sub_a_bf(s, bftmp1);
  898.                 }
  899.             }
  900.         k = !k; /* toggle */
  901. #if CALCULATING_BIG_PI
  902.         printf("."); /* lets you know it's doing something */
  903. #endif
  904.  
  905.         } while (!cos_done || !sin_done);
  906.  
  907. #if !CALCULATING_BIG_PI
  908.         /* now need to undo what was done by cutting angles in half */
  909.         for (i = 0; i < halves; i++)
  910.             {
  911.             unsafe_mult_bf(bftmp2, s, c); /* no need for safe mult */
  912.             double_bf(s, bftmp2); /* sin(2x) = 2*sin(x)*cos(x) */
  913.             unsafe_square_bf(bftmp2,c);
  914.             double_a_bf(bftmp2);
  915.             inttobf(bftmp1, 1);
  916.             unsafe_sub_bf(c, bftmp2, bftmp1); /* cos(2x) = 2*cos(x)*cos(x) - 1 */
  917.             }
  918.  
  919.         if (switch_sincos)
  920.         {
  921.             copy_bf(bftmp1, s);
  922.             copy_bf(s, c);
  923.             copy_bf(c, bftmp1);
  924.         }
  925.     if (signsin)
  926.             neg_a_bf(s);
  927.     if (signcos)
  928.             neg_a_bf(c);
  929. #endif
  930.  
  931.     return s; /* return sine I guess */
  932.     }
  933.  
  934. /********************************************************************/
  935. /* atan(r)                                                          */
  936. /* uses bftmp1 - bftmp5 - global temp bigfloats                    */
  937. /*  SIDE-EFFECTS:                                                   */
  938. /*      n ends up as |n| or 1/|n|                                   */
  939. bf_t unsafe_atan_bf(bf_t r, bf_t n)
  940.     {
  941.     int i, comp, almost_match=0, signflag=0;
  942.     LDBL f;
  943.     bf_t orig_r, orig_n, orig_bf_pi, orig_bftmp3;
  944.     int  orig_bflength,
  945.          orig_bnlength,
  946.          orig_padding,
  947.          orig_rlength,
  948.          orig_shiftfactor,
  949.          orig_rbflength,
  950.          orig_bfshiftfactor;
  951.     int large_arg;
  952.  
  953. /* use Newton's recursive method for zeroing in on atan(n): r=r-cos(r)(sin(r)-n*cos(r)) */
  954.  
  955.     if (is_bf_neg(n))
  956.         {
  957.         signflag = 1;
  958.         neg_a_bf(n);
  959.         }
  960.  
  961. /* If n is very large, atanl() won't give enough decimal places to be a */
  962. /* good enough initial guess for Newton's Method.  If it is larger than */
  963. /* say, 1, atan(n) = pi/2 - acot(n) = pi/2 - atan(1/n).                 */
  964.  
  965.     f = bftofloat(n);
  966.     large_arg = f > 1.0;
  967.     if (large_arg)
  968.         {
  969.         unsafe_inv_bf(bftmp3, n);
  970.         copy_bf(n, bftmp3);
  971.         f = bftofloat(n);
  972.         }
  973.  
  974.     /* With Newton's Method, there is no need to calculate all the digits */
  975.     /* every time.  The precision approximately doubles each iteration.   */
  976.     /* Save original values. */
  977.     orig_bflength      = bflength;
  978.     orig_bnlength      = bnlength;
  979.     orig_padding       = padding;
  980.     orig_rlength       = rlength;
  981.     orig_shiftfactor   = shiftfactor;
  982.     orig_rbflength     = rbflength;
  983.     orig_bfshiftfactor = bfshiftfactor;
  984.     orig_bf_pi         = bf_pi;
  985.     orig_r             = r;
  986.     orig_n             = n;
  987.     orig_bftmp3        = bftmp3;
  988.  
  989.     /* calculate new starting values */
  990.     bnlength = intlength + (int)(LDBL_DIG/LOG10_256) + 1; /* round up */
  991.     if (bnlength > orig_bnlength)
  992.         bnlength = orig_bnlength;
  993.     calc_lengths();
  994.  
  995.     /* adjust pointers */
  996.     r = orig_r + orig_bflength - bflength;
  997.     n = orig_n + orig_bflength - bflength;
  998.     bf_pi = orig_bf_pi + orig_bflength - bflength;
  999.     bftmp3 = orig_bftmp3 + orig_bflength - bflength;
  1000.  
  1001.     f = atanl(f); /* approximate arctangent */
  1002.     /* no need to check overflow */
  1003.  
  1004.     floattobf(r, f); /* start with approximate atan */
  1005.     copy_bf(bftmp3, r);
  1006.  
  1007.  
  1008. #if CALCULATING_BIG_PI
  1009.     setvideomode(3,0,0,0); /* put in text mode */
  1010. #endif
  1011.  
  1012.     for (i=0; i<25; i++) /* safety net, this shouldn't ever be needed */
  1013.         {
  1014.         /* adjust lengths */
  1015.         bnlength <<= 1; /* double precision */
  1016.         if (bnlength > orig_bnlength)
  1017.            bnlength = orig_bnlength;
  1018.         calc_lengths();
  1019.         r = orig_r + orig_bflength - bflength;
  1020.         n = orig_n + orig_bflength - bflength;
  1021.         bf_pi = orig_bf_pi + orig_bflength - bflength;
  1022.         bftmp3 = orig_bftmp3 + orig_bflength - bflength;
  1023.  
  1024. #if CALCULATING_BIG_PI
  1025.         printf("\natan() loop #%i, bflength=%i\nsincos() loops\n", i, bflength);
  1026. #endif
  1027.         unsafe_sincos_bf(bftmp4, bftmp5, bftmp3);   /* sin(r), cos(r) */
  1028.         copy_bf(bftmp3, r); /* restore bftmp3 from sincos_bf() */
  1029.         copy_bf(bftmp1, bftmp5);
  1030.         unsafe_mult_bf(bftmp2, n, bftmp1);     /* n*cos(r) */
  1031.         unsafe_sub_a_bf(bftmp4, bftmp2); /* sin(r) - n*cos(r) */
  1032.         unsafe_mult_bf(bftmp1, bftmp5, bftmp4); /* cos(r) * (sin(r) - n*cos(r)) */
  1033.         copy_bf(bftmp3, r);
  1034.         unsafe_sub_a_bf(r, bftmp1); /* r - cos(r) * (sin(r) - n*cos(r)) */
  1035.         if (bflength == orig_bflength && (comp=abs(cmp_bf(r, bftmp3))) <= 2 ) /* if match or almost match */
  1036.             {
  1037. #if CALCULATING_BIG_PI
  1038.             printf("\natan() loop comp=%i\n", comp);
  1039. #endif
  1040.             if (comp == 0)  /* perfect match */
  1041.                 break;
  1042.             if (almost_match == 2) /* did they almost match before? */
  1043.                 { /* yes, this must be the third time */
  1044.                 break;
  1045.                 }
  1046.             else /* this is the first or second time they almost matched */
  1047.                 almost_match++;
  1048.             }
  1049. #if CALCULATING_BIG_PI
  1050.         if (bflength == orig_bflength && comp > 2)
  1051.             printf("\natan() loop comp=%i\n", comp);
  1052. #endif
  1053.  
  1054.         copy_bf(bftmp3, r); /* make a copy for later comparison */
  1055.         }
  1056.  
  1057.     /* restore original values */
  1058.     bflength      = orig_bflength;
  1059.     bnlength      = orig_bnlength;
  1060.     padding       = orig_padding;
  1061.     rlength       = orig_rlength;
  1062.     shiftfactor   = orig_shiftfactor;
  1063.     rbflength     = orig_rbflength;
  1064.     bfshiftfactor = orig_bfshiftfactor;
  1065.     bf_pi         = orig_bf_pi;
  1066.     r             = orig_r;
  1067.     n             = orig_n;
  1068.     bftmp3        = orig_bftmp3;
  1069.  
  1070.     if (large_arg)
  1071.         {
  1072.         half_bf(bftmp3, bf_pi);  /* pi/2 */
  1073.         sub_a_bf(bftmp3, r);     /* pi/2 - atan(1/n) */
  1074.         copy_bf(r, bftmp3);
  1075.         }
  1076.  
  1077.     if (signflag)
  1078.         neg_a_bf(r);
  1079.     return r;
  1080.     }
  1081.  
  1082. /********************************************************************/
  1083. /* atan2(r,ny,nx)                                                     */
  1084. /* uses bftmp1 - bftmp6 - global temp bigfloats                     */
  1085. bf_t unsafe_atan2_bf(bf_t r, bf_t ny, bf_t nx)
  1086.    {
  1087.    int signx, signy;
  1088.  
  1089.    signx = sign_bf(nx);
  1090.    signy = sign_bf(ny);
  1091.  
  1092.    if (signy == 0)
  1093.       {
  1094.       if (signx < 0)
  1095.          copy_bf(r, bf_pi); /* negative x axis, 180 deg */
  1096.       else    /* signx >= 0    positive x axis, 0 */
  1097.          clear_bf(r);
  1098.       return(r);
  1099.       }
  1100.    if (signx == 0)
  1101.       {
  1102.       copy_bf(r, bf_pi); /* y axis */
  1103.       half_a_bf(r);      /* +90 deg */
  1104.       if (signy < 0)
  1105.          neg_a_bf(r);    /* -90 deg */
  1106.       return(r);
  1107.       }
  1108.  
  1109.    if (signy < 0)
  1110.       neg_a_bf(ny);
  1111.    if (signx < 0)
  1112.       neg_a_bf(nx);
  1113.    unsafe_div_bf(bftmp6,ny,nx);
  1114.    unsafe_atan_bf(r, bftmp6);
  1115.    if (signx < 0)
  1116.       sub_bf(r,bf_pi,r);
  1117.    if(signy < 0)
  1118.       neg_a_bf(r);
  1119.    return(r);
  1120.    }
  1121.  
  1122. /**********************************************************************/
  1123. /* The rest of the functions are "safe" versions of the routines that */
  1124. /* have side effects which alter the parameters.                      */
  1125. /* Most bf routines change values of parameters, not just the sign.   */
  1126. /**********************************************************************/
  1127.  
  1128. /**********************************************************************/
  1129. bf_t add_bf(bf_t r, bf_t n1, bf_t n2)
  1130.     {
  1131.     copy_bf(bftmpcpy1, n1);
  1132.     copy_bf(bftmpcpy2, n2);
  1133.     unsafe_add_bf(r, bftmpcpy1, bftmpcpy2);
  1134.     return r;
  1135.     }
  1136.  
  1137. /**********************************************************************/
  1138. bf_t add_a_bf(bf_t r, bf_t n)
  1139.     {
  1140.     copy_bf(bftmpcpy1, n);
  1141.     unsafe_add_a_bf(r, bftmpcpy1);
  1142.     return r;
  1143.     }
  1144.  
  1145. /**********************************************************************/
  1146. bf_t sub_bf(bf_t r, bf_t n1, bf_t n2)
  1147.     {
  1148.     copy_bf(bftmpcpy1, n1);
  1149.     copy_bf(bftmpcpy2, n2);
  1150.     unsafe_sub_bf(r, bftmpcpy1, bftmpcpy2);
  1151.     return r;
  1152.     }
  1153.  
  1154. /**********************************************************************/
  1155. bf_t sub_a_bf(bf_t r, bf_t n)
  1156.     {
  1157.     copy_bf(bftmpcpy1, n);
  1158.     unsafe_sub_a_bf(r, bftmpcpy1);
  1159.     return r;
  1160.     }
  1161.  
  1162. /**********************************************************************/
  1163. /* mult and div only change sign */
  1164. bf_t full_mult_bf(bf_t r, bf_t n1, bf_t n2)
  1165.     {
  1166.     copy_bf(bftmpcpy1, n1);
  1167.     copy_bf(bftmpcpy2, n2);
  1168.     unsafe_full_mult_bf(r, bftmpcpy1, bftmpcpy2);
  1169.     return r;
  1170.     }
  1171.  
  1172. /**********************************************************************/
  1173. bf_t mult_bf(bf_t r, bf_t n1, bf_t n2)
  1174.     {
  1175.     copy_bf(bftmpcpy1, n1);
  1176.     copy_bf(bftmpcpy2, n2);
  1177.     unsafe_mult_bf(r, bftmpcpy1, bftmpcpy2);
  1178.     return r;
  1179.     }
  1180.  
  1181. /**********************************************************************/
  1182. bf_t full_square_bf(bf_t r, bf_t n)
  1183.     {
  1184.     copy_bf(bftmpcpy1, n);
  1185.     unsafe_full_square_bf(r, bftmpcpy1);
  1186.     return r;
  1187.     }
  1188.  
  1189. /**********************************************************************/
  1190. bf_t square_bf(bf_t r, bf_t n)
  1191.     {
  1192.     copy_bf(bftmpcpy1, n);
  1193.     unsafe_square_bf(r, bftmpcpy1);
  1194.     return r;
  1195.     }
  1196.  
  1197. /**********************************************************************/
  1198. bf_t mult_bf_int(bf_t r, bf_t n, U16 u)
  1199.     {
  1200.     copy_bf(bftmpcpy1, n);
  1201.     unsafe_mult_bf_int(r, bftmpcpy1, u);
  1202.     return r;
  1203.     }
  1204.  
  1205. /**********************************************************************/
  1206. bf_t div_bf_int(bf_t r, bf_t n,  U16 u)
  1207.     {
  1208.     copy_bf(bftmpcpy1, n);
  1209.     unsafe_div_bf_int(r, bftmpcpy1, u);
  1210.     return r;
  1211.     }
  1212.  
  1213. #if (_MSC_VER >= 700)
  1214. #pragma code_seg ("bigflt1_text")     /* place following in an overlay */
  1215. #endif
  1216.  
  1217. /**********************************************************************/
  1218. char *bftostr(char *s, int dec, bf_t r)
  1219.     {
  1220.     copy_bf(bftmpcpy1, r);
  1221.     unsafe_bftostr(s, dec, bftmpcpy1);
  1222.     return s;
  1223.     }
  1224.  
  1225. /**********************************************************************/
  1226. char *bftostr_e(char *s, int dec, bf_t r)
  1227.     {
  1228.     copy_bf(bftmpcpy1, r);
  1229.     unsafe_bftostr_e(s, dec, bftmpcpy1);
  1230.     return s;
  1231.     }
  1232.  
  1233. /**********************************************************************/
  1234. char *bftostr_f(char *s, int dec, bf_t r)
  1235.     {
  1236.     copy_bf(bftmpcpy1, r);
  1237.     unsafe_bftostr_f(s, dec, bftmpcpy1);
  1238.     return s;
  1239.     }
  1240.  
  1241. #if (_MSC_VER >= 700)
  1242. #pragma code_seg ( )       /* back to normal segment */
  1243. #endif
  1244.  
  1245. /**********************************************************************/
  1246. bf_t inv_bf(bf_t r, bf_t n)
  1247.     {
  1248.     copy_bf(bftmpcpy1, n);
  1249.     unsafe_inv_bf(r, bftmpcpy1);
  1250.     return r;
  1251.     }
  1252.  
  1253. /**********************************************************************/
  1254. bf_t div_bf(bf_t r, bf_t n1, bf_t n2)
  1255.     {
  1256.     copy_bf(bftmpcpy1, n1);
  1257.     copy_bf(bftmpcpy2, n2);
  1258.     unsafe_div_bf(r, bftmpcpy1, bftmpcpy2);
  1259.     return r;
  1260.     }
  1261.  
  1262. /**********************************************************************/
  1263. bf_t sqrt_bf(bf_t r, bf_t n)
  1264.     {
  1265.     copy_bf(bftmpcpy1, n);
  1266.     unsafe_sqrt_bf(r, bftmpcpy1);
  1267.     return r;
  1268.     }
  1269.  
  1270. /**********************************************************************/
  1271. bf_t ln_bf(bf_t r, bf_t n)
  1272.     {
  1273.     copy_bf(bftmpcpy1, n);
  1274.     unsafe_ln_bf(r, bftmpcpy1);
  1275.     return r;
  1276.     }
  1277.  
  1278. /**********************************************************************/
  1279. bf_t sincos_bf(bf_t s, bf_t c, bf_t n)
  1280.     {
  1281.     copy_bf(bftmpcpy1, n);
  1282.     return unsafe_sincos_bf(s, c, bftmpcpy1);
  1283.     }
  1284.  
  1285. /**********************************************************************/
  1286. bf_t atan_bf(bf_t r, bf_t n)
  1287.     {
  1288.     copy_bf(bftmpcpy1, n);
  1289.     unsafe_atan_bf(r, bftmpcpy1);
  1290.     return r;
  1291.     }
  1292.  
  1293. /**********************************************************************/
  1294. bf_t atan2_bf(bf_t r, bf_t ny, bf_t nx)
  1295.    {
  1296.     copy_bf(bftmpcpy1, ny);
  1297.     copy_bf(bftmpcpy2, nx);
  1298.     unsafe_atan2_bf(r, bftmpcpy1, bftmpcpy2);
  1299.     return r;
  1300.     }
  1301.  
  1302. /**********************************************************************/
  1303. int is_bf_zero(bf_t n)
  1304.     {
  1305.     return !is_bf_not_zero(n);
  1306.     }
  1307.  
  1308. /************************************************************************/
  1309. /* convert_bf  -- convert bigfloat numbers from old to new lengths      */
  1310. int convert_bf(bf_t new, bf_t old, int newbflength, int oldbflength)
  1311.    {
  1312.    int savebflength;
  1313.  
  1314.    /* save lengths so not dependent on external environment */
  1315.    savebflength  = bflength;
  1316.    bflength      = newbflength;
  1317.    clear_bf(new);
  1318.    bflength      = savebflength;
  1319.  
  1320.    if(newbflength > oldbflength)
  1321.       _fmemcpy(new+newbflength-oldbflength,old,oldbflength+2);
  1322.    else
  1323.       _fmemcpy(new,old+oldbflength-newbflength,newbflength+2);
  1324.    return(0);
  1325.    }
  1326.  
  1327. /* big10flt.c - C routines for base 10 big floating point numbers */
  1328.  
  1329. /**********************************************************
  1330. (Just when you thought it was safe to go back in the water.)
  1331. Just when you thought you seen every type of format possible,
  1332. 16 bit integer, 32 bit integer, double, long double, mpmath,
  1333. bn_t, bf_t, I now give you bf10_t (big float base 10)!
  1334.  
  1335. Why, because this is the only way (I can think of) to properly do a
  1336. bftostr() without rounding errors.  With out this, then
  1337.    -1.9999999999( > LDBL_DIG of 9's)9999999123456789...
  1338. will round to -2.0.  The good news is that we only need to do two
  1339. mathematical operations: multiplication and division by integers
  1340.  
  1341. bf10_t format: (notice the position of the MSB and LSB)
  1342.  
  1343. MSB                                         LSB
  1344.   _  _  _  _  _  _  _  _  _  _  _  _ _ _ _ _
  1345. n <><------------- dec --------------><> <->
  1346.   1 byte pad            1 byte rounding   2 byte exponent.
  1347.  
  1348.   total length = dec + 4
  1349.  
  1350. ***********************************************************/
  1351.  
  1352. /**********************************************************************/
  1353. /* unsafe_bftobf10() - converts a bigfloat into a bigfloat10          */
  1354. /*   n - pointer to a bigfloat                                        */
  1355. /*   r - result array of BYTE big enough to hold the bf10_t number    */
  1356. /* dec - number of decimals, not including the one extra for rounding */
  1357. /*  SIDE-EFFECTS: n is changed to |n|.  Make copy of n if necessary.  */
  1358.  
  1359. bf10_t unsafe_bftobf10(bf10_t r, int dec, bf_t n)
  1360.    {
  1361.    int d;
  1362.    int power256;
  1363.    int p;
  1364.    int bnl;
  1365.    bf_t onesbyte;
  1366.    bf10_t power10;
  1367.  
  1368.    if (is_bf_zero(n))
  1369.       { /* in scientific notation, the leading digit can't be zero */
  1370.       r[1] = (BYTE)0; /* unless the number is zero */
  1371.       return r;
  1372.       }
  1373.  
  1374.    onesbyte = n + bflength - 1;           /* really it's n+bflength-2 */
  1375.    power256 = (S16)access16(n + bflength) + 1; /* so adjust power256 by 1  */
  1376.  
  1377.    if (dec == 0)
  1378.        dec = decimals;
  1379.    dec++;  /* one extra byte for rounding */
  1380.    power10 = r + dec + 1;
  1381.  
  1382.    if (is_bf_neg(n))
  1383.       {
  1384.       neg_a_bf(n);
  1385.       r[0] = 1; /* sign flag */
  1386.       }
  1387.    else
  1388.       r[0] = 0;
  1389.  
  1390.    p = -1;  /* multiply by 10 right away */
  1391.    bnl = bnlength;
  1392.    bnlength = bflength;
  1393.    for (d=1; d<=dec; d++)
  1394.       {
  1395.       /* pretend it's a bn_t instead of a bf_t */
  1396.       /* this leaves n un-normalized, which is what we want here  */
  1397.       mult_a_bn_int(n, 10);
  1398.  
  1399.       r[d] = *onesbyte;
  1400.       if (d == 1 && r[d] == 0)
  1401.          {
  1402.          d = 0; /* back up a digit */
  1403.          p--; /* and decrease by a factor of 10 */
  1404.          }
  1405.       *onesbyte = 0;
  1406.       }
  1407.    bnlength = bnl;
  1408.    set16(power10, p); /* save power of ten */
  1409.  
  1410.    /* the digits are all read in, now scale it by 256^power256 */
  1411.    if (power256 > 0)
  1412.       for (d=0; d<power256; d++)
  1413.          mult_a_bf10_int(r, dec, 256);
  1414.  
  1415.    else if (power256 < 0)
  1416.       for (d=0; d>power256; d--)
  1417.          div_a_bf10_int(r, dec, 256);
  1418.  
  1419.    /* else power256 is zero, don't do anything */
  1420.  
  1421.    /* round the last digit */
  1422.    if (r[dec] >= 5)
  1423.       {
  1424.       d = dec-1;
  1425.       while (d > 0) /* stop before you get to the sign flag */
  1426.          {
  1427.          r[d]++;  /* round up */
  1428.          if (r[d] < 10)
  1429.             {
  1430.             d = -1; /* flag for below */
  1431.             break; /* finished rounding */
  1432.             }
  1433.          r[d] = 0;
  1434.          d--;
  1435.          }
  1436.       if (d == 0) /* rounding went back to the first digit and it overflowed */
  1437.          {
  1438.          r[1] = 0;
  1439.          _fmemmove(r+2, r+1, dec-1);
  1440.          r[1] = 1;
  1441.          p = (S16)access16(power10);
  1442.          set16(power10, p+1);
  1443.          }
  1444.       }
  1445.    r[dec] = 0; /* truncate the rounded digit */
  1446.  
  1447.    return r;
  1448.    }
  1449.  
  1450.  
  1451. /**********************************************************************/
  1452. /* mult_a_bf10_int()                                                  */
  1453. /* r *= n                                                             */
  1454. /* dec - number of decimals, including the one extra for rounding */
  1455.  
  1456. bf10_t mult_a_bf10_int(bf10_t r, int dec, U16 n)
  1457.    {
  1458.    int signflag;
  1459.    int d, p;
  1460.    unsigned value, overflow;
  1461.    bf10_t power10;
  1462.  
  1463.    if (r[1] == 0 || n == 0)
  1464.       {
  1465.       r[1] = 0;
  1466.       return r;
  1467.       }
  1468.  
  1469.    power10 = r + dec + 1;
  1470.    p = (S16)access16(power10);
  1471.  
  1472.    signflag = r[0];  /* r[0] to be used as a padding */
  1473.    overflow = 0;
  1474.    for (d = dec; d>0; d--)
  1475.       {
  1476.       value = r[d] * n + overflow;
  1477.       r[d] = (BYTE)(value % 10);
  1478.       overflow = value / 10;
  1479.       }
  1480.    while (overflow)
  1481.       {
  1482.       p++;
  1483.       _fmemmove(r+2, r+1, dec-1);
  1484.       r[1] = (BYTE)(overflow % 10);
  1485.       overflow = overflow / 10;
  1486.       }
  1487.    set16(power10, p); /* save power of ten */
  1488.    r[0] = (BYTE)signflag; /* restore sign flag */
  1489.    return r;
  1490.    }
  1491.  
  1492. /**********************************************************************/
  1493. /* div_a_bf10_int()                                                   */
  1494. /* r /= n                                                             */
  1495. /* dec - number of decimals, including the one extra for rounding */
  1496.  
  1497. bf10_t div_a_bf10_int (bf10_t r, int dec, U16 n)
  1498.    {
  1499.    int src, dest, p;
  1500.    unsigned value, remainder;
  1501.    bf10_t power10;
  1502.  
  1503.    if (r[1] == 0 || n == 0)
  1504.       {
  1505.       r[1] = 0;
  1506.       return r;
  1507.       }
  1508.  
  1509.    power10 = r + dec + 1;
  1510.    p = (S16)access16(power10);
  1511.  
  1512.    remainder = 0;
  1513.    for (src=dest=1; src<=dec; dest++, src++)
  1514.       {
  1515.       value = 10*remainder + r[src];
  1516.       r[dest] = (BYTE)(value / n);
  1517.       remainder = value % n;
  1518.       if (dest == 1 && r[dest] == 0)
  1519.          {
  1520.          dest = 0; /* back up a digit */
  1521.          p--;      /* and decrease by a factor of 10 */
  1522.          }
  1523.       }
  1524.    for (; dest<=dec; dest++)
  1525.       {
  1526.       value = 10*remainder;
  1527.       r[dest] = (BYTE)(value / n);
  1528.       remainder = value % n;
  1529.       if (dest == 1 && r[dest] == 0)
  1530.          {
  1531.          dest = 0; /* back up a digit */
  1532.          p--;      /* and decrease by a factor of 10 */
  1533.          }
  1534.       }
  1535.  
  1536.    set16(power10, p); /* save power of ten */
  1537.    return r;
  1538.    }
  1539.  
  1540.  
  1541. /*************************************************************************/
  1542. /* bf10tostr_e()                                                         */
  1543. /* Takes a bf10 number and converts it to an ascii string, sci. notation */
  1544. /* dec - number of decimals, not including the one extra for rounding    */
  1545.  
  1546. char *bf10tostr_e(char *s, int dec, bf10_t n)
  1547.    {
  1548.    int d, p;
  1549.    bf10_t power10;
  1550.  
  1551.    if (n[1] == 0)
  1552.       {
  1553.       strcpy(s, "0.0");
  1554.       return s;
  1555.       }
  1556.  
  1557.    if (dec == 0)
  1558.        dec = decimals;
  1559.    dec++;  /* one extra byte for rounding */
  1560.    power10 = n + dec + 1;
  1561.    p = (S16)access16(power10);
  1562.  
  1563.    /* if p is negative, it is not necessary to show all the decimal places */
  1564.    if (p < 0 && dec > 8) /* 8 sounds like a reasonable value */
  1565.       {
  1566.       dec = dec + p;
  1567.       if (dec < 8) /* let's keep at least a few */
  1568.          dec = 8;
  1569.       }
  1570.  
  1571.    if (n[0] == 1) /* sign flag */
  1572.       *(s++) = '-';
  1573.    *(s++) = (char)(n[1] + '0');
  1574.    *(s++) = '.';
  1575.    for (d=2; d<=dec; d++)
  1576.       {
  1577.       *(s++) = (char)(n[d] + '0');
  1578.       }
  1579.    /* clean up trailing 0's */
  1580.    while (*(s-1) == '0')
  1581.       s--;
  1582.    if (*(s-1) == '.') /* put at least one 0 after the decimal */
  1583.       *(s++) = '0';
  1584.    sprintf(s, "e%d", p);
  1585.    return s;
  1586.    }
  1587.  
  1588. /****************************************************************************/
  1589. /* bf10tostr_f()                                                            */
  1590. /* Takes a bf10 number and converts it to an ascii string, decimal notation */
  1591.  
  1592. char *bf10tostr_f(char *s, int dec, bf10_t n)
  1593.    {
  1594.    int d, p;
  1595.    bf10_t power10;
  1596.  
  1597.    if (n[1] == 0)
  1598.       {
  1599.       strcpy(s, "0.0");
  1600.       return s;
  1601.       }
  1602.  
  1603.    if (dec == 0)
  1604.        dec = decimals;
  1605.    dec++;  /* one extra byte for rounding */
  1606.    power10 = n + dec + 1;
  1607.    p = (S16)access16(power10);
  1608.  
  1609.    /* if p is negative, it is not necessary to show all the decimal places */
  1610.    if (p < 0 && dec > 8) /* 8 sounds like a reasonable value */
  1611.       {
  1612.       dec = dec + p;
  1613.       if (dec < 8) /* let's keep at least a few */
  1614.          dec = 8;
  1615.       }
  1616.  
  1617.    if (n[0] == 1) /* sign flag */
  1618.       *(s++) = '-';
  1619.    if (p >= 0)
  1620.       {
  1621.       for (d=1; d<=p+1; d++)
  1622.          *(s++) = (char)(n[d] + '0');
  1623.       *(s++) = '.';
  1624.       for (; d<=dec; d++)
  1625.          *(s++) = (char)(n[d] + '0');
  1626.       }
  1627.    else
  1628.       {
  1629.       *(s++) = '0';
  1630.       *(s++) = '.';
  1631.       for (d=0; d>p+1; d--)
  1632.          *(s++) = '0';
  1633.       for (d=1; d<=dec; d++)
  1634.          *(s++) = (char)(n[d] + '0');
  1635.       }
  1636.  
  1637.    /* clean up trailing 0's */
  1638.    while (*(s-1) == '0')
  1639.       s--;
  1640.    if (*(s-1) == '.') /* put at least one 0 after the decimal */
  1641.       *(s++) = '0';
  1642.    *s = '\0'; /* terminating nul */
  1643.    return s;
  1644.    }
  1645.  
  1646.