home *** CD-ROM | disk | FTP | other *** search
/ Usenet 1994 October / usenetsourcesnewsgroupsinfomagicoctober1994disk2.iso / unix / volume26 / calc / part11 / zmul.c
C/C++ Source or Header  |  1992-05-09  |  27KB  |  1,040 lines

  1. /*
  2.  * Copyright (c) 1992 David I. Bell
  3.  * Permission is granted to use, distribute, or modify this source,
  4.  * provided that this copyright notice remains intact.
  5.  *
  6.  * Faster than usual multiplying and squaring routines.
  7.  * The algorithm used is the reasonably simple one from Knuth, volume 2,
  8.  * section 4.3.3.  These recursive routines are of speed O(N^1.585)
  9.  * instead of O(N^2).  The usual multiplication and (almost usual) squaring
  10.  * algorithms are used for small numbers.  On a 386 with its compiler, the
  11.  * two algorithms are equal in speed at about 100 decimal digits.
  12.  */
  13.  
  14. #include <stdio.h>
  15. #include "math.h"
  16.  
  17.  
  18. LEN _mul2_ = MUL_ALG2;        /* size of number to use multiply algorithm 2 */
  19. LEN _sq2_ = SQ_ALG2;        /* size of number to use square algorithm 2 */
  20.  
  21.  
  22. #if 0
  23. static char *abortmsg = "Calculation aborted";
  24. static char *memmsg = "Not enough memory";
  25. #endif
  26.  
  27. static HALF *tempbuf;        /* temporary buffer for multiply and square */
  28.  
  29. static LEN domul proto((HALF *v1, LEN size1, HALF *v2, LEN size2, HALF *ans));
  30. static LEN dosquare proto((HALF *vp, LEN size, HALF *ans));
  31.  
  32.  
  33. /*
  34.  * Multiply two numbers using the following formula recursively:
  35.  *    (A*S+B)*(C*S+D) = (S^2+S)*A*C + S*(A-B)*(D-C) + (S+1)*B*D
  36.  * where S is a power of 2^16, and so multiplies by it are shifts, and
  37.  * A,B,C,D are the left and right halfs of the numbers to be multiplied.
  38.  */
  39. void
  40. zmul(z1, z2, res)
  41.     ZVALUE z1, z2;        /* numbers to multiply */
  42.     ZVALUE *res;        /* result of multiplication */
  43. {
  44.     LEN len;        /* size of array */
  45.  
  46.     if (iszero(z1) || iszero(z2)) {
  47.         *res = _zero_;
  48.         return;
  49.     }
  50.     if (isunit(z1)) {
  51.         zcopy(z2, res);
  52.         res->sign = (z1.sign != z2.sign);
  53.         return;
  54.     }
  55.     if (isunit(z2)) {
  56.         zcopy(z1, res);
  57.         res->sign = (z1.sign != z2.sign);
  58.         return;
  59.     }
  60.  
  61.     /*
  62.      * Allocate a temporary buffer for the recursion levels to use.
  63.      * An array needs to be allocated large enough for all of the
  64.      * temporary results to fit in.  This size is about twice the size
  65.      * of the largest original number, since each recursion level uses
  66.      * the size of its given number, and whose size is 1/2 the size of
  67.      * the previous level.  The sum of the infinite series is 2.
  68.      * Add some extra words because of rounding when dividing by 2
  69.      * and also because of the extra word that each multiply needs.
  70.      */
  71.     len = z1.len;
  72.     if (len < z2.len)
  73.         len = z2.len;
  74.     len = 2 * len + 64;
  75.     tempbuf = zalloctemp(len);
  76.  
  77.     res->sign = (z1.sign != z2.sign);
  78.     res->v = alloc(z1.len + z2.len + 1);
  79.     res->len = domul(z1.v, z1.len, z2.v, z2.len, res->v);
  80. }
  81.  
  82.  
  83. /*
  84.  * Recursive routine to multiply two numbers by splitting them up into
  85.  * two numbers of half the size, and using the results of multiplying the
  86.  * subpieces.  The result is placed in the indicated array, which must be
  87.  * large enough for the result plus one extra word (size1 + size2 + 1).
  88.  * Returns the actual size of the result with leading zeroes stripped.
  89.  * This also uses a temporary array which must be twice as large as
  90.  * one more than the size of the number at the top level recursive call.
  91.  */
  92. static LEN
  93. domul(v1, size1, v2, size2, ans)
  94.     HALF *v1;        /* first number */
  95.     LEN size1;        /* size of first number */
  96.     HALF *v2;        /* second number */
  97.     LEN size2;        /* size of second number */
  98.     HALF *ans;        /* location for result */
  99. {
  100.     LEN shift;        /* amount numbers are shifted by */
  101.     LEN sizeA;        /* size of left half of first number */
  102.     LEN sizeB;        /* size of right half of first number */
  103.     LEN sizeC;        /* size of left half of second number */
  104.     LEN sizeD;        /* size of right half of second number */
  105.     LEN sizeAB;        /* size of subtraction of A and B */
  106.     LEN sizeDC;        /* size of subtraction of D and C */
  107.     LEN sizeABDC;        /* size of product of above two results */
  108.     LEN subsize;        /* size of difference of halfs */
  109.     LEN copysize;        /* size of number left to copy */
  110.     LEN sizetotal;        /* total size of product */
  111.     LEN len;        /* temporary length */
  112.     HALF *baseA;        /* base of left half of first number */
  113.     HALF *baseB;        /* base of right half of first number */
  114.     HALF *baseC;        /* base of left half of second number */
  115.     HALF *baseD;        /* base of right half of second number */
  116.     HALF *baseAB;        /* base of result of subtraction of A and B */
  117.     HALF *baseDC;        /* base of result of subtraction of D and C */
  118.     HALF *baseABDC;        /* base of product of above two results */
  119.     HALF *baseAC;        /* base of product of A and C */
  120.     HALF *baseBD;        /* base of product of B and D */
  121.     FULL carry;        /* carry digit for small multiplications */
  122.     FULL carryACBD;        /* carry from addition of A*C and B*D */
  123.     FULL digit;        /* single digit multiplying by */
  124.     HALF *temp;        /* base for temporary calculations */
  125.     BOOL neg;        /* whether imtermediate term is negative */
  126.     register HALF *hd, *h1, *h2;    /* for inner loops */
  127.     SIUNION sival;        /* for addition of digits */
  128.  
  129.     /*
  130.      * Trim the numbers of leading zeroes and initialize the
  131.      * estimated size of the result.
  132.      */
  133.     hd = &v1[size1 - 1];
  134.     while ((*hd == 0) && (size1 > 1)) {
  135.         hd--;
  136.         size1--;
  137.     }
  138.     hd = &v2[size2 - 1];
  139.     while ((*hd == 0) && (size2 > 1)) {
  140.         hd--;
  141.         size2--;
  142.     }
  143.     sizetotal = size1 + size2;
  144.  
  145.     /*
  146.      * First check for zero answer.
  147.      */
  148.     if (((size1 == 1) && (*v1 == 0)) || ((size2 == 1) && (*v2 == 0))) {
  149.         *ans = 0;
  150.         return 1;
  151.     }
  152.  
  153.     /*
  154.      * Exchange the two numbers if necessary to make the number of
  155.      * digits of the first number be greater than or equal to the
  156.      * second number.
  157.      */
  158.     if (size1 < size2) {
  159.         len = size1; size1 = size2; size2 = len;
  160.         hd = v1; v1 = v2; v2 = hd;
  161.     }
  162.  
  163.     /*
  164.      * If the smaller number has only a few digits, then calculate
  165.      * the result in the normal manner in order to avoid the overhead
  166.      * of the recursion for small numbers.  The number of digits where
  167.      * the algorithm changes is settable from 2 to maxint.
  168.      */
  169.     if (size2 < _mul2_) {
  170.         /*
  171.          * First clear the top part of the result, and then multiply
  172.          * by the lowest digit to get the first partial sum.  Later
  173.          * products will then add into this result.
  174.          */
  175.         hd = &ans[size1];
  176.         len = size2;
  177.         while (len--)
  178.             *hd++ = 0;
  179.  
  180.         digit = *v2++;
  181.         h1 = v1;
  182.         hd = ans;
  183.         carry = 0;
  184.         len = size1;
  185.         while (len >= 4) {    /* expand the loop some */
  186.             len -= 4;
  187.             sival.ivalue = ((FULL) *h1++) * digit + carry;
  188.             *hd++ = sival.silow;
  189.             carry = sival.sihigh;
  190.             sival.ivalue = ((FULL) *h1++) * digit + carry;
  191.             *hd++ = sival.silow;
  192.             carry = sival.sihigh;
  193.             sival.ivalue = ((FULL) *h1++) * digit + carry;
  194.             *hd++ = sival.silow;
  195.             carry = sival.sihigh;
  196.             sival.ivalue = ((FULL) *h1++) * digit + carry;
  197.             *hd++ = sival.silow;
  198.             carry = sival.sihigh;
  199.         }
  200.         while (len--) {
  201.             sival.ivalue = ((FULL) *h1++) * digit + carry;
  202.             *hd++ = sival.silow;
  203.             carry = sival.sihigh;
  204.         }
  205.         *hd = (HALF)carry;
  206.  
  207.         /*
  208.          * Now multiply by the remaining digits of the second number,
  209.          * adding each product into the final result.
  210.          */
  211.         h2 = ans;
  212.         while (--size2 > 0) {
  213.             digit = *v2++;
  214.             h1 = v1;
  215.             hd = ++h2;
  216.             if (digit == 0)
  217.                 continue;
  218.             carry = 0;
  219.             len = size1;
  220.             while (len >= 4) {    /* expand the loop some */
  221.                 len -= 4;
  222.                 sival.ivalue = ((FULL) *h1++) * digit
  223.                     + ((FULL) *hd) + carry;
  224.                 *hd++ = sival.silow;
  225.                 carry = sival.sihigh;
  226.                 sival.ivalue = ((FULL) *h1++) * digit
  227.                     + ((FULL) *hd) + carry;
  228.                 *hd++ = sival.silow;
  229.                 carry = sival.sihigh;
  230.                 sival.ivalue = ((FULL) *h1++) * digit
  231.                     + ((FULL) *hd) + carry;
  232.                 *hd++ = sival.silow;
  233.                 carry = sival.sihigh;
  234.                 sival.ivalue = ((FULL) *h1++) * digit
  235.                     + ((FULL) *hd) + carry;
  236.                 *hd++ = sival.silow;
  237.                 carry = sival.sihigh;
  238.             }
  239.             while (len--) {
  240.                 sival.ivalue = ((FULL) *h1++) * digit
  241.                     + ((FULL) *hd) + carry;
  242.                 *hd++ = sival.silow;
  243.                 carry = sival.sihigh;
  244.             }
  245.             while (carry) {
  246.                 sival.ivalue = ((FULL) *hd) + carry;
  247.                 *hd++ = sival.silow;
  248.                 carry = sival.sihigh;
  249.             }
  250.         }
  251.  
  252.         /*
  253.          * Now return the true size of the number.
  254.          */
  255.         len = sizetotal;
  256.         hd = &ans[len - 1];
  257.         while ((*hd == 0) && (len > 1)) {
  258.             hd--;
  259.             len--;
  260.         }
  261.         return len;
  262.     }
  263.  
  264.     /*
  265.      * Need to multiply by a large number.
  266.      * Allocate temporary space for calculations, and calculate the
  267.      * value for the shift.  The shift value is 1/2 the size of the
  268.      * larger (first) number (rounded up).  The amount of temporary
  269.      * space needed is twice the size of the shift, plus one more word
  270.      * for the multiply to use.
  271.      */
  272.     shift = (size1 + 1) / 2;
  273.     temp = tempbuf;
  274.     tempbuf += (2 * shift) + 1;
  275.  
  276.     /*
  277.      * Determine the sizes and locations of all the numbers.
  278.      * The value of sizeC can be negative, and this is checked later.
  279.      * The value of sizeD is limited by the full size of the number.
  280.      */
  281.     baseA = v1 + shift;
  282.     baseB = v1;
  283.     baseC = v2 + shift;
  284.     baseD = v2;
  285.     baseAB = ans;
  286.     baseDC = ans + shift;
  287.     baseAC = ans + shift * 2;
  288.     baseBD = ans;
  289.  
  290.     sizeA = size1 - shift;
  291.     sizeC = size2 - shift;
  292.  
  293.     sizeB = shift;
  294.     hd = &baseB[shift - 1];
  295.     while ((*hd == 0) && (sizeB > 1)) {
  296.         hd--;
  297.         sizeB--;
  298.     }
  299.  
  300.     sizeD = shift;
  301.     if (sizeD > size2)
  302.         sizeD = size2;
  303.     hd = &baseD[sizeD - 1];
  304.     while ((*hd == 0) && (sizeD > 1)) {
  305.         hd--;
  306.         sizeD--;
  307.     }
  308.  
  309.     /*
  310.      * If the smaller number has a high half of zero, then calculate
  311.      * the result by breaking up the first number into two numbers
  312.      * and combining the results using the obvious formula:
  313.      *    (A*S+B) * D = (A*D)*S + B*D
  314.      */
  315.     if (sizeC <= 0) {
  316.         len = domul(baseB, sizeB, baseD, sizeD, ans);
  317.         hd = &ans[len];
  318.         len = sizetotal - len;
  319.         while (len--)
  320.             *hd++ = 0;
  321.  
  322.         /*
  323.          * Add the second number into the first number, shifted
  324.          * over at the correct position.
  325.          */
  326.         len = domul(baseA, sizeA, baseD, sizeD, temp);
  327.         h1 = temp;
  328.         hd = ans + shift;
  329.         carry = 0;
  330.         while (len--) {
  331.             sival.ivalue = ((FULL) *h1++) + ((FULL) *hd) + carry;
  332.             *hd++ = sival.silow;
  333.             carry = sival.sihigh;
  334.         }
  335.         while (carry) {
  336.             sival.ivalue = ((FULL) *hd) + carry;
  337.             *hd++ = sival.silow;
  338.             carry = sival.sihigh;
  339.         }
  340.  
  341.         /*
  342.          * Determine the final size of the number and return it.
  343.          */
  344.         len = sizetotal;
  345.         hd = &ans[len - 1];
  346.         while ((*hd == 0) && (len > 1)) {
  347.             hd--;
  348.             len--;
  349.         }
  350.         tempbuf = temp;
  351.         return len;
  352.     }
  353.  
  354.     /*
  355.      * Now we know that the high halfs of the numbers are nonzero,
  356.      * so we can use the complete formula.
  357.      *    (A*S+B)*(C*S+D) = (S^2+S)*A*C + S*(A-B)*(D-C) + (S+1)*B*D.
  358.      * The steps are done in the following order:
  359.      *    A-B
  360.      *    D-C
  361.      *    (A-B)*(D-C)
  362.      *    S^2*A*C + B*D
  363.      *    (S^2+S)*A*C + (S+1)*B*D                (*)
  364.      *    (S^2+S)*A*C + S*(A-B)*(D-C) + (S+1)*B*D
  365.      *
  366.      * Note: step (*) above can produce a result which is larger than
  367.      * the final product will be, and this is where the extra word
  368.      * needed in the product comes from.  After the final subtraction is
  369.      * done, the result fits in the expected size.  Using the extra word
  370.      * is easier than suppressing the carries and borrows everywhere.
  371.      *
  372.      * Begin by forming the product (A-B)*(D-C) into a temporary
  373.      * location that we save until the final step.  Do each subtraction
  374.      * at positions 0 and S.  Be very careful about the relative sizes
  375.      * of the numbers since this result can be negative.  For the first
  376.      * step calculate the absolute difference of A and B into a temporary
  377.      * location at position 0 of the result.  Negate the sign if A is
  378.      * smaller than B.
  379.      */
  380.     neg = FALSE;
  381.     if (sizeA == sizeB) {
  382.         len = sizeA;
  383.         h1 = &baseA[len - 1];
  384.         h2 = &baseB[len - 1];
  385.         while ((len > 1) && (*h1 == *h2)) {
  386.             len--;
  387.             h1--;
  388.             h2--;
  389.         }
  390.     }
  391.     if ((sizeA > sizeB) || ((sizeA == sizeB) && (*h1 > *h2)))
  392.     {
  393.         h1 = baseA;
  394.         h2 = baseB;
  395.         sizeAB = sizeA;
  396.         subsize = sizeB;
  397.     } else {
  398.         neg = !neg;
  399.         h1 = baseB;
  400.         h2 = baseA;
  401.         sizeAB = sizeB;
  402.         subsize = sizeA;
  403.     }
  404.     copysize = sizeAB - subsize;
  405.  
  406.     hd = baseAB;
  407.     carry = 0;
  408.     while (subsize--) {
  409.         sival.ivalue = BASE1 - ((FULL) *h1++) + ((FULL) *h2++) + carry;
  410.         *hd++ = BASE1 - sival.silow;
  411.         carry = sival.sihigh;
  412.     }
  413.     while (copysize--) {
  414.         sival.ivalue = (BASE1 - ((FULL) *h1++)) + carry;
  415.         *hd++ = BASE1 - sival.silow;
  416.         carry = sival.sihigh;
  417.     }
  418.  
  419.     hd = &baseAB[sizeAB - 1];
  420.     while ((*hd == 0) && (sizeAB > 1)) {
  421.         hd--;
  422.         sizeAB--;
  423.     }
  424.  
  425.     /*
  426.      * This completes the calculation of abs(A-B).  For the next step
  427.      * calculate the absolute difference of D and C into a temporary
  428.      * location at position S of the result.  Negate the sign if C is
  429.      * larger than D.
  430.      */
  431.     if (sizeC == sizeD) {
  432.         len = sizeC;
  433.         h1 = &baseC[len - 1];
  434.         h2 = &baseD[len - 1];
  435.         while ((len > 1) && (*h1 == *h2)) {
  436.             len--;
  437.             h1--;
  438.             h2--;
  439.         }
  440.     }
  441.     if ((sizeC > sizeD) || ((sizeC == sizeD) && (*h1 > *h2)))
  442.     {
  443.         neg = !neg;
  444.         h1 = baseC;
  445.         h2 = baseD;
  446.         sizeDC = sizeC;
  447.         subsize = sizeD;
  448.     } else {
  449.         h1 = baseD;
  450.         h2 = baseC;
  451.         sizeDC = sizeD;
  452.         subsize = sizeC;
  453.     }
  454.     copysize = sizeDC - subsize;
  455.  
  456.     hd = baseDC;
  457.     carry = 0;
  458.     while (subsize--) {
  459.         sival.ivalue = BASE1 - ((FULL) *h1++) + ((FULL) *h2++) + carry;
  460.         *hd++ = BASE1 - sival.silow;
  461.         carry = sival.sihigh;
  462.     }
  463.     while (copysize--) {
  464.         sival.ivalue = (BASE1 - ((FULL) *h1++)) + carry;
  465.         *hd++ = BASE1 - sival.silow;
  466.         carry = sival.sihigh;
  467.     }
  468.     hd = &baseDC[sizeDC - 1];
  469.     while ((*hd == 0) && (sizeDC > 1)) {
  470.         hd--;
  471.         sizeDC--;
  472.     }
  473.  
  474.     /*
  475.      * This completes the calculation of abs(D-C).  Now multiply
  476.      * together abs(A-B) and abs(D-C) into a temporary location,
  477.      * which is preserved until the final steps.
  478.      */
  479.     baseABDC = temp;
  480.     sizeABDC = domul(baseAB, sizeAB, baseDC, sizeDC, baseABDC);
  481.  
  482.     /*
  483.      * Now calculate B*D and A*C into one of their two final locations.
  484.      * Make sure the high order digits of the products are zeroed since
  485.      * this initializes the final result.  Be careful about this zeroing
  486.      * since the size of the high order words might be smaller than
  487.      * the shift size.  Do B*D first since the multiplies use one more
  488.      * word than the size of the product.  Also zero the final extra
  489.      * word in the result for possible carries to use.
  490.      */
  491.     len = domul(baseB, sizeB, baseD, sizeD, baseBD);
  492.     hd = &baseBD[len];
  493.     len = shift * 2 - len;
  494.     while (len--)
  495.         *hd++ = 0;
  496.  
  497.     len = domul(baseA, sizeA, baseC, sizeC, baseAC);
  498.     hd = &baseAC[len];
  499.     len = sizetotal - shift * 2 - len + 1;
  500.     while (len--)
  501.         *hd++ = 0;
  502.  
  503.     /*
  504.      * Now add in A*C and B*D into themselves at the other shifted
  505.      * position that they need.  This addition is tricky in order to
  506.      * make sure that the two additions cannot interfere with each other.
  507.      * Therefore we first add in the top half of B*D and the lower half
  508.      * of A*C.  The sources and destinations of these two additions
  509.      * overlap, and so the same answer results from the two additions,
  510.      * thus only two pointers suffice for both additions.  Keep the
  511.      * final carry from these additions for later use since we cannot
  512.      * afford to change the top half of A*C yet.
  513.      */
  514.     h1 = baseBD + shift;
  515.     h2 = baseAC;
  516.     carryACBD = 0;
  517.     len = shift;
  518.     while (len--) {
  519.         sival.ivalue = ((FULL) *h1) + ((FULL) *h2) + carryACBD;
  520.         *h1++ = sival.silow;
  521.         *h2++ = sival.silow;
  522.         carryACBD = sival.sihigh;
  523.     }
  524.  
  525.     /*
  526.      * Now add in the bottom half of B*D and the top half of A*C.
  527.      * These additions are straightforward, except that A*C should
  528.      * be done first because of possible carries from B*D, and the
  529.      * top half of A*C might not exist.  Add in one of the carries
  530.      * from the previous addition while we are at it.
  531.      */
  532.     h1 = baseAC + shift;
  533.     hd = baseAC;
  534.     carry = carryACBD;
  535.     len = sizetotal - 3 * shift;
  536.     while (len--) {
  537.         sival.ivalue = ((FULL) *h1++) + ((FULL) *hd) + carry;
  538.         *hd++ = sival.silow;
  539.         carry = sival.sihigh;
  540.     }
  541.     while (carry) {
  542.         sival.ivalue = ((FULL) *hd) + carry;
  543.         *hd++ = sival.silow;
  544.         carry = sival.sihigh;
  545.     }
  546.  
  547.     h1 = baseBD;
  548.     hd = baseBD + shift;
  549.     carry = 0;
  550.     len = shift;
  551.     while (len--) {
  552.         sival.ivalue = ((FULL) *h1++) + ((FULL) *hd) + carry;
  553.         *hd++ = sival.silow;
  554.         carry = sival.sihigh;
  555.     }
  556.     while (carry) {
  557.         sival.ivalue = ((FULL) *hd) + carry;
  558.         *hd++ = sival.silow;
  559.         carry = sival.sihigh;
  560.     }
  561.  
  562.     /*
  563.      * Now finally add in the other delayed carry from the
  564.      * above addition.
  565.      */
  566.     hd = baseAC + shift;
  567.     while (carryACBD) {
  568.         sival.ivalue = ((FULL) *hd) + carryACBD;
  569.         *hd++ = sival.silow;
  570.         carryACBD = sival.sihigh;
  571.     }
  572.  
  573.     /*
  574.      * Now finally add or subtract (A-B)*(D-C) into the final result at
  575.      * the correct position (S), according to whether it is positive or
  576.      * negative.  When subtracting, the answer cannot go negative.
  577.      */
  578.     h1 = baseABDC;
  579.     hd = ans + shift;
  580.     carry = 0;
  581.     len = sizeABDC;
  582.     if (neg) {
  583.         while (len--) {
  584.             sival.ivalue = BASE1 - ((FULL) *hd) +
  585.                 ((FULL) *h1++) + carry;
  586.             *hd++ = BASE1 - sival.silow;
  587.             carry = sival.sihigh;
  588.         }
  589.         while (carry) {
  590.             sival.ivalue = BASE1 - ((FULL) *hd) + carry;
  591.             *hd++ = BASE1 - sival.silow;
  592.             carry = sival.sihigh;
  593.         }
  594.     } else {
  595.         while (len--) {
  596.             sival.ivalue = ((FULL) *h1++) + ((FULL) *hd) + carry;
  597.             *hd++ = sival.silow;
  598.             carry = sival.sihigh;
  599.         }
  600.         while (carry) {
  601.             sival.ivalue = ((FULL) *hd) + carry;
  602.             *hd++ = sival.silow;
  603.             carry = sival.sihigh;
  604.         }
  605.     }
  606.  
  607.     /*
  608.      * Finally determine the size of the final result and return that.
  609.      */
  610.     len = sizetotal;
  611.     hd = &ans[len - 1];
  612.     while ((*hd == 0) && (len > 1)) {
  613.         hd--;
  614.         len--;
  615.     }
  616.     tempbuf = temp;
  617.     return len;
  618. }
  619.  
  620.  
  621. /*
  622.  * Square a number by using the following formula recursively:
  623.  *    (A*S+B)^2 = (S^2+S)*A^2 + (S+1)*B^2 - S*(A-B)^2
  624.  * where S is a power of 2^16, and so multiplies by it are shifts,
  625.  * and A and B are the left and right halfs of the number to square.
  626.  */
  627. void
  628. zsquare(z, res)
  629.     ZVALUE z, *res;
  630. {
  631.     LEN len;
  632.  
  633.     if (iszero(z)) {
  634.         *res = _zero_;
  635.         return;
  636.     }
  637.     if (isunit(z)) {
  638.         *res = _one_;
  639.         return;
  640.     }
  641.  
  642.     /*
  643.      * Allocate a temporary array if necessary for the recursion to use.
  644.      * The array needs to be allocated large enough for all of the
  645.      * temporary results to fit in.  This size is about 3 times the
  646.      * size of the original number, since each recursion level uses 3/2
  647.      * of the size of its given number, and whose size is 1/2 the size
  648.      * of the previous level.  The sum of the infinite series is 3.
  649.      * Allocate some extra words for rounding up the sizes.
  650.      */
  651.     len = 3 * z.len + 32;
  652.     tempbuf = zalloctemp(len);
  653.  
  654.     res->sign = 0;
  655.     res->v = alloc(z.len * 2);
  656.     res->len = dosquare(z.v, z.len, res->v);
  657. }
  658.  
  659.  
  660. /*
  661.  * Recursive routine to square a number by splitting it up into two numbers
  662.  * of half the size, and using the results of squaring the subpieces.
  663.  * The result is placed in the indicated array, which must be large
  664.  * enough for the result (size * 2).  Returns the size of the result.
  665.  * This uses a temporary array which must be 3 times as large as the
  666.  * size of the number at the top level recursive call.
  667.  */
  668. static LEN
  669. dosquare(vp, size, ans)
  670.     HALF *vp;        /* value to be squared */
  671.     LEN size;        /* length of value to square */
  672.     HALF *ans;        /* location for result */
  673. {
  674.     LEN shift;        /* amount numbers are shifted by */
  675.     LEN sizeA;        /* size of left half of number to square */
  676.     LEN sizeB;        /* size of right half of number to square */
  677.     LEN sizeAA;        /* size of square of left half */
  678.     LEN sizeBB;        /* size of square of right half */
  679.     LEN sizeAABB;        /* size of sum of squares of A and B */
  680.     LEN sizeAB;        /* size of difference of A and B */
  681.     LEN sizeABAB;        /* size of square of difference of A and B */
  682.     LEN subsize;        /* size of difference of halfs */
  683.     LEN copysize;        /* size of number left to copy */
  684.     LEN sumsize;        /* size of sum */
  685.     LEN sizetotal;        /* total size of square */
  686.     LEN len;        /* temporary length */
  687.     LEN len1;        /* another temporary length */
  688.     FULL carry;        /* carry digit for small multiplications */
  689.     FULL digit;        /* single digit multiplying by */
  690.     HALF *temp;        /* base for temporary calculations */
  691.     HALF *baseA;        /* base of left half of number */
  692.     HALF *baseB;        /* base of right half of number */
  693.     HALF *baseAA;        /* base of square of left half of number */
  694.     HALF *baseBB;        /* base of square of right half of number */
  695.     HALF *baseAABB;        /* base of sum of squares of A and B */
  696.     HALF *baseAB;        /* base of difference of A and B */
  697.     HALF *baseABAB;        /* base of square of difference of A and B */
  698.     register HALF *hd, *h1, *h2, *h3;    /* for inner loops */
  699.     SIUNION sival;        /* for addition of digits */
  700.  
  701.     /*
  702.      * First trim the number of leading zeroes.
  703.      */
  704.     hd = &vp[size - 1];
  705.     while ((*hd == 0) && (size > 1)) {
  706.         size--;
  707.         hd--;
  708.     }
  709.     sizetotal = size + size;
  710.  
  711.     /*
  712.      * If the number has only a small number of digits, then use the
  713.      * (almost) normal multiplication method.  Multiply each halfword
  714.      * only by those halfwards further on in the number.  Missed terms
  715.      * will then be the same pairs of products repeated, and the squares
  716.      * of each halfword.  The first case is handled by doubling the
  717.      * result.  The second case is handled explicitly.  The number of
  718.      * digits where the algorithm changes is settable from 2 to maxint.
  719.      */
  720.     if (size < _sq2_) {
  721.         hd = ans;
  722.         len = sizetotal;
  723.         while (len--)
  724.             *hd++ = 0;
  725.  
  726.         h2 = vp;
  727.         hd = ans + 1;
  728.         for (len = size; len--; hd += 2) {
  729.             digit = (FULL) *h2++;
  730.             if (digit == 0)
  731.                 continue;
  732.             h3 = h2;
  733.             h1 = hd;
  734.             carry = 0;
  735.             len1 = len;
  736.             while (len1 >= 4) {    /* expand the loop some */
  737.                 len1 -= 4;
  738.                 sival.ivalue = (digit * ((FULL) *h3++))
  739.                     + ((FULL) *h1) + carry;
  740.                 *h1++ = sival.silow;
  741.                 sival.ivalue = (digit * ((FULL) *h3++))
  742.                     + ((FULL) *h1) + ((FULL) sival.sihigh);
  743.                 *h1++ = sival.silow;
  744.                 sival.ivalue = (digit * ((FULL) *h3++))
  745.                     + ((FULL) *h1) + ((FULL) sival.sihigh);
  746.                 *h1++ = sival.silow;
  747.                 sival.ivalue = (digit * ((FULL) *h3++))
  748.                     + ((FULL) *h1) + ((FULL) sival.sihigh);
  749.                 *h1++ = sival.silow;
  750.                 carry = sival.sihigh;
  751.             }
  752.             while (len1--) {
  753.                 sival.ivalue = (digit * ((FULL) *h3++))
  754.                     + ((FULL) *h1) + carry;
  755.                 *h1++ = sival.silow;
  756.                 carry = sival.sihigh;
  757.             }
  758.             while (carry) {
  759.                 sival.ivalue = ((FULL) *h1) + carry;
  760.                 *h1++ = sival.silow;
  761.                 carry = sival.sihigh;
  762.             }
  763.         }
  764.  
  765.         /*
  766.          * Now double the result.
  767.          * There is no final carry to worry about because we
  768.          * handle all digits of the result which must fit.
  769.          */
  770.         carry = 0;
  771.         hd = ans;
  772.         len = sizetotal;
  773.         while (len--) {
  774.             digit = ((FULL) *hd);
  775.             sival.ivalue = digit + digit + carry;
  776.             *hd++ = sival.silow;
  777.             carry = sival.sihigh;
  778.         }
  779.  
  780.         /*
  781.          * Now add in the squares of each halfword.
  782.          */
  783.         carry = 0;
  784.         hd = ans;
  785.         h3 = vp;
  786.         len = size;
  787.         while (len--) {
  788.             digit = ((FULL) *h3++);
  789.             sival.ivalue = digit * digit + ((FULL) *hd) + carry;
  790.             *hd++ = sival.silow;
  791.             carry = sival.sihigh;
  792.             sival.ivalue = ((FULL) *hd) + carry;
  793.             *hd++ = sival.silow;
  794.             carry = sival.sihigh;
  795.         }
  796.         while (carry) {
  797.             sival.ivalue = ((FULL) *hd) + carry;
  798.             *hd++ = sival.silow;
  799.             carry = sival.sihigh;
  800.         }
  801.  
  802.         /*
  803.          * Finally return the size of the result.
  804.          */
  805.         len = sizetotal;
  806.         hd = &ans[len - 1];
  807.         while ((*hd == 0) && (len > 1)) {
  808.             len--;
  809.             hd--;
  810.         }
  811.         return len;
  812.     }
  813.  
  814.     /*
  815.      * The number to be squared is large.
  816.      * Allocate temporary space and determine the sizes and
  817.      * positions of the values to be calculated.
  818.      */
  819.     temp = tempbuf;
  820.     tempbuf += (3 * (size + 1) / 2);
  821.  
  822.     sizeA = size / 2;
  823.     sizeB = size - sizeA;
  824.     shift = sizeB;
  825.     baseA = vp + sizeB;
  826.     baseB = vp;
  827.     baseAA = &ans[shift * 2];
  828.     baseBB = ans;
  829.     baseAABB = temp;
  830.     baseAB = temp;
  831.     baseABAB = &temp[shift];
  832.  
  833.     /*
  834.      * Trim the second number of leading zeroes.
  835.      */
  836.     hd = &baseB[sizeB - 1];
  837.     while ((*hd == 0) && (sizeB > 1)) {
  838.         sizeB--;
  839.         hd--;
  840.     }
  841.  
  842.     /*
  843.      * Now to proceed to calculate the result using the formula.
  844.      *    (A*S+B)^2 = (S^2+S)*A^2 + (S+1)*B^2 - S*(A-B)^2.
  845.      * The steps are done in the following order:
  846.      *    S^2*A^2 + B^2
  847.      *    A^2 + B^2
  848.      *    (S^2+S)*A^2 + (S+1)*B^2
  849.      *    (A-B)^2
  850.      *    (S^2+S)*A^2 + (S+1)*B^2 - S*(A-B)^2.
  851.      *
  852.      * Begin by forming the squares of two the halfs concatenated
  853.      * together in the final result location.  Make sure that the
  854.      * highest words of the results are zero.
  855.      */
  856.     sizeBB = dosquare(baseB, sizeB, baseBB);
  857.     hd = &baseBB[sizeBB];
  858.     len = shift * 2 - sizeBB;
  859.     while (len--)
  860.         *hd++ = 0;
  861.  
  862.     sizeAA = dosquare(baseA, sizeA, baseAA);
  863.     hd = &baseAA[sizeAA];
  864.     len = sizetotal - shift * 2 - sizeAA;
  865.     while (len--)
  866.         *hd++ = 0;
  867.  
  868.     /*
  869.      * Sum the two squares into a temporary location.
  870.      */
  871.     if (sizeAA >= sizeBB) {
  872.         h1 = baseAA;
  873.         h2 = baseBB;
  874.         sizeAABB = sizeAA;
  875.         sumsize = sizeBB;
  876.     } else {
  877.         h1 = baseBB;
  878.         h2 = baseAA;
  879.         sizeAABB = sizeBB;
  880.         sumsize = sizeAA;
  881.     }
  882.     copysize = sizeAABB - sumsize;
  883.  
  884.     hd = baseAABB;
  885.     carry = 0;
  886.     while (sumsize--) {
  887.         sival.ivalue = ((FULL) *h1++) + ((FULL) *h2++) + carry;
  888.         *hd++ = sival.silow;
  889.         carry = sival.sihigh;
  890.     }
  891.     while (copysize--) {
  892.         sival.ivalue = ((FULL) *h1++) + carry;
  893.         *hd++ = sival.silow;
  894.         carry = sival.sihigh;
  895.     }
  896.     if (carry) {
  897.         *hd = (HALF)carry;
  898.         sizeAABB++;
  899.     }
  900.  
  901.     /*
  902.      * Add the sum back into the previously calculated squares
  903.      * shifted over to the proper location.
  904.      */
  905.     h1 = baseAABB;
  906.     hd = ans + shift;
  907.     carry = 0;
  908.     len = sizeAABB;
  909.     while (len--) {
  910.         sival.ivalue = ((FULL) *hd) + ((FULL) *h1++) + carry;
  911.         *hd++ = sival.silow;
  912.         carry = sival.sihigh;
  913.     }
  914.     while (carry) {
  915.         sival.ivalue = ((FULL) *hd) + carry;
  916.         *hd++ = sival.silow;
  917.         carry = sival.sihigh;
  918.     }
  919.  
  920.     /*
  921.      * Calculate the absolute value of the difference of the two halfs
  922.      * into a temporary location.
  923.      */
  924.     if (sizeA == sizeB) {
  925.         len = sizeA;
  926.         h1 = &baseA[len - 1];
  927.         h2 = &baseB[len - 1];
  928.         while ((len > 1) && (*h1 == *h2)) {
  929.             len--;
  930.             h1--;
  931.             h2--;
  932.         }
  933.     }
  934.     if ((sizeA > sizeB) || ((sizeA == sizeB) && (*h1 > *h2)))
  935.     {
  936.         h1 = baseA;
  937.         h2 = baseB;
  938.         sizeAB = sizeA;
  939.         subsize = sizeB;
  940.     } else {
  941.         h1 = baseB;
  942.         h2 = baseA;
  943.         sizeAB = sizeB;
  944.         subsize = sizeA;
  945.     }
  946.     copysize = sizeAB - subsize;
  947.  
  948.     hd = baseAB;
  949.     carry = 0;
  950.     while (subsize--) {
  951.         sival.ivalue = BASE1 - ((FULL) *h1++) + ((FULL) *h2++) + carry;
  952.         *hd++ = BASE1 - sival.silow;
  953.         carry = sival.sihigh;
  954.     }
  955.     while (copysize--) {
  956.         sival.ivalue = (BASE1 - ((FULL) *h1++)) + carry;
  957.         *hd++ = BASE1 - sival.silow;
  958.         carry = sival.sihigh;
  959.     }
  960.  
  961.     hd = &baseAB[sizeAB - 1];
  962.     while ((*hd == 0) && (sizeAB > 1)) {
  963.         sizeAB--;
  964.         hd--;
  965.     }
  966.  
  967.     /*
  968.      * Now square the number into another temporary location,
  969.      * and subtract that from the final result.
  970.      */
  971.     sizeABAB = dosquare(baseAB, sizeAB, baseABAB);
  972.  
  973.     h1 = baseABAB;
  974.     hd = ans + shift;
  975.     carry = 0;
  976.     while (sizeABAB--) {
  977.         sival.ivalue = BASE1 - ((FULL) *hd) + ((FULL) *h1++) + carry;
  978.         *hd++ = BASE1 - sival.silow;
  979.         carry = sival.sihigh;
  980.     }
  981.     while (carry) {
  982.         sival.ivalue = BASE1 - ((FULL) *hd) + carry;
  983.         *hd++ = BASE1 - sival.silow;
  984.         carry = sival.sihigh;
  985.     }
  986.  
  987.     /*
  988.      * Return the size of the result.
  989.      */
  990.     len = sizetotal;
  991.     hd = &ans[len - 1];
  992.     while ((*hd == 0) && (len > 1)) {
  993.         len--;
  994.         hd--;
  995.     }
  996.     tempbuf = temp;
  997.     return len;
  998. }
  999.  
  1000.  
  1001. /*
  1002.  * Return a pointer to a buffer to be used for holding a temporary number.
  1003.  * The buffer will be at least as large as the specified number of HALFs,
  1004.  * and remains valid until the next call to this routine.  The buffer cannot
  1005.  * be freed by the caller.  There is only one temporary buffer, and so as to
  1006.  * avoid possible conflicts this is only used by the lowest level routines
  1007.  * such as divide, multiply, and square.
  1008.  */
  1009. HALF *
  1010. zalloctemp(len)
  1011.     LEN len;        /* required number of HALFs in buffer */
  1012. {
  1013.     HALF *hp;
  1014.     static LEN buflen;    /* current length of temp buffer */
  1015.     static HALF *bufptr;    /* pointer to current temp buffer */
  1016.  
  1017.     if (len <= buflen)
  1018.         return bufptr;
  1019.  
  1020.     /*
  1021.      * We need to grow the temporary buffer.
  1022.      * First free any existing buffer, and then allocate the new one.
  1023.      * While we are at it, make the new buffer bigger than necessary
  1024.      * in order to reduce the number of reallocations.
  1025.      */
  1026.     len += 100;
  1027.     if (buflen) {
  1028.         buflen = 0;
  1029.         free(bufptr);
  1030.     }
  1031.     hp = (HALF *) malloc(len * sizeof(HALF));
  1032.     if (hp == NULL)
  1033.         error("No memory for temp buffer");
  1034.     bufptr = hp;
  1035.     buflen = len;
  1036.     return hp;
  1037. }
  1038.  
  1039. /* END CODE */
  1040.