home *** CD-ROM | disk | FTP | other *** search
/ Computerworld 1996 March / Computerworld_1996-03_cd.bin / idg_cd3 / grafika / fraktaly / frasr192 / bigfltc.c < prev    next >
C/C++ Source or Header  |  1995-01-12  |  24KB  |  903 lines
  1. /* bigfltc.c -  C version of routines to be eventually put in bigflta.asm */
  2.  
  3. /* Wesley Loewer's Big Floating Point Numbers. (C) 1994, Wesley B. Loewer */
  4.  
  5.  
  6. #include <memory.h>
  7. #include <string.h>
  8. #include <float.h>
  9. #include <math.h>
  10. #include "prototyp.h"
  11. #include "bignum.h"
  12. #include "bigflt.h"
  13.  
  14. /********************************************************************/
  15. /* normalize big float */
  16. bf_t norm_bf(bf_t r)
  17.     {
  18.     int scale;
  19.     BYTE hi_byte;
  20.     S16 far *rexp;
  21.  
  22.     rexp  = (S16 far *)(r+bflength);
  23.  
  24.     /* check for overflow */
  25.     hi_byte = r[bflength-1];
  26.     if (hi_byte != 0x00 && hi_byte != 0xFF)
  27.         {
  28.         _fmemmove(r, r+1, bflength-1);
  29.         r[bflength-1] = (BYTE)(hi_byte & 0x80 ? 0xFF : 0x00);
  30.         setS16(rexp,accessS16(rexp)+1);   /* exp */
  31.         }
  32.  
  33.     /* check for underflow */
  34.     else
  35.         {
  36.         for (scale = 2; scale < bflength && r[bflength-scale] == hi_byte; scale++)
  37.             ; /* do nothing */
  38.         if (scale == bflength && hi_byte == 0) /* zero */
  39.             *rexp = 0;
  40.         else
  41.             {
  42.             scale -= 2;
  43.             if (scale > 0) /* it did underflow */
  44.                 {
  45.                 _fmemmove(r+scale, r, bflength-scale-1);
  46.                 _fmemset(r, 0, scale);
  47.                 setS16(rexp,accessS16(rexp)-(S16)scale);    /* exp */
  48.                 }
  49.             }
  50.         }
  51.  
  52.     return r;
  53.     }
  54.  
  55. /********************************************************************/
  56. /* normalize big float with forced sign */
  57. /* positive = 1, force to be positive   */
  58. /*          = 0, force to be negative   */
  59. void norm_sign_bf(bf_t r, int positive)
  60.     {
  61.     norm_bf(r);
  62.     r[bflength-1] = (BYTE)(positive ? 0x00 : 0xFF);
  63.     }
  64. /******************************************************/
  65. /* adjust n1, n2 for before addition or subtraction   */
  66. /* by forcing exp's to match.                         */
  67. /* returns the value of the adjusted exponents        */
  68. S16 adjust_bf_add(bf_t n1, bf_t n2)
  69.     {
  70.     int scale, fill_byte;
  71.     S16 rexp;
  72.     S16 far *n1exp, far *n2exp;
  73.  
  74.     /* scale n1 or n2 */
  75.     /* compare exp's */
  76.     n1exp = (S16 far *)(n1+bflength);
  77.     n2exp = (S16 far *)(n2+bflength);
  78.     if (accessS16(n1exp) > accessS16(n2exp))
  79.         { /* scale n2 */
  80.         scale = accessS16(n1exp) - accessS16(n2exp); /* n1exp - n2exp */
  81.         if (scale < bflength)
  82.             {
  83.             fill_byte = is_bf_neg(n2) ? 0xFF : 0x00;
  84.             _fmemmove(n2, n2+scale, bflength-scale);
  85.             _fmemset(n2+bflength-scale, fill_byte, scale);
  86.             }
  87.         else
  88.             clear_bf(n2);
  89.         *n2exp = *n1exp; /* set exp's = */
  90.     rexp = accessS16(n2exp);
  91.         }
  92.     else if (accessS16(n1exp) < accessS16(n2exp))
  93.         { /* scale n1 */
  94.         scale = accessS16(n2exp) - accessS16(n1exp);  /* n2exp - n1exp */
  95.         if (scale < bflength)
  96.             {
  97.             fill_byte = is_bf_neg(n1) ? 0xFF : 0x00;
  98.             _fmemmove(n1, n1+scale, bflength-scale);
  99.             _fmemset(n1+bflength-scale, fill_byte, scale);
  100.             }
  101.         else
  102.             clear_bf(n1);
  103.         *n1exp = *n2exp; /* set exp's = */
  104.     rexp = accessS16(n2exp);
  105.         }
  106.     else
  107.     rexp = accessS16(n1exp);
  108.     return rexp;
  109.     }
  110.  
  111. /********************************************************************/
  112. /* r = 0 */
  113. bf_t clear_bf(bf_t r)
  114.     {
  115.     _fmemset( r, 0, bflength+2); /* set array to zero */
  116.     return r;
  117.     }
  118.  
  119. /********************************************************************/
  120. /* r = max positive value */
  121. bf_t max_bf(bf_t r)
  122.     {
  123.     inttobf(r, 1);
  124.     set16(r+bflength, (S16)(LDBL_MAX_EXP/8));
  125.     return r;
  126.     }
  127.  
  128. /********************************************************************/
  129. /* r = n */
  130. bf_t copy_bf(bf_t r, bf_t n)
  131.     {
  132.     _fmemcpy( r, n, bflength+2);
  133.     return r;
  134.     }
  135.  
  136. /***************************************************************************/
  137. /* n1 != n2 ?                                                              */
  138. /* RETURNS:                                                                */
  139. /*  if n1 == n2 returns 0                                                  */
  140. /*  if n1 > n2 returns a positive (steps left to go when mismatch occurred) */
  141. /*  if n1 < n2 returns a negative (steps left to go when mismatch occurred) */
  142.  
  143. int cmp_bf(bf_t n1, bf_t n2)
  144.     {
  145.     int i;
  146.     int sign1, sign2;
  147.     S16 far *n1exp, far *n2exp;
  148.  
  149.     /* compare signs */
  150.     sign1 = sign_bf(n1);
  151.     sign2 = sign_bf(n2);
  152.     if (sign1 > sign2)
  153.         return bflength;
  154.     else if (sign1 < sign2)
  155.         return -bflength;
  156.     /* signs are the same */
  157.  
  158.     /* compare exponents, using signed comparisons */
  159.     n1exp = (S16 far *)(n1+bflength);
  160.     n2exp = (S16 far *)(n2+bflength);
  161.     if (accessS16(n1exp) > accessS16(n2exp))
  162.         return sign1*(bflength);
  163.     else if (accessS16(n1exp) < accessS16(n2exp))
  164.         return -sign1*(bflength);
  165.  
  166.     /* To get to this point, the signs must match */
  167.     /* so unsigned comparison is ok. */
  168.     /* two bytes at a time */
  169.     for (i=bflength-2; i>=0; i-=2)
  170.         {
  171.         if (access16(n1+i) > access16(n2+i))
  172.             return i+2;
  173.         else if (access16(n1+i) < access16(n2+i))
  174.             return -(i+2);
  175.         }
  176.     return 0;
  177.     }
  178.  
  179. /********************************************************************/
  180. /* r < 0 ?                                      */
  181. /* returns 1 if negative, 0 if positive or zero */
  182. int is_bf_neg(bf_t n)
  183.     {
  184.     return (S8)n[bflength-1] < 0;
  185.     }
  186.  
  187. /********************************************************************/
  188. /* n != 0 ?                      */
  189. /* RETURNS: if n != 0 returns 1  */
  190. /*          else returns 0       */
  191. int is_bf_not_zero(bf_t n)
  192.     {
  193.     int bnl;
  194.     int retval;
  195.  
  196.     bnl = bnlength;
  197.     bnlength = bflength;
  198.     retval = is_bn_not_zero(n);
  199.     bnlength = bnl;
  200.  
  201.     return retval;
  202.     }
  203.  
  204. /********************************************************************/
  205. /* r = n1 + n2 */
  206. /* SIDE-EFFECTS: n1 and n2 can be "de-normalized" and lose precision */
  207. bf_t unsafe_add_bf(bf_t r, bf_t n1, bf_t n2)
  208.     {
  209.     int bnl;
  210.     S16 far *rexp;
  211.  
  212.     if (is_bf_zero(n1))
  213.         {
  214.         copy_bf(r, n2);
  215.         return r;
  216.         }
  217.     if (is_bf_zero(n2))
  218.         {
  219.         copy_bf(r, n1);
  220.         return r;
  221.         }
  222.  
  223.     rexp = (S16 far *)(r+bflength);
  224.     setS16(rexp,adjust_bf_add(n1, n2));
  225.  
  226.     bnl = bnlength;
  227.     bnlength = bflength;
  228.     add_bn(r, n1, n2);
  229.     bnlength = bnl;
  230.  
  231.     norm_bf(r);
  232.     return r;
  233.     }
  234.  
  235. /********************************************************************/
  236. /* r += n */
  237. bf_t unsafe_add_a_bf(bf_t r, bf_t n)
  238.     {
  239.     int bnl;
  240.  
  241.     if (is_bf_zero(r))
  242.         {
  243.         copy_bf(r, n);
  244.         return r;
  245.         }
  246.     if (is_bf_zero(n))
  247.         {
  248.         return r;
  249.         }
  250.  
  251.     adjust_bf_add(r, n);
  252.  
  253.     bnl = bnlength;
  254.     bnlength = bflength;
  255.     add_a_bn(r, n);
  256.     bnlength = bnl;
  257.  
  258.     norm_bf(r);
  259.  
  260.     return r;
  261.     }
  262.  
  263. /********************************************************************/
  264. /* r = n1 - n2 */
  265. /* SIDE-EFFECTS: n1 and n2 can be "de-normalized" and lose precision */
  266. bf_t unsafe_sub_bf(bf_t r, bf_t n1, bf_t n2)
  267.     {
  268.     int bnl;
  269.     S16 far *rexp;
  270.  
  271.     if (is_bf_zero(n1))
  272.         {
  273.         neg_bf(r, n2);
  274.         return r;
  275.         }
  276.     if (is_bf_zero(n2))
  277.         {
  278.         copy_bf(r, n1);
  279.         return r;
  280.         }
  281.  
  282.     rexp = (S16 far *)(r+bflength);
  283.     setS16(rexp,adjust_bf_add(n1, n2));
  284.  
  285.     bnl = bnlength;
  286.     bnlength = bflength;
  287.     sub_bn(r, n1, n2);
  288.     bnlength = bnl;
  289.  
  290.     norm_bf(r);
  291.     return r;
  292.     }
  293.  
  294. /********************************************************************/
  295. /* r -= n */
  296. bf_t unsafe_sub_a_bf(bf_t r, bf_t n)
  297.     {
  298.     int bnl;
  299.  
  300.     if (is_bf_zero(r))
  301.         {
  302.         neg_bf(r,n);
  303.         return r;
  304.         }
  305.     if (is_bf_zero(n))
  306.         {
  307.         return r;
  308.         }
  309.  
  310.     adjust_bf_add(r, n);
  311.  
  312.     bnl = bnlength;
  313.     bnlength = bflength;
  314.     sub_a_bn(r, n);
  315.     bnlength = bnl;
  316.  
  317.     norm_bf(r);
  318.     return r;
  319.     }
  320.  
  321. /********************************************************************/
  322. /* r = -n */
  323. bf_t neg_bf(bf_t r, bf_t n)
  324.     {
  325.     int bnl;
  326.     S16 far *rexp, far *nexp;
  327.  
  328.     rexp = (S16 far *)(r+bflength);
  329.     nexp = (S16 far *)(n+bflength);
  330.     *rexp = *nexp;
  331.  
  332.     bnl = bnlength;
  333.     bnlength = bflength;
  334.     neg_bn(r, n);
  335.     bnlength = bnl;
  336.  
  337.     norm_bf(r);
  338.     return r;
  339.     }
  340.  
  341. /********************************************************************/
  342. /* r *= -1 */
  343. bf_t neg_a_bf(bf_t r)
  344.     {
  345.     int bnl;
  346.  
  347.     bnl = bnlength;
  348.     bnlength = bflength;
  349.     neg_a_bn(r);
  350.     bnlength = bnl;
  351.  
  352.     norm_bf(r);
  353.     return r;
  354.     }
  355.  
  356. /********************************************************************/
  357. /* r = 2*n */
  358. bf_t double_bf(bf_t r, bf_t n)
  359.     {
  360.     int bnl;
  361.     S16 far *rexp, far *nexp;
  362.  
  363.     rexp = (S16 far *)(r+bflength);
  364.     nexp = (S16 far *)(n+bflength);
  365.     *rexp = *nexp;
  366.  
  367.     bnl = bnlength;
  368.     bnlength = bflength;
  369.     double_bn(r, n);
  370.     bnlength = bnl;
  371.  
  372.     norm_bf(r);
  373.     return r;
  374.     }
  375.  
  376. /********************************************************************/
  377. /* r *= 2 */
  378. bf_t double_a_bf(bf_t r)
  379.     {
  380.     int bnl;
  381.  
  382.     bnl = bnlength;
  383.     bnlength = bflength;
  384.     double_a_bn(r);
  385.     bnlength = bnl;
  386.  
  387.     norm_bf(r);
  388.     return r;
  389.     }
  390.  
  391. /********************************************************************/
  392. /* r = n/2 */
  393. bf_t half_bf(bf_t r, bf_t n)
  394.     {
  395.     int bnl;
  396.     S16 far *rexp, far *nexp;
  397.  
  398.     rexp = (S16 far *)(r+bflength);
  399.     nexp = (S16 far *)(n+bflength);
  400.     *rexp = *nexp;
  401.  
  402.     bnl = bnlength;
  403.     bnlength = bflength;
  404.     half_bn(r, n);
  405.     bnlength = bnl;
  406.  
  407.     norm_bf(r);
  408.     return r;
  409.     }
  410.  
  411. /********************************************************************/
  412. /* r /= 2 */
  413. bf_t half_a_bf(bf_t r)
  414.     {
  415.     int bnl;
  416.  
  417.     bnl = bnlength;
  418.     bnlength = bflength;
  419.     half_a_bn(r);
  420.     bnlength = bnl;
  421.  
  422.     norm_bf(r);
  423.     return r;
  424.     }
  425.  
  426. /************************************************************************/
  427. /* r = n1 * n2                                                          */
  428. /* Note: r will be a double wide result, 2*bflength                     */
  429. /*       n1 and n2 can be the same pointer                              */
  430. /* SIDE-EFFECTS: n1 and n2 are changed to their absolute values         */
  431. bf_t unsafe_full_mult_bf(bf_t r, bf_t n1, bf_t n2)
  432.     {
  433.     int bnl, dbfl;
  434.     S16 far *rexp, far *n1exp, far *n2exp;
  435.  
  436.     if (is_bf_zero(n1) || is_bf_zero(n2) )
  437.         {
  438.         bflength <<= 1;
  439.         clear_bf(r);
  440.         bflength >>= 1;
  441.         return r;
  442.         }
  443.  
  444.     dbfl = 2*bflength; /* double width bflength */
  445.     rexp  = (S16 far *)(r+dbfl); /* note: 2*bflength */
  446.     n1exp = (S16 far *)(n1+bflength);
  447.     n2exp = (S16 far *)(n2+bflength);
  448.     /* add exp's */
  449.     *rexp = (S16)(accessS16(n1exp) + accessS16(n2exp));
  450.  
  451.     bnl = bnlength;
  452.     bnlength = bflength;
  453.     unsafe_full_mult_bn(r, n1, n2);
  454.     bnlength = bnl;
  455.  
  456.     /* handle normalizing full mult on individual basis */
  457.  
  458.     return r;
  459.     }
  460.  
  461. /************************************************************************/
  462. /* r = n1 * n2 calculating only the top rlength bytes                   */
  463. /* Note: r will be of length rlength                                    */
  464. /*       2*bflength <= rlength < bflength                               */
  465. /*       n1 and n2 can be the same pointer                              */
  466. /* SIDE-EFFECTS: n1 and n2 are changed to their absolute values         */
  467. bf_t unsafe_mult_bf(bf_t r, bf_t n1, bf_t n2)
  468.     {
  469.     int positive;
  470.     int bnl, bfl, rl;
  471.     int rexp;
  472.     S16 far *n1exp, far *n2exp;
  473.  
  474.     if (is_bf_zero(n1) || is_bf_zero(n2) )
  475.         {
  476.         clear_bf(r);
  477.         return r;
  478.         }
  479.  
  480.     n1exp = (S16 far *)(n1+bflength);
  481.     n2exp = (S16 far *)(n2+bflength);
  482.     /* add exp's */
  483.     rexp = accessS16(n1exp) + accessS16(n2exp);
  484.  
  485.     positive = is_bf_neg(n1) == is_bf_neg(n2); /* are they the same sign? */
  486.  
  487.     bnl = bnlength;
  488.     bnlength = bflength;
  489.     rl = rlength;
  490.     rlength = rbflength;
  491.     unsafe_mult_bn(r, n1, n2);
  492.     bnlength = bnl;
  493.     rlength = rl;
  494.  
  495.     bfl = bflength;
  496.     bflength = rbflength;
  497.     set16(r+bflength, (S16)(rexp+2)); /* adjust after mult */
  498.     norm_sign_bf(r, positive);
  499.     bflength = bfl;
  500.     _fmemmove(r, r+padding, bflength+2); /* shift back */
  501.  
  502.     return r;
  503.     }
  504.  
  505. /************************************************************************/
  506. /* r = n^2                                                              */
  507. /*   because of the symmetry involved, n^2 is much faster than n*n      */
  508. /*   for a bignumber of length l                                        */
  509. /*      n*n takes l^2 multiplications                                   */
  510. /*      n^2 takes (l^2+l)/2 multiplications                             */
  511. /*          which is about 1/2 n*n as l gets large                      */
  512. /*  uses the fact that (a+b+c+...)^2 = (a^2+b^2+c^2+...)+2(ab+ac+bc+...)*/
  513. /*                                                                      */
  514. /* SIDE-EFFECTS: n is changed to its absolute value                     */
  515. bf_t unsafe_full_square_bf(bf_t r, bf_t n)
  516.     {
  517.     int bnl, dbfl;
  518.     S16 far *rexp, far *nexp;
  519.  
  520.     if (is_bf_zero(n))
  521.         {
  522.         bflength <<= 1;
  523.         clear_bf(r);
  524.         bflength >>= 1;
  525.         return r;
  526.         }
  527.  
  528.     dbfl = 2*bflength; /* double width bflength */
  529.     rexp  = (S16 far *)(r+dbfl); /* note: 2*bflength */
  530.     nexp = (S16 far *)(n+bflength);
  531.     setS16(rexp, 2 * accessS16(nexp));
  532.  
  533.     bnl = bnlength;
  534.     bnlength = bflength;
  535.     unsafe_full_square_bn(r, n);
  536.     bnlength = bnl;
  537.  
  538.     /* handle normalizing full mult on individual basis */
  539.  
  540.     return r;
  541.     }
  542.  
  543.  
  544. /************************************************************************/
  545. /* r = n^2                                                              */
  546. /*   because of the symmetry involved, n^2 is much faster than n*n      */
  547. /*   for a bignumber of length l                                        */
  548. /*      n*n takes l^2 multiplications                                   */
  549. /*      n^2 takes (l^2+l)/2 multiplications                             */
  550. /*          which is about 1/2 n*n as l gets large                      */
  551. /*  uses the fact that (a+b+c+...)^2 = (a^2+b^2+c^2+...)+2(ab+ac+bc+...)*/
  552. /*                                                                      */
  553. /* Note: r will be of length rlength                                    */
  554. /*       2*bflength >= rlength > bflength                               */
  555. /* SIDE-EFFECTS: n is changed to its absolute value                     */
  556. bf_t unsafe_square_bf(bf_t r, bf_t n)
  557.     {
  558.     int bnl, bfl, rl;
  559.     int rexp;
  560.     S16 far *nexp;
  561.  
  562.     if (is_bf_zero(n))
  563.         {
  564.         clear_bf(r);
  565.         return r;
  566.         }
  567.  
  568.     nexp = (S16 far *)(n+bflength);
  569.     rexp = (S16)(2 * accessS16(nexp));
  570.  
  571.     bnl = bnlength;
  572.     bnlength = bflength;
  573.     rl = rlength;
  574.     rlength = rbflength;
  575.     unsafe_square_bn(r, n);
  576.     bnlength = bnl;
  577.     rlength = rl;
  578.  
  579.     bfl = bflength;
  580.     bflength = rbflength;
  581.     set16(r+bflength, (S16)(rexp+2)); /* adjust after mult */
  582.  
  583.     norm_sign_bf(r, 1);
  584.     bflength = bfl;
  585.     _fmemmove(r, r+padding, bflength+2); /* shift back */
  586.  
  587.     return r;
  588.     }
  589.  
  590. /********************************************************************/
  591. /* r = n * u  where u is an unsigned integer */
  592. /* SIDE-EFFECTS: n can be "de-normalized" and lose precision */
  593. bf_t unsafe_mult_bf_int(bf_t r, bf_t n, U16 u)
  594.     {
  595.     int positive;
  596.     int bnl;
  597.     S16 far *rexp, far *nexp;
  598.  
  599.     rexp = (S16 far *)(r+bflength);
  600.     nexp = (S16 far *)(n+bflength);
  601.     *rexp = *nexp;
  602.  
  603.     positive = !is_bf_neg(n);
  604.  
  605. /*
  606. if u > 0x00FF, then the integer part of the mantissa will overflow the
  607. 2 byte (16 bit) integer size.  Therefore, make adjustment before
  608. multiplication is performed.
  609. */
  610.     if (u > 0x00FF)
  611.         { /* un-normalize n */
  612.         _fmemmove(n, n+1, bflength-1);  /* this sign extends as well */
  613.     setS16(rexp,accessS16(rexp)+1);
  614.         }
  615.  
  616.     bnl = bnlength;
  617.     bnlength = bflength;
  618.     mult_bn_int(r, n, u);
  619.     bnlength = bnl;
  620.  
  621.     norm_sign_bf(r, positive);
  622.     return r;
  623.     }
  624.  
  625. /********************************************************************/
  626. /* r *= u  where u is an unsigned integer */
  627. bf_t mult_a_bf_int(bf_t r, U16 u)
  628.     {
  629.     int positive;
  630.     int bnl;
  631.     S16 far *rexp;
  632.  
  633.     rexp = (S16 far *)(r+bflength);
  634.     positive = !is_bf_neg(r);
  635.  
  636. /*
  637. if u > 0x00FF, then the integer part of the mantissa will overflow the
  638. 2 byte (16 bit) integer size.  Therefore, make adjustment before
  639. multiplication is performed.
  640. */
  641.     if (u > 0x00FF)
  642.         { /* un-normalize n */
  643.         _fmemmove(r, r+1, bflength-1);  /* this sign extends as well */
  644.     setS16(rexp,accessS16(rexp)+1);
  645.         }
  646.  
  647.     bnl = bnlength;
  648.     bnlength = bflength;
  649.     mult_a_bn_int(r, u);
  650.     bnlength = bnl;
  651.  
  652.     norm_sign_bf(r, positive);
  653.     return r;
  654.     }
  655.  
  656. /********************************************************************/
  657. /* r = n / u  where u is an unsigned integer */
  658. bf_t unsafe_div_bf_int(bf_t r, bf_t n,  U16 u)
  659.     {
  660.     int bnl;
  661.     S16 far *rexp, far *nexp;
  662.  
  663.     if (u == 0) /* division by zero */
  664.         {
  665.         max_bf(r);
  666.         if (is_bf_neg(n))
  667.             neg_a_bf(r);
  668.         return r;
  669.         }
  670.  
  671.     rexp = (S16 far *)(r+bflength);
  672.     nexp = (S16 far *)(n+bflength);
  673.     *rexp = *nexp;
  674.  
  675.     bnl = bnlength;
  676.     bnlength = bflength;
  677.     unsafe_div_bn_int(r, n, u);
  678.     bnlength = bnl;
  679.  
  680.     norm_bf(r);
  681.     return r;
  682.     }
  683.  
  684. /********************************************************************/
  685. /* r /= u  where u is an unsigned integer */
  686. bf_t div_a_bf_int(bf_t r, U16 u)
  687.     {
  688.     int bnl;
  689.  
  690.     if (u == 0) /* division by zero */
  691.         {
  692.         if (is_bf_neg(r))
  693.             {
  694.             max_bf(r);
  695.             neg_a_bf(r);
  696.             }
  697.         else
  698.             {
  699.             max_bf(r);
  700.             }
  701.         return r;
  702.         }
  703.  
  704.     bnl = bnlength;
  705.     bnlength = bflength;
  706.     div_a_bn_int(r, u);
  707.     bnlength = bnl;
  708.  
  709.     norm_bf(r);
  710.     return r;
  711.     }
  712.  
  713. /********************************************************************/
  714. /* extracts the mantissa and exponent of f                          */
  715. /* finds m and n such that 1<=|m|<b and f = m*b^n                   */
  716. /* n is stored in *exp_ptr and m is returned, sort of like frexp()  */
  717. LDBL extract_value(LDBL f, LDBL b, int *exp_ptr)
  718.    {
  719.    int n;
  720.    LDBL af, ff, orig_b;
  721.    LDBL value[15];
  722.    unsigned powertwo;
  723.  
  724.    if (b <= 0 || f == 0)
  725.       {
  726.       *exp_ptr = 0;
  727.       return 0;
  728.       }
  729.  
  730.    orig_b = b;
  731.    af = f >= 0 ? f: -f;     /* abs value */
  732.    ff = af > 1 ? af : 1/af;
  733.    n = 0;
  734.    powertwo = 1;
  735.    while (b < ff)
  736.       {
  737.       value[n] = b;
  738.       n++;
  739.       powertwo <<= 1;
  740.       b *= b;
  741.       }
  742.  
  743.    *exp_ptr = 0;
  744.    for (; n > 0; n--)
  745.       {
  746.       powertwo >>= 1;
  747.       if (value[n-1] < ff)
  748.          {
  749.          ff /= value[n-1];
  750.          *exp_ptr += powertwo;
  751.          }
  752.       }
  753.    if (f < 0)
  754.       ff = -ff;
  755.    if (af < 1)
  756.       {
  757.       ff = orig_b/ff;
  758.       *exp_ptr = -*exp_ptr - 1;
  759.       }
  760.  
  761.    return ff;
  762.    }
  763.  
  764. /********************************************************************/
  765. /* calculates and returns the value of f*b^n                        */
  766. /* sort of like ldexp()                                             */
  767. LDBL scale_value( LDBL f, LDBL b , int n )
  768.    {
  769.    LDBL total=1;
  770.    int an;
  771.  
  772.    if (b == 0 || f == 0)
  773.       return 0;
  774.  
  775.    if (n == 0)
  776.       return f;
  777.  
  778.    an = abs(n);
  779.  
  780.    while (an != 0)
  781.       {
  782.       if (an & 0x0001)
  783.          total *= b;
  784.       b *= b;
  785.       an >>= 1;
  786.       }
  787.  
  788.    if (n > 0)
  789.       f *= total;
  790.    else /* n < 0 */
  791.       f /= total;
  792.    return f;
  793.    }
  794.  
  795. /********************************************************************/
  796. /* extracts the mantissa and exponent of f                          */
  797. /* finds m and n such that 1<=|m|<10 and f = m*10^n                 */
  798. /* n is stored in *exp_ptr and m is returned, sort of like frexp()  */
  799. LDBL extract_10(LDBL f, int *exp_ptr)
  800.    {
  801.    return extract_value(f, 10, exp_ptr);
  802.    }
  803.  
  804. /********************************************************************/
  805. /* calculates and returns the value of f*10^n                       */
  806. /* sort of like ldexp()                                             */
  807. LDBL scale_10( LDBL f, int n )
  808.    {
  809.    return scale_value( f, 10, n );
  810.    }
  811.  
  812.  
  813. #ifdef USE_BIGNUM_C_CODE
  814.  
  815. /********************************************************************/
  816. /* extracts the mantissa and exponent of f                          */
  817. /* finds m and n such that 1<=|m|<256 and f = m*256^n               */
  818. /* n is stored in *exp_ptr and m is returned, sort of like frexp()  */
  819. LDBL extract_256(LDBL f, int *exp_ptr)
  820.    {
  821.    return extract_value(f, 256, exp_ptr);
  822.    }
  823.  
  824. /********************************************************************/
  825. /* calculates and returns the value of f*256^n                      */
  826. /* sort of like ldexp()                                             */
  827. /*                                                                  */
  828. /* n must be in the range -2^12 <= n < 2^12 (2^12=4096),            */
  829. /* which should not be a problem                                    */
  830. LDBL scale_256( LDBL f, int n )
  831.    {
  832.    return scale_value( f, 256, n );
  833.    }
  834.  
  835. /*********************************************************************/
  836. /*  b = f                                                            */
  837. /*  Converts a double to a bigfloat                                  */
  838. bf_t floattobf(bf_t r, LDBL f)
  839.     {
  840.     int power;
  841.     int bnl, il;
  842.     if (f == 0)
  843.         {
  844.         clear_bf(r);
  845.         return r;
  846.         }
  847.  
  848.     /* remove the exp part */
  849.     f = extract_256(f, &power);
  850.  
  851.     bnl = bnlength;
  852.     bnlength = bflength;
  853.     il = intlength;
  854.     intlength = 2;
  855.     floattobn(r, f);
  856.     bnlength = bnl;
  857.     intlength = il;
  858.  
  859.     set16(r + bflength, (S16)power); /* exp */
  860.  
  861.     return r;
  862.     }
  863.  
  864. /*********************************************************************/
  865. /*  b = f                                                            */
  866. /*  Converts a double to a bigfloat                                  */
  867. bf_t floattobf1(bf_t r, LDBL f)
  868.     {
  869.     char msg[80];
  870. #ifdef USE_LONG_DOUBLE
  871.     sprintf(msg,"%-.22Le",f);
  872. #else
  873.     sprintf(msg,"%-.22le",f);
  874. #endif    
  875.     strtobf(r,msg); 
  876.     return r;
  877.     }
  878.  
  879. /*********************************************************************/
  880. /*  f = b                                                            */
  881. /*  Converts a bigfloat to a double                                 */
  882. LDBL bftofloat(bf_t n)
  883.     {
  884.     int power;
  885.     int bnl, il;
  886.     LDBL f;
  887.  
  888.     bnl = bnlength;
  889.     bnlength = bflength;
  890.     il = intlength;
  891.     intlength = 2;
  892.     f = bntofloat(n);
  893.     bnlength = bnl;
  894.     intlength = il;
  895.  
  896.     power = (S16)access16(n + bflength);
  897.     f = scale_256(f,power);
  898.  
  899.     return f;
  900.     }
  901.  
  902. #endif /* USE_BIGNUM_C_CODE */
  903.