home *** CD-ROM | disk | FTP | other *** search
/ Computerworld 1996 March / Computerworld_1996-03_cd.bin / idg_cd3 / grafika / fraktaly / frasr192 / bignumc.c < prev    next >
C/C++ Source or Header  |  1995-01-12  |  29KB  |  825 lines
  1. /* bignumc.c - C routines equivalent to ASM routines in bignuma.asm */
  2. /* Wesley Loewer's Big Numbers.        (C) 1994, Wesley B. Loewer   */
  3.  
  4. #include <memory.h>
  5. #include "prototyp.h"
  6.  
  7. #ifndef _BIGNUM_H
  8. #include "bignum.h"
  9. #endif
  10.  
  11. #ifdef sign
  12. #undef sign
  13. #endif
  14.  
  15.  
  16. /********************************************************************
  17.  The following code contains the C versions of the routines from the
  18.  file BIGNUMA.ASM.  It is provided here for portibility and for clarity.
  19. *********************************************************************/
  20.  
  21. /********************************************************************/
  22. /* r = 0 */
  23. bn_t clear_bn(bn_t r)
  24.     {
  25.     return _fmemset( r, 0, bnlength); /* set array to zero */
  26.     }
  27.  
  28. /********************************************************************/
  29. /* r = max positive value */
  30. bn_t max_bn(bn_t r)
  31.     {
  32.     _fmemset( r, 0xFF, bnlength-1); /* set to max values */
  33.     r[bnlength-1] = 0x7F;  /* turn off the sign bit */
  34.     return r;
  35.     }
  36.  
  37. /********************************************************************/
  38. /* r = n */
  39. bn_t copy_bn(bn_t r, bn_t n)
  40.     {
  41.     return _fmemcpy( r, n, bnlength);
  42.     }
  43.  
  44. /***************************************************************************/
  45. /* n1 != n2 ?                                                              */
  46. /* RETURNS:                                                                */
  47. /*  if n1 == n2 returns 0                                                  */
  48. /*  if n1 > n2 returns a positive (steps left to go when mismatch occured) */
  49. /*  if n1 < n2 returns a negative (steps left to go when mismatch occured) */
  50. int cmp_bn(bn_t n1, bn_t n2)
  51.     {
  52.     int i;
  53.  
  54.     /* two bytes at a time */
  55.     /* signed comparison for msb */
  56.     if (accessS16((S16 far *)(n1+bnlength-2)) >
  57.         accessS16((S16 far *)(n2+bnlength-2)))
  58.         return bnlength;
  59.     else if (accessS16((S16 far *)(n1+bnlength-2)) <
  60.         accessS16((S16 far *)(n2+bnlength-2)))
  61.         return -(bnlength);
  62.  
  63.     /* unsigned comparison for the rest */
  64.     for (i=bnlength-4; i>=0; i-=2)
  65.         {
  66.         if (access16(n1+i) > access16(n2+i))
  67.             return i+2;
  68.         else if (access16(n1+i) < access16(n2+i))
  69.             return -(i+2);
  70.         }
  71.     return 0;
  72.     }
  73.  
  74. /********************************************************************/
  75. /* r < 0 ?                                      */
  76. /* returns 1 if negative, 0 if positive or zero */
  77. int is_bn_neg(bn_t n)
  78.     {
  79.     return (S8)n[bnlength-1] < 0;
  80.     }
  81.  
  82. /********************************************************************/
  83. /* n != 0 ?                      */
  84. /* RETURNS: if n != 0 returns 1  */
  85. /*          else returns 0       */
  86. int is_bn_not_zero(bn_t n)
  87.     {
  88.     int i;
  89.  
  90.     /* two bytes at a time */
  91.     for (i=0; i<bnlength; i+=2)
  92.         if (access16(n+i) != 0)
  93.             return 1;
  94.     return 0;
  95.     }
  96.  
  97. /********************************************************************/
  98. /* r = n1 + n2 */
  99. bn_t add_bn(bn_t r, bn_t n1, bn_t n2)
  100.     {
  101.     int i;
  102.     U32 sum=0;
  103.  
  104.     /* two bytes at a time */
  105.     for (i=0; i<bnlength; i+=2)
  106.         {
  107.         sum += (U32)access16(n1+i) + (U32)access16(n2+i); /* add 'em up */
  108.         set16(r+i, (U16)sum);   /* store the lower 2 bytes */
  109.         sum >>= 16; /* shift the overflow for next time */
  110.         }
  111.     return r;
  112.     }
  113.  
  114. /********************************************************************/
  115. /* r += n */
  116. bn_t add_a_bn(bn_t r, bn_t n)
  117.     {
  118.     int i;
  119.     U32 sum=0;
  120.  
  121.     /* two bytes at a time */
  122.     for (i=0; i<bnlength; i+=2)
  123.         {
  124.         sum += (U32)access16(r+i) + (U32)access16(n+i); /* add 'em up */
  125.         set16(r+i, (U16)sum);   /* store the lower 2 bytes */
  126.         sum >>= 16; /* shift the overflow for next time */
  127.         }
  128.     return r;
  129.     }
  130.  
  131. /********************************************************************/
  132. /* r = n1 - n2 */
  133. bn_t sub_bn(bn_t r, bn_t n1, bn_t n2)
  134.     {
  135.     int i;
  136.     U32 diff=0;
  137.  
  138.     /* two bytes at a time */
  139.     for (i=0; i<bnlength; i+=2)
  140.         {
  141.         diff = (U32)access16(n1+i) - ((U32)access16(n2+i)-(S32)(S16)diff); /* subtract with borrow */
  142.         set16(r+i, (U16)diff);   /* store the lower 2 bytes */
  143.         diff >>= 16; /* shift the underflow for next time */
  144.         }
  145.     return r;
  146.     }
  147.  
  148. /********************************************************************/
  149. /* r -= n */
  150. bn_t sub_a_bn(bn_t r, bn_t n)
  151.     {
  152.     int i;
  153.     U32 diff=0;
  154.  
  155.     /* two bytes at a time */
  156.     for (i=0; i<bnlength; i+=2)
  157.         {
  158.         diff = (U32)access16(r+i) - ((U32)access16(n+i)-(S32)(S16)diff); /* subtract with borrow */
  159.         set16(r+i, (U16)diff);   /* store the lower 2 bytes */
  160.         diff >>= 16; /* shift the underflow for next time */
  161.         }
  162.     return r;
  163.     }
  164.  
  165. /********************************************************************/
  166. /* r = -n */
  167. bn_t neg_bn(bn_t r, bn_t n)
  168.     {
  169.     int i;
  170.     U16 t_short;
  171.     U32 neg=1; /* to get the 2's complement started */
  172.  
  173.     /* two bytes at a time */
  174.     for (i=0; neg != 0 && i<bnlength; i+=2)
  175.         {
  176.         t_short = ~access16(n+i);
  177.         neg += ((U32)t_short); /* two's complement */
  178.         set16(r+i, (U16)neg);   /* store the lower 2 bytes */
  179.         neg >>= 16; /* shift the sign bit for next time */
  180.         }
  181.     /* if neg was 0, then just "not" the rest */
  182.     for (; i<bnlength; i+=2)
  183.         { /* notice that access16() and set16() are not needed here */
  184.         *(U16 far *)(r+i) = ~*(U16 far *)(n+i); /* toggle all the bits */
  185.         }
  186.     return r;
  187.     }
  188.  
  189. /********************************************************************/
  190. /* r *= -1 */
  191. bn_t neg_a_bn(bn_t r)
  192.     {
  193.     int i;
  194.     U16 t_short;
  195.     U32 neg=1; /* to get the 2's complement started */
  196.  
  197.     /* two bytes at a time */
  198.     for (i=0; neg != 0 && i<bnlength; i+=2)
  199.         {
  200.     t_short = ~access16(r+i);
  201.         neg += ((U32)t_short); /* two's complement */
  202.         set16(r+i, (U16)neg);   /* store the lower 2 bytes */
  203.         neg >>= 16; /* shift the sign bit for next time */
  204.         }
  205.     /* if neg was 0, then just "not" the rest */
  206.     for (; i<bnlength; i+=2)
  207.         { /* notice that access16() and set16() are not needed here */
  208.         *(U16 far *)(r+i) = ~*(U16 far *)(r+i); /* toggle all the bits */
  209.         }
  210.     return r;
  211.     }
  212.  
  213. /********************************************************************/
  214. /* r = 2*n */
  215. bn_t double_bn(bn_t r, bn_t n)
  216.     {
  217.     int i;
  218.     U32 prod=0;
  219.  
  220.     /* two bytes at a time */
  221.     for (i=0; i<bnlength; i+=2)
  222.         {
  223.         prod += (U32)access16(n+i)<<1 ; /* double it */
  224.         set16(r+i, (U16)prod);   /* store the lower 2 bytes */
  225.         prod >>= 16; /* shift the overflow for next time */
  226.         }
  227.     return r;
  228.     }
  229.  
  230. /********************************************************************/
  231. /* r *= 2 */
  232. bn_t double_a_bn(bn_t r)
  233.     {
  234.     int i;
  235.     U32 prod=0;
  236.  
  237.     /* two bytes at a time */
  238.     for (i=0; i<bnlength; i+=2)
  239.         {
  240.         prod += (U32)access16(r+i)<<1 ; /* double it */
  241.         set16(r+i, (U16)prod);   /* store the lower 2 bytes */
  242.         prod >>= 16; /* shift the overflow for next time */
  243.         }
  244.     return r;
  245.     }
  246.  
  247. /********************************************************************/
  248. /* r = n/2 */
  249. bn_t half_bn(bn_t r, bn_t n)
  250.     {
  251.     int i;
  252.     U32 quot=0;
  253.  
  254.     /* two bytes at a time */
  255.  
  256.     /* start with an arithmetic shift */
  257.     i=bnlength-2;
  258.     quot += (U32)(((S32)(S16)access16(n+i)<<16)>>1) ; /* shift to upper 2 bytes and half it */
  259.     set16(r+i, (U16)(quot>>16));   /* store the upper 2 bytes */
  260.     quot <<= 16; /* shift the underflow for next time */
  261.  
  262.     for (i=bnlength-4; i>=0; i-=2)
  263.         {
  264.         /* looks wierd, but properly sign extends argument */
  265.         quot += (U32)(((U32)access16(n+i)<<16)>>1) ; /* shift to upper 2 bytes and half it */
  266.         set16(r+i, (U16)(quot>>16));   /* store the upper 2 bytes */
  267.         quot <<= 16; /* shift the underflow for next time */
  268.         }
  269.  
  270.     return r;
  271.     }
  272.  
  273. /********************************************************************/
  274. /* r /= 2 */
  275. bn_t half_a_bn(bn_t r)
  276.     {
  277.     int i;
  278.     U32 quot=0;
  279.  
  280.     /* two bytes at a time */
  281.  
  282.     /* start with an arithmetic shift */
  283.     i=bnlength-2;
  284.     quot += (U32)(((S32)(S16)access16(r+i)<<16)>>1) ; /* shift to upper 2 bytes and half it */
  285.     set16(r+i, (U16)(quot>>16));   /* store the upper 2 bytes */
  286.     quot <<= 16; /* shift the underflow for next time */
  287.  
  288.     for (i=bnlength-4; i>=0; i-=2)
  289.         {
  290.         /* looks wierd, but properly sign extends argument */
  291.         quot += (U32)(((U32)(U16)access16(r+i)<<16)>>1) ; /* shift to upper 2 bytes and half it */
  292.         set16(r+i, (U16)(quot>>16));   /* store the upper 2 bytes */
  293.         quot <<= 16; /* shift the underflow for next time */
  294.         }
  295.     return r;
  296.     }
  297.  
  298. /************************************************************************/
  299. /* r = n1 * n2                                                          */
  300. /* Note: r will be a double wide result, 2*bnlength                     */
  301. /*       n1 and n2 can be the same pointer                              */
  302. /* SIDE-EFFECTS: n1 and n2 are changed to their absolute values         */
  303. bn_t unsafe_full_mult_bn(bn_t r, bn_t n1, bn_t n2)
  304.     {
  305.     int sign1, sign2, samevar;
  306.     int i, j, k, steps, doublesteps, carry_steps;
  307.     bn_t n1p, n2p;      /* pointers for n1, n2 */
  308.     bn_t rp1, rp2, rp3; /* pointers for r */
  309.     U32 prod, sum;
  310.  
  311.     if ((sign1 = is_bn_neg(n1)) != 0) /* =, not == */
  312.         neg_a_bn(n1);
  313.     if (!(samevar = (n1 == n2))) /* check to see if they're the same pointer */
  314.         if ((sign2 = is_bn_neg(n2)) != 0) /* =, not == */
  315.             neg_a_bn(n2);
  316.  
  317.     n1p = n1;
  318.     steps = bnlength>>1; /* two bytes at a time */
  319.     carry_steps = doublesteps = (steps<<1) - 2;
  320.     bnlength <<= 1;
  321.     clear_bn(r);        /* double width */
  322.     bnlength >>= 1;
  323.     rp1 = rp2 = r;
  324.     for (i = 0; i < steps; i++)
  325.         {
  326.         n2p = n2;
  327.         for (j = 0; j < steps; j++)
  328.             {
  329.             prod = (U32)access16(n1p) * (U32)access16(n2p); /* U16*U16=U32 */
  330.             sum = (U32)access16(rp2) + prod; /* add to previous, including overflow */
  331.             set16(rp2, (U16)sum); /* save the lower 2 bytes */
  332.             sum >>= 16;             /* keep just the upper 2 bytes */
  333.             rp3 = rp2 + 2;          /* move over 2 bytes */
  334.             sum += access16(rp3);     /* add what was the upper two bytes */
  335.             set16(rp3 ,(U16)sum); /* save what was the upper two bytes */
  336.             sum >>= 16;             /* keep just the overflow */
  337.             for (k=0; sum != 0 && k<carry_steps; k++)
  338.                 {
  339.                 rp3 += 2;               /* move over 2 bytes */
  340.                 sum += access16(rp3);     /* add to what was the overflow */
  341.                 set16(rp3, (U16)sum); /* save what was the overflow */
  342.                 sum >>= 16;             /* keep just the new overflow */
  343.                 }
  344.             n2p += 2;       /* to next word */
  345.             rp2 += 2;
  346.             carry_steps--;  /* use one less step */
  347.             }
  348.         n1p += 2;           /* to next word */
  349.         rp2 = rp1 += 2;
  350.         carry_steps = --doublesteps; /* decrease doubles steps and reset carry_steps */
  351.         }
  352.  
  353.     /* if they were the same or same sign, the product must be positive */
  354.     if (!samevar && sign1 != sign2)
  355.         {
  356.         bnlength <<= 1;         /* for a double wide number */
  357.         neg_a_bn(r);
  358.         bnlength >>= 1; /* restore bnlength */
  359.         }
  360.     return r;
  361.     }
  362.  
  363. /************************************************************************/
  364. /* r = n1 * n2 calculating only the top rlength bytes                   */
  365. /* Note: r will be of length rlength                                    */
  366. /*       2*bnlength <= rlength < bnlength                               */
  367. /*       n1 and n2 can be the same pointer                              */
  368. /* SIDE-EFFECTS: n1 and n2 are changed to their absolute values         */
  369. bn_t unsafe_mult_bn(bn_t r, bn_t n1, bn_t n2)
  370.     {
  371.     int sign1, sign2, samevar;
  372.     int i, j, k, steps, doublesteps, carry_steps, skips;
  373.     bn_t n1p, n2p;      /* pointers for n1, n2 */
  374.     bn_t rp1, rp2, rp3; /* pointers for r */
  375.     U32 prod, sum;
  376.     int bnl; /* temp bnlength holder */
  377.  
  378.     bnl = bnlength;
  379.     if ((sign1 = is_bn_neg(n1)) != 0) /* =, not == */
  380.         neg_a_bn(n1);
  381.     if (!(samevar = (n1 == n2))) /* check to see if they're the same pointer */
  382.         if ((sign2 = is_bn_neg(n2)) != 0) /* =, not == */
  383.             neg_a_bn(n2);
  384.     n1p = n1;
  385.     n2 += (bnlength<<1) - rlength;  /* shift n2 over to where it is needed */
  386.  
  387.     bnlength = rlength;
  388.     clear_bn(r);        /* zero out r, rlength width */
  389.     bnlength = bnl;
  390.  
  391.     steps = (rlength-bnlength)>>1;
  392.     skips = (bnlength>>1) - steps;
  393.     carry_steps = doublesteps = (rlength>>1)-2;
  394.     rp2 = rp1 = r;
  395.     for (i=bnlength>>1; i>0; i--)
  396.         {
  397.         n2p = n2;
  398.         for (j=0; j<steps; j++)
  399.             {
  400.             prod = (U32)access16(n1p) * (U32)access16(n2p); /* U16*U16=U32 */
  401.             sum = (U32)access16(rp2) + prod; /* add to previous, including overflow */
  402.             set16(rp2, (U16)sum); /* save the lower 2 bytes */
  403.             sum >>= 16;             /* keep just the upper 2 bytes */
  404.             rp3 = rp2 + 2;          /* move over 2 bytes */
  405.             sum += access16(rp3);     /* add what was the upper two bytes */
  406.             set16(rp3, (U16)sum); /* save what was the upper two bytes */
  407.             sum >>= 16;             /* keep just the overflow */
  408.             for (k=0; sum != 0 && k<carry_steps; k++)
  409.                 {
  410.                 rp3 += 2;               /* move over 2 bytes */
  411.                 sum += access16(rp3);     /* add to what was the overflow */
  412.                 set16(rp3, (U16)sum); /* save what was the overflow */
  413.                 sum >>= 16;             /* keep just the new overflow */
  414.                 }
  415.             n2p += 2;                   /* increase by two bytes */
  416.             rp2 += 2;
  417.             carry_steps--;
  418.             }
  419.         n1p += 2;   /* increase by two bytes */
  420.  
  421.         if (skips != 0)
  422.             {
  423.             n2 -= 2;    /* shift n2 back a word */
  424.             steps++;    /* one more step this time */
  425.             skips--;    /* keep track of how many times we've done this */
  426.             }
  427.         else
  428.             {
  429.             rp1 += 2;           /* shift forward a word */
  430.             doublesteps--;      /* reduce the carry steps needed next time */
  431.             }
  432.         rp2 = rp1;
  433.         carry_steps = doublesteps;
  434.         }
  435.  
  436.     /* if they were the same or same sign, the product must be positive */
  437.     if (!samevar && sign1 != sign2)
  438.         {
  439.         bnlength = rlength;
  440.         neg_a_bn(r);            /* wider bignumber */
  441.         bnlength = bnl;
  442.         }
  443.     return r;
  444.     }
  445.  
  446. /************************************************************************/
  447. /* r = n^2                                                              */
  448. /*   because of the symetry involved, n^2 is much faster than n*n       */
  449. /*   for a bignumber of length l                                        */
  450. /*      n*n takes l^2 multiplications                                   */
  451. /*      n^2 takes (l^2+l)/2 multiplications                             */
  452. /*          which is about 1/2 n*n as l gets large                      */
  453. /*  uses the fact that (a+b+c+...)^2 = (a^2+b^2+c^2+...)+2(ab+ac+bc+...)*/
  454. /*                                                                      */
  455. /* SIDE-EFFECTS: n is changed to its absolute value                     */
  456. bn_t unsafe_full_square_bn(bn_t r, bn_t n)
  457.     {
  458.     int i, j, k, steps, doublesteps, carry_steps;
  459.     bn_t n1p, n2p;
  460.     bn_t rp1, rp2, rp3;
  461.     U32 prod, sum;
  462.  
  463.     if (is_bn_neg(n))  /* don't need to keep track of sign since the */
  464.         neg_a_bn(n);   /* answer must be positive. */
  465.  
  466.     bnlength <<= 1;
  467.     clear_bn(r);        /* zero out r, double width */
  468.     bnlength >>= 1;
  469.  
  470.     steps = (bnlength>>1)-1;
  471.     carry_steps = doublesteps = (steps<<1) - 1;
  472.     rp2 = rp1 = r + 2;  /* start with second two-byte word */
  473.     n1p = n;
  474.     if (steps != 0) /* if zero, then skip all the middle term calculations */
  475.         {
  476.         for (i=steps; i>0; i--) /* steps gets altered, count backwards */
  477.             {
  478.             n2p = n1p + 2;  /* set n2p pointer to 1 step beyond n1p */
  479.             for (j=0; j<steps; j++)
  480.                 {
  481.                 prod = (U32)access16(n1p) * (U32)access16(n2p); /* U16*U16=U32 */
  482.                 sum = (U32)access16(rp2) + prod; /* add to previous, including overflow */
  483.                 set16(rp2, (U16)sum); /* save the lower 2 bytes */
  484.                 sum >>= 16;             /* keep just the upper 2 bytes */
  485.                 rp3 = rp2 + 2;          /* move over 2 bytes */
  486.                 sum += access16(rp3);     /* add what was the upper two bytes */
  487.                 set16(rp3, (U16)sum); /* save what was the upper two bytes */
  488.                 sum >>= 16;             /* keep just the overflow */
  489.                 for (k=0; sum != 0 && k<carry_steps; k++)
  490.                     {
  491.                     rp3 += 2;               /* move over 2 bytes */
  492.                     sum += access16(rp3);     /* add to what was the overflow */
  493.                     set16(rp3, (U16)sum); /* save what was the overflow */
  494.                     sum >>= 16;             /* keep just the new overflow */
  495.                     }
  496.                 n2p += 2;       /* increase by two bytes */
  497.                 rp2 += 2;
  498.                 carry_steps--;
  499.                 }
  500.             n1p += 2;           /* increase by two bytes */
  501.             rp2 = rp1 += 4;     /* increase by 2 * two bytes */
  502.             carry_steps = doublesteps -= 2;   /* reduce the carry steps needed */
  503.             steps--;
  504.             }
  505.         /* All the middle terms have been multiplied.  Now double it. */
  506.         bnlength <<= 1;     /* double wide bignumber */
  507.         double_a_bn(r);
  508.         bnlength >>= 1;
  509.         /* finished with middle terms */
  510.         }
  511.  
  512.     /* Now go back and add in the squared terms. */
  513.     n1p = n;
  514.     steps = (bnlength>>1);
  515.     carry_steps = doublesteps = (steps<<1) - 2;
  516.     rp1 = r;
  517.     for (i=0; i<steps; i++)
  518.         {
  519.         /* square it */
  520.         prod = (U32)access16(n1p) * (U32)access16(n1p); /* U16*U16=U32 */
  521.         sum = (U32)access16(rp1) + prod; /* add to previous, including overflow */
  522.         set16(rp1, (U16)sum); /* save the lower 2 bytes */
  523.         sum >>= 16;             /* keep just the upper 2 bytes */
  524.         rp3 = rp1 + 2;          /* move over 2 bytes */
  525.         sum += access16(rp3);     /* add what was the upper two bytes */
  526.         set16(rp3, (U16)sum); /* save what was the upper two bytes */
  527.         sum >>= 16;             /* keep just the overflow */
  528.         for (k=0; sum != 0 && k<carry_steps; k++)
  529.             {
  530.             rp3 += 2;               /* move over 2 bytes */
  531.             sum += access16(rp3);     /* add to what was the overflow */
  532.             set16(rp3, (U16)sum); /* save what was the overflow */
  533.             sum >>= 16;             /* keep just the new overflow */
  534.             }
  535.         n1p += 2;       /* increase by 2 bytes */
  536.         rp1 += 4;       /* increase by 4 bytes */
  537.         carry_steps = doublesteps -= 2;
  538.         }
  539.     return r;
  540.     }
  541.  
  542.  
  543. /************************************************************************/
  544. /* r = n^2                                                              */
  545. /*   because of the symetry involved, n^2 is much faster than n*n       */
  546. /*   for a bignumber of length l                                        */
  547. /*      n*n takes l^2 multiplications                                   */
  548. /*      n^2 takes (l^2+l)/2 multiplications                             */
  549. /*          which is about 1/2 n*n as l gets large                      */
  550. /*  uses the fact that (a+b+c+...)^2 = (a^2+b^2+c^2+...)+2(ab+ac+bc+...)*/
  551. /*                                                                      */
  552. /* Note: r will be of length rlength                                    */
  553. /*       2*bnlength >= rlength > bnlength                               */
  554. /* SIDE-EFFECTS: n is changed to its absolute value                     */
  555. bn_t unsafe_square_bn(bn_t r, bn_t n)
  556.     {
  557.     int i, j, k, steps, doublesteps, carry_steps;
  558.     int skips, rodd;
  559.     bn_t n1p, n2p, n3p;
  560.     bn_t rp1, rp2, rp3;
  561.     U32 prod, sum;
  562.     int bnl;
  563.  
  564. /* This whole procedure would be a great deal simpler if we could assume that */
  565. /* rlength < 2*bnlength (that is, not =).  Therefore, we will take the        */
  566. /* easy way out and call full_square_bn() if it is.                           */
  567.     if (rlength == bnlength<<1) /* rlength == 2*bnlength */
  568.         return unsafe_full_square_bn(r, n);    /* call full_square_bn() and quit */
  569.  
  570.     if (is_bn_neg(n))  /* don't need to keep track of sign since the */
  571.         neg_a_bn(n);   /* answer must be positive. */
  572.  
  573.     bnl = bnlength;
  574.     bnlength = rlength;
  575.     clear_bn(r);        /* zero out r, of width rlength */
  576.     bnlength = bnl;
  577.  
  578.     /* determine whether r is on an odd or even two-byte word in the number */
  579.     rodd = (U16)(((bnlength<<1)-rlength)>>1) & 0x0001;
  580.     i = (bnlength>>1)-1;
  581.     steps = (rlength-bnlength)>>1;
  582.     carry_steps = doublesteps = (bnlength>>1)+steps-2;
  583.     skips = (i - steps)>>1;     /* how long to skip over pointer shifts */
  584.     rp2 = rp1 = r;
  585.     n1p = n;
  586.     n3p = n2p = n1p + (((bnlength>>1)-steps)<<1);    /* n2p = n1p + 2*(bnlength/2 - steps) */
  587.     if (i != 0) /* if zero, skip middle term calculations */
  588.         {
  589.         /* i is already set */
  590.         for (; i>0; i--)
  591.             {
  592.             for (j=0; j<steps; j++)
  593.                 {
  594.                 prod = (U32)access16(n1p) * (U32)access16(n2p); /* U16*U16=U32 */
  595.                 sum = (U32)access16(rp2) + prod; /* add to previous, including overflow */
  596.                 set16(rp2, (U16)sum); /* save the lower 2 bytes */
  597.                 sum >>= 16;             /* keep just the upper 2 bytes */
  598.                 rp3 = rp2 + 2;          /* move over 2 bytes */
  599.                 sum += access16(rp3);     /* add what was the upper two bytes */
  600.                 set16(rp3, (U16)sum); /* save what was the upper two bytes */
  601.                 sum >>= 16;             /* keep just the overflow */
  602.                 for (k=0; sum != 0 && k<carry_steps; k++)
  603.                     {
  604.                     rp3 += 2;               /* move over 2 bytes */
  605.                     sum += access16(rp3);     /* add to what was the overflow */
  606.                     set16(rp3, (U16)sum); /* save what was the overflow */
  607.                     sum >>= 16;             /* keep just the new overflow */
  608.                     }
  609.                 n2p += 2;       /* increase by 2-byte word size */
  610.                 rp2 += 2;
  611.                 carry_steps--;
  612.                 }
  613.             n1p += 2;       /* increase by 2-byte word size */
  614.             if (skips > 0)
  615.                 {
  616.                 n2p = n3p -= 2;
  617.                 steps++;
  618.                 skips--;
  619.                 }
  620.             else if (skips == 0)    /* only gets executed once */
  621.                 {
  622.                 steps -= rodd;  /* rodd is 1 or 0 */
  623.                 doublesteps -= rodd+1;
  624.                 rp1 += (rodd+1)<<1;
  625.                 n2p = n1p+2;
  626.                 skips--;
  627.                 }
  628.             else /* skips < 0 */
  629.                 {
  630.                 steps--;
  631.                 doublesteps -= 2;
  632.                 rp1 += 4;           /* add two 2-byte words */
  633.                 n2p = n1p + 2;
  634.                 }
  635.             rp2 = rp1;
  636.             carry_steps = doublesteps;
  637.             }
  638.         /* All the middle terms have been multiplied.  Now double it. */
  639.         bnlength = rlength;
  640.         double_a_bn(r);
  641.         bnlength = bnl;
  642.         }
  643.     /* Now go back and add in the squared terms. */
  644.  
  645.     /* be careful, the next dozen or so lines are confusing!       */
  646.     /* determine whether r is on an odd or even word in the number */
  647.     /* using i as a temporary variable here */
  648.     i = (bnlength<<1)-rlength;
  649.     rp1 = r + ((U16)i & (U16)0x0002);
  650.     i = (U16)((i>>1)+1) & (U16)0xFFFE;
  651.     n1p = n + i;
  652.     /* i here is no longer a temp var., but will be used as a loop counter */
  653.     i = (bnlength - i)>>1;
  654.     carry_steps = doublesteps = (i<<1)-2;
  655.     /* i is already set */
  656.     for (; i>0; i--)
  657.         {
  658.         /* square it */
  659.         prod = (U32)access16(n1p) * (U32)access16(n1p); /* U16*U16=U32 */
  660.         sum = (U32)access16(rp1) + prod; /* add to previous, including overflow */
  661.         set16(rp1, (U16)sum); /* save the lower 2 bytes */
  662.         sum >>= 16;             /* keep just the upper 2 bytes */
  663.         rp3 = rp1 + 2;          /* move over 2 bytes */
  664.         sum += access16(rp3);     /* add what was the upper two bytes */
  665.         set16(rp3, (U16)sum); /* save what was the upper two bytes */
  666.         sum >>= 16;             /* keep just the overflow */
  667.         for (k=0; sum != 0 && k<carry_steps; k++)
  668.             {
  669.             rp3 += 2;               /* move over 2 bytes */
  670.             sum += access16(rp3);     /* add to what was the overflow */
  671.             set16(rp3, (U16)sum); /* save what was the overflow */
  672.             sum >>= 16;             /* keep just the new overflow */
  673.             }
  674.         n1p += 2;
  675.         rp1 += 4;
  676.         carry_steps = doublesteps -= 2;
  677.         }
  678.     return r;
  679.     }
  680.  
  681. /********************************************************************/
  682. /* r = n * u  where u is an unsigned integer */
  683. bn_t mult_bn_int(bn_t r, bn_t n, U16 u)
  684.     {
  685.     int i;
  686.     U32 prod=0;
  687.  
  688.     /* two bytes at a time */
  689.     for (i=0; i<bnlength; i+=2)
  690.         {
  691.         prod += (U32)access16(n+i) * u ; /* n*u */
  692.         set16(r+i, (U16)prod);   /* store the lower 2 bytes */
  693.         prod >>= 16; /* shift the overflow for next time */
  694.         }
  695.     return r;
  696.     }
  697.  
  698. /********************************************************************/
  699. /* r *= u  where u is an unsigned integer */
  700. bn_t mult_a_bn_int(bn_t r, U16 u)
  701.     {
  702.     int i;
  703.     U32 prod=0;
  704.  
  705.     /* two bytes at a time */
  706.     for (i=0; i<bnlength; i+=2)
  707.         {
  708.         prod += (U32)access16(r+i) * u ; /* r*u */
  709.         set16(r+i, (U16)prod);   /* store the lower 2 bytes */
  710.         prod >>= 16; /* shift the overflow for next time */
  711.         }
  712.     return r;
  713.     }
  714.  
  715. /********************************************************************/
  716. /* r = n / u  where u is an unsigned integer */
  717. bn_t unsafe_div_bn_int(bn_t r, bn_t n,  U16 u)
  718.     {
  719.     int i, sign;
  720.     U32 full_number;
  721.     U16 quot, rem=0;
  722.  
  723.     sign = is_bn_neg(n);
  724.     if (sign)
  725.         neg_a_bn(n);
  726.  
  727.     if (u == 0) /* division by zero */
  728.         {
  729.         max_bn(r);
  730.         if (sign)
  731.             neg_a_bn(r);
  732.         return r;
  733.         }
  734.  
  735.     /* two bytes at a time */
  736.     for (i=bnlength-2; i>=0; i-=2)
  737.         {
  738.         full_number = ((U32)rem<<16) + (U32)access16(n+i);
  739.         quot = (U16)(full_number / u);
  740.         rem  = (U16)(full_number % u);
  741.         set16(r+i, quot);
  742.         }
  743.  
  744.     if (sign)
  745.         neg_a_bn(r);
  746.     return r;
  747.     }
  748.  
  749. /********************************************************************/
  750. /* r /= u  where u is an unsigned integer */
  751. bn_t div_a_bn_int(bn_t r, U16 u)
  752.     {
  753.     int i, sign;
  754.     U32 full_number;
  755.     U16 quot, rem=0;
  756.  
  757.     sign = is_bn_neg(r);
  758.     if (sign)
  759.         neg_a_bn(r);
  760.  
  761.     if (u == 0) /* division by zero */
  762.         {
  763.         max_bn(r);
  764.         if (sign)
  765.             neg_a_bn(r);
  766.         return r;
  767.         }
  768.  
  769.     /* two bytes at a time */
  770.     for (i=bnlength-2; i>=0; i-=2)
  771.         {
  772.         full_number = ((U32)rem<<16) + (U32)access16(r+i);
  773.         quot = (U16)(full_number / u);
  774.         rem  = (U16)(full_number % u);
  775.         set16(r+i, quot);
  776.         }
  777.  
  778.     if (sign)
  779.         neg_a_bn(r);
  780.     return r;
  781.     }
  782.  
  783. /*********************************************************************/
  784. /*  f = b                                                            */
  785. /*  Converts a bignumber to a double                                 */
  786. LDBL bntofloat(bn_t n)
  787.     {
  788.     int i;
  789.     int signflag=0;
  790.     int expon;
  791.     bn_t getbyte;
  792.     LDBL f=0;
  793.  
  794.     if (is_bn_neg(n))
  795.         {
  796.         signflag = 1;
  797.         neg_a_bn(n);
  798.         }
  799.  
  800.     expon = intlength - 1;
  801.     getbyte = n + bnlength - 1;
  802.     while (*getbyte == 0 && getbyte >= n)
  803.       {
  804.       getbyte--;
  805.       expon--;
  806.       }
  807.  
  808.     /* There is no need to use all bnlength bytes.  To get the full */
  809.     /* precision of LDBL, all you need is LDBL_MANT_DIG/8+1.        */
  810.     for (i = 0; i < (LDBL_MANT_DIG/8+1) && getbyte >= n; i++, getbyte--)
  811.         {
  812.         f += scale_256(*getbyte,-i);
  813.         }
  814.  
  815.     f = scale_256(f,expon);
  816.  
  817.     if (signflag)
  818.         {
  819.         f = -f;
  820.         neg_a_bn(n);
  821.         }
  822.     return f;
  823.     }
  824.  
  825.