home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 2 / Apprentice-Release2.iso / Source Code / C / Utilities / Calc 1.24.7 / zmod.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-10-10  |  31.7 KB  |  1,185 lines  |  [TEXT/R*ch]
  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.  * Routines to do modulo arithmetic both normally and also using the REDC
  7.  * algorithm given by Peter L. Montgomery in Mathematics of Computation,
  8.  * volume 44, number 170 (April, 1985).  For multiple multiplies using
  9.  * the same large modulus, the REDC algorithm avoids the usual division
  10.  * by the modulus, instead replacing it with two multiplies or else a
  11.  * special algorithm.  When these two multiplies or the special algorithm
  12.  * are faster then the division, then the REDC algorithm is better.
  13.  */
  14.  
  15. #include "xmath.h"
  16.  
  17.  
  18. #define    POWBITS    4        /* bits for power chunks (must divide BASEB) */
  19. #define    POWNUMS    (1<<POWBITS)    /* number of powers needed in table */
  20.  
  21.  
  22. LEN _pow2_ = POW_ALG2;        /* modulo size to use REDC for powers */
  23. LEN _redc2_ = REDC_ALG2;    /* modulo size to use second REDC algorithm */
  24.  
  25. static REDC *powermodredc = NULL;    /* REDC info for raising to power */
  26.  
  27. #if 0
  28. extern void zaddmod proto((ZVALUE z1, ZVALUE z2, ZVALUE z3, ZVALUE *res));
  29. extern void znegmod proto((ZVALUE z1, ZVALUE z2, ZVALUE *res));
  30.  
  31. /*
  32.  * Multiply two numbers together and then mod the result with a third number.
  33.  * The two numbers to be multiplied can be negative or out of modulo range.
  34.  * The result will be in the range 0 to the modulus - 1.
  35.  */
  36. void
  37. zmulmod(z1, z2, z3, res)
  38.     ZVALUE z1;        /* first number to be multiplied */
  39.     ZVALUE z2;        /* second number to be multiplied */
  40.     ZVALUE z3;        /* number to take mod with */
  41.     ZVALUE *res;        /* result */
  42. {
  43.     ZVALUE tmp;
  44.     FULL prod;
  45.     FULL digit;
  46.     BOOL neg;
  47.  
  48.     if (iszero(z3) || isneg(z3))
  49.         error("Mod of non-positive integer");
  50.     if (iszero(z1) || iszero(z2) || isunit(z3)) {
  51.         *res = _zero_;
  52.         return;
  53.     }
  54.  
  55.     /*
  56.      * If the modulus is a single digit number, then do the result
  57.      * cheaply.  Check especially for a small power of two.
  58.      */
  59.     if (istiny(z3)) {
  60.         neg = (z1.sign != z2.sign);
  61.         digit = z3.v[0];
  62.         if ((digit & -digit) == digit) {    /* NEEDS 2'S COMP */
  63.             prod = ((FULL) z1.v[0]) * ((FULL) z2.v[0]);
  64.             prod &= (digit - 1);
  65.         } else {
  66.             z1.sign = 0;
  67.             z2.sign = 0;
  68.             prod = (FULL) zmodi(z1, (long) digit);
  69.             prod *= (FULL) zmodi(z2, (long) digit);
  70.             prod %= digit;
  71.         }
  72.         if (neg && prod)
  73.             prod = digit - prod;
  74.         itoz((long) prod, res);
  75.         return;
  76.     }
  77.  
  78.     /*
  79.      * The modulus is more than one digit.
  80.      * Actually do the multiply and divide if necessary.
  81.      */
  82.     zmul(z1, z2, &tmp);
  83.     if (ispos(tmp) && ((tmp.len < z3.len) || ((tmp.len == z3.len) &&
  84.         (tmp.v[tmp.len-1] < z2.v[z3.len-1]))))
  85.     {
  86.         *res = tmp;
  87.         return;
  88.     }
  89.     zmod(tmp, z3, res);
  90.     freeh(tmp.v);
  91. }
  92.  
  93.  
  94. /*
  95.  * Square a number and then mod the result with a second number.
  96.  * The number to be squared can be negative or out of modulo range.
  97.  * The result will be in the range 0 to the modulus - 1.
  98.  */
  99. void
  100. zsquaremod(z1, z2, res)
  101.     ZVALUE z1;        /* number to be squared */
  102.     ZVALUE z2;        /* number to take mod with */
  103.     ZVALUE *res;        /* result */
  104. {
  105.     ZVALUE tmp;
  106.     FULL prod;
  107.     FULL digit;
  108.  
  109.     if (iszero(z2) || isneg(z2))
  110.         error("Mod of non-positive integer");
  111.     if (iszero(z1) || isunit(z2)) {
  112.         *res = _zero_;
  113.         return;
  114.     }
  115.  
  116.     /*
  117.      * If the modulus is a single digit number, then do the result
  118.      * cheaply.  Check especially for a small power of two.
  119.      */
  120.     if (istiny(z2)) {
  121.         digit = z2.v[0];
  122.         if ((digit & -digit) == digit) {    /* NEEDS 2'S COMP */
  123.             prod = (FULL) z1.v[0];
  124.             prod = (prod * prod) & (digit - 1);
  125.         } else {
  126.             z1.sign = 0;
  127.             prod = (FULL) zmodi(z1, (long) digit);
  128.             prod = (prod * prod) % digit;
  129.         }
  130.         itoz((long) prod, res);
  131.         return;
  132.     }
  133.  
  134.     /*
  135.      * The modulus is more than one digit.
  136.      * Actually do the square and divide if necessary.
  137.      */
  138.     zsquare(z1, &tmp);
  139.     if ((tmp.len < z2.len) ||
  140.         ((tmp.len == z2.len) && (tmp.v[tmp.len-1] < z2.v[z2.len-1]))) {
  141.             *res = tmp;
  142.             return;
  143.     }
  144.     zmod(tmp, z2, res);
  145.     freeh(tmp.v);
  146. }
  147.  
  148.  
  149. /*
  150.  * Add two numbers together and then mod the result with a third number.
  151.  * The two numbers to be added can be negative or out of modulo range.
  152.  * The result will be in the range 0 to the modulus - 1.
  153.  */
  154. static void
  155. zaddmod(z1, z2, z3, res)
  156.     ZVALUE z1;        /* first number to be added */
  157.     ZVALUE z2;        /* second number to be added */
  158.     ZVALUE z3;        /* number to take mod with */
  159.     ZVALUE *res;        /* result */
  160. {
  161.     ZVALUE tmp;
  162.     FULL sumdigit;
  163.     FULL moddigit;
  164.  
  165.     if (iszero(z3) || isneg(z3))
  166.         error("Mod of non-positive integer");
  167.     if ((iszero(z1) && iszero(z2)) || isunit(z3)) {
  168.         *res = _zero_;
  169.         return;
  170.     }
  171.     if (istwo(z2)) {
  172.         if ((z1.v[0] + z2.v[0]) & 0x1)
  173.             *res = _one_;
  174.         else
  175.             *res = _zero_;
  176.         return;
  177.     }
  178.     zadd(z1, z2, &tmp);
  179.     if (isneg(tmp) || (tmp.len > z3.len)) {
  180.         zmod(tmp, z3, res);
  181.         freeh(tmp.v);
  182.         return;
  183.     }
  184.     sumdigit = tmp.v[tmp.len - 1];
  185.     moddigit = z3.v[z3.len - 1];
  186.     if ((tmp.len < z3.len) || (sumdigit < moddigit)) {
  187.         *res = tmp;
  188.         return;
  189.     }
  190.     if (sumdigit < 2 * moddigit) {
  191.         zsub(tmp, z3, res);
  192.         freeh(tmp.v);
  193.         return;
  194.     }
  195.     zmod(tmp, z2, res);
  196.     freeh(tmp.v);
  197. }
  198.  
  199.  
  200. /*
  201.  * Subtract two numbers together and then mod the result with a third number.
  202.  * The two numbers to be subtract can be negative or out of modulo range.
  203.  * The result will be in the range 0 to the modulus - 1.
  204.  */
  205. void
  206. zsubmod(z1, z2, z3, res)
  207.     ZVALUE z1;        /* number to be subtracted from */
  208.     ZVALUE z2;        /* number to be subtracted */
  209.     ZVALUE z3;        /* number to take mod with */
  210.     ZVALUE *res;        /* result */
  211. {
  212.     if (iszero(z3) || isneg(z3))
  213.         error("Mod of non-positive integer");
  214.     if (iszero(z2)) {
  215.         zmod(z1, z3, res);
  216.         return;
  217.     }
  218.     if (iszero(z1)) {
  219.         znegmod(z2, z3, res);
  220.         return;
  221.     }
  222.     if ((z1.sign == z2.sign) && (z1.len == z2.len) &&
  223.         (z1.v[0] == z2.v[0]) && (zcmp(z1, z2) == 0)) {
  224.             *res = _zero_;
  225.             return;
  226.     }
  227.     z2.sign = !z2.sign;
  228.     zaddmod(z1, z2, z3, res);
  229. }
  230.  
  231.  
  232. /*
  233.  * Calculate the negative of a number modulo another number.
  234.  * The number to be negated can be negative or out of modulo range.
  235.  * The result will be in the range 0 to the modulus - 1.
  236.  */
  237. static void
  238. znegmod(z1, z2, res)
  239.     ZVALUE z1;        /* number to take negative of */
  240.     ZVALUE z2;        /* number to take mod with */
  241.     ZVALUE *res;        /* result */
  242. {
  243.     int sign;
  244.     int cv;
  245.  
  246.     if (iszero(z2) || isneg(z2))
  247.         error("Mod of non-positive integer");
  248.     if (iszero(z1) || isunit(z2)) {
  249.         *res = _zero_;
  250.         return;
  251.     }
  252.     if (istwo(z2)) {
  253.         if (z1.v[0] & 0x1)
  254.             *res = _one_;
  255.         else
  256.             *res = _zero_;
  257.         return;
  258.     }
  259.  
  260.     /*
  261.      * If the absolute value of the number is within the modulo range,
  262.      * then the result is just a copy or a subtraction.  Otherwise go
  263.      * ahead and negate and reduce the result.
  264.      */
  265.     sign = z1.sign;
  266.     z1.sign = 0;
  267.     cv = zrel(z1, z2);
  268.     if (cv == 0) {
  269.         *res = _zero_;
  270.         return;
  271.     }
  272.     if (cv < 0) {
  273.         if (sign)
  274.             zcopy(z1, res);
  275.         else
  276.             zsub(z2, z1, res);
  277.         return;
  278.     }
  279.     z1.sign = !sign;
  280.     zmod(z1, z2, res);
  281. }
  282. #endif
  283.  
  284.  
  285. /*
  286.  * Calculate the number congruent to the given number whose absolute
  287.  * value is minimal.  The number to be reduced can be negative or out of
  288.  * modulo range.  The result will be within the range -int((modulus-1)/2)
  289.  * to int(modulus/2) inclusive.  For example, for modulus 7, numbers are
  290.  * reduced to the range [-3, 3], and for modulus 8, numbers are reduced to
  291.  * the range [-3, 4].
  292.  */
  293. void
  294. zminmod(z1, z2, res)
  295.     ZVALUE z1;        /* number to find minimum congruence of */
  296.     ZVALUE z2;        /* number to take mod with */
  297.     ZVALUE *res;        /* result */
  298. {
  299.     ZVALUE tmp1, tmp2;
  300.     int sign;
  301.     int cv;
  302.  
  303.     if (iszero(z2) || isneg(z2))
  304.         error("Mod of non-positive integer");
  305.     if (iszero(z1) || isunit(z2)) {
  306.         *res = _zero_;
  307.         return;
  308.     }
  309.     if (istwo(z2)) {
  310.         if (isodd(z1))
  311.             *res = _one_;
  312.         else
  313.             *res = _zero_;
  314.         return;
  315.     }
  316.  
  317.     /*
  318.      * Do a quick check to see if the number is very small compared
  319.      * to the modulus.  If so, then the result is obvious.
  320.      */
  321.     if (z1.len < z2.len - 1) {
  322.         zcopy(z1, res);
  323.         return;
  324.     }
  325.  
  326.     /*
  327.      * Now make sure the input number is within the modulo range.
  328.      * If not, then reduce it to be within range and make the
  329.      * quick check again.
  330.      */
  331.     sign = z1.sign;
  332.     z1.sign = 0;
  333.     cv = zrel(z1, z2);
  334.     if (cv == 0) {
  335.         *res = _zero_;
  336.         return;
  337.     }
  338.     tmp1 = z1;
  339.     if (cv > 0) {
  340.         z1.sign = (BOOL)sign;
  341.         zmod(z1, z2, &tmp1);
  342.         if (tmp1.len < z2.len - 1) {
  343.             *res = tmp1;
  344.             return;
  345.         }
  346.         sign = 0;
  347.     }
  348.  
  349.     /*
  350.      * Now calculate the difference of the modulus and the absolute
  351.      * value of the original number.  Compare the original number with
  352.      * the difference, and return the one with the smallest absolute
  353.      * value, with the correct sign.  If the two values are equal, then
  354.      * return the positive result.
  355.      */
  356.     zsub(z2, tmp1, &tmp2);
  357.     cv = zrel(tmp1, tmp2);
  358.     if (cv < 0) {
  359.         freeh(tmp2.v);
  360.         tmp1.sign = (BOOL)sign;
  361.         if (tmp1.v == z1.v)
  362.             zcopy(tmp1, res);
  363.         else
  364.             *res = tmp1;
  365.     } else {
  366.         if (cv)
  367.             tmp2.sign = !sign;
  368.         if (tmp1.v != z1.v)
  369.             freeh(tmp1.v);
  370.         *res = tmp2;
  371.     }
  372. }
  373.  
  374.  
  375. /*
  376.  * Compare two numbers for equality modulo a third number.
  377.  * The two numbers to be compared can be negative or out of modulo range.
  378.  * Returns TRUE if the numbers are not congruent, and FALSE if they are
  379.  * congruent.
  380.  */
  381. BOOL
  382. zcmpmod(z1, z2, z3)
  383.     ZVALUE z1;        /* first number to be compared */
  384.     ZVALUE z2;        /* second number to be compared */
  385.     ZVALUE z3;        /* modulus */
  386. {
  387.     ZVALUE tmp1, tmp2, tmp3;
  388.     FULL digit;
  389.     LEN len;
  390.     int cv;
  391.  
  392.     if (isneg(z3) || iszero(z3))
  393.         error("Non-positive modulus in zcmpmod");
  394.     if (istwo(z3))
  395.         return (((z1.v[0] + z2.v[0]) & 0x1) != 0);
  396.  
  397.     /*
  398.      * If the two numbers are equal, then their mods are equal.
  399.      */
  400.     if ((z1.sign == z2.sign) && (z1.len == z2.len) &&
  401.         (z1.v[0] == z2.v[0]) && (zcmp(z1, z2) == 0))
  402.             return FALSE;
  403.  
  404.     /*
  405.      * If both numbers are negative, then we can make them positive.
  406.      */
  407.     if (isneg(z1) && isneg(z2)) {
  408.         z1.sign = 0;
  409.         z2.sign = 0;
  410.     }
  411.  
  412.     /*
  413.      * For small negative numbers, make them positive before comparing.
  414.      * In any case, the resulting numbers are in tmp1 and tmp2.
  415.      */
  416.     tmp1 = z1;
  417.     tmp2 = z2;
  418.     len = z3.len;
  419.     digit = z3.v[len - 1];
  420.  
  421.     if (isneg(z1) && ((z1.len < len) ||
  422.         ((z1.len == len) && (z1.v[z1.len - 1] < digit))))
  423.             zadd(z1, z3, &tmp1);
  424.  
  425.     if (isneg(z2) && ((z2.len < len) ||
  426.         ((z2.len == len) && (z2.v[z2.len - 1] < digit))))
  427.             zadd(z2, z3, &tmp2);
  428.  
  429.     /*
  430.      * Now compare the two numbers for equality.
  431.      * If they are equal we are all done.
  432.      */
  433.     if (zcmp(tmp1, tmp2) == 0) {
  434.         if (tmp1.v != z1.v)
  435.             freeh(tmp1.v);
  436.         if (tmp2.v != z2.v)
  437.             freeh(tmp2.v);
  438.         return FALSE;
  439.     }
  440.  
  441.     /*
  442.      * They are not identical.  Now if both numbers are positive
  443.      * and less than the modulus, then they are definitely not equal.
  444.      */
  445.     if ((tmp1.sign == tmp2.sign) &&
  446.         ((tmp1.len < len) || (zrel(tmp1, z3) < 0)) &&
  447.         ((tmp2.len < len) || (zrel(tmp2, z3) < 0)))
  448.     {
  449.         if (tmp1.v != z1.v)
  450.             freeh(tmp1.v);
  451.         if (tmp2.v != z2.v)
  452.             freeh(tmp2.v);
  453.         return TRUE;
  454.     }
  455.  
  456.     /*
  457.      * Either one of the numbers is negative or is large.
  458.      * So do the standard thing and subtract the two numbers.
  459.      * Then they are equal if the result is 0 (mod z3).
  460.      */
  461.     zsub(tmp1, tmp2, &tmp3);
  462.     if (tmp1.v != z1.v)
  463.         freeh(tmp1.v);
  464.     if (tmp2.v != z2.v)
  465.         freeh(tmp2.v);
  466.  
  467.     /*
  468.      * Compare the result with the modulus to see if it is equal to
  469.      * or less than the modulus.  If so, we know the mod result.
  470.      */
  471.     tmp3.sign = 0;
  472.     cv = zrel(tmp3, z3);
  473.     if (cv == 0) {
  474.         freeh(tmp3.v);
  475.         return FALSE;
  476.     }
  477.     if (cv < 0) {
  478.         freeh(tmp3.v);
  479.         return TRUE;
  480.     }
  481.  
  482.     /*
  483.      * We are forced to actually do the division.
  484.      * The numbers are congruent if the result is zero.
  485.      */
  486.     zmod(tmp3, z3, &tmp1);
  487.     freeh(tmp3.v);
  488.     if (iszero(tmp1)) {
  489.         freeh(tmp1.v);
  490.         return FALSE;
  491.     } else {
  492.         freeh(tmp1.v);
  493.         return TRUE;
  494.     }
  495. }
  496.  
  497.  
  498. /*
  499.  * Compute the result of raising one number to a power modulo another number.
  500.  * That is, this computes:  a^b (modulo c).
  501.  * This calculates the result by examining the power POWBITS bits at a time,
  502.  * using a small table of POWNUMS low powers to calculate powers for those bits,
  503.  * and repeated squaring and multiplying by the partial powers to generate
  504.  * the complete power.  If the power being raised to is high enough, then
  505.  * this uses the REDC algorithm to avoid doing many divisions.  When using
  506.  * REDC, multiple calls to this routine using the same modulus will be
  507.  * slightly faster.
  508.  */
  509. void
  510. zpowermod(z1, z2, z3, res)
  511.     ZVALUE z1, z2, z3, *res;
  512. {
  513.     HALF *hp;        /* pointer to current word of the power */
  514.     REDC *rp;        /* REDC information to be used */
  515.     ZVALUE *pp;        /* pointer to low power table */
  516.     ZVALUE ans, temp;    /* calculation values */
  517.     ZVALUE modpow;        /* current small power */
  518.     ZVALUE lowpowers[POWNUMS];    /* low powers */
  519.     int sign;        /* original sign of number */
  520.     int curshift;        /* shift value for word of power */
  521.     HALF curhalf;        /* current word of power */
  522.     unsigned int curpow;    /* current low power */
  523.     unsigned int curbit;    /* current bit of low power */
  524.     int i;
  525.  
  526.     if (isneg(z3) || iszero(z3))
  527.         error("Non-positive modulus in zpowermod");
  528.     if (isneg(z2))
  529.         error("Negative power in zpowermod");
  530.  
  531.     sign = z1.sign;
  532.     z1.sign = 0;
  533.  
  534.     /*
  535.      * Check easy cases first.
  536.      */
  537.     if (iszero(z1) || isunit(z3)) {        /* 0^x or mod 1 */
  538.         *res = _zero_;
  539.         return;
  540.     }
  541.     if (istwo(z3)) {            /* mod 2 */
  542.         if (isodd(z1))
  543.             *res = _one_;
  544.         else
  545.             *res = _zero_;
  546.         return;
  547.     }
  548.     if (isunit(z1) && (!sign || iseven(z2))) {
  549.         /* 1^x or (-1)^(2x) */
  550.         *res = _one_;
  551.         return;
  552.     }
  553.  
  554.     /*
  555.      * Normalize the number being raised to be non-negative and to lie
  556.      * within the modulo range.  Then check for zero or one specially.
  557.      */
  558.     zmod(z1, z3, &temp);
  559.     if (iszero(temp)) {
  560.         freeh(temp.v);
  561.         *res = _zero_;
  562.         return;
  563.     }
  564.     z1 = temp;
  565.     if (sign) {
  566.         zsub(z3, z1, &temp);
  567.         freeh(z1.v);
  568.         z1 = temp;
  569.     }
  570.     if (isunit(z1)) {
  571.         freeh(z1.v);
  572.         *res = _one_;
  573.         return;
  574.     }
  575.  
  576.     /*
  577.      * If the modulus is odd, large enough, is not one less than an
  578.      * exact power of two, and if the power is large enough, then use
  579.      * the REDC algorithm.  The size where this is done is configurable.
  580.      */
  581.     if ((z2.len > 1) && (z3.len >= _pow2_) && isodd(z3)
  582.         && !zisallbits(z3))
  583.     {
  584.         if (powermodredc && zcmp(powermodredc->mod, z3)) {
  585.             zredcfree(powermodredc);
  586.             powermodredc = NULL;
  587.         }
  588.         if (powermodredc == NULL)
  589.             powermodredc = zredcalloc(z3);
  590.         rp = powermodredc;
  591.         zredcencode(rp, z1, &temp);
  592.         zredcpower(rp, temp, z2, &z1);
  593.         freeh(temp.v);
  594.         zredcdecode(rp, z1, res);
  595.         freeh(z1.v);
  596.         return;
  597.     }
  598.  
  599.     /*
  600.      * Modulus or power is small enough to perform the power raising
  601.      * directly.  Initialize the table of powers.
  602.      */
  603.     for (pp = &lowpowers[2]; pp < &lowpowers[POWNUMS]; pp++)
  604.         pp->len = 0;
  605.     lowpowers[0] = _one_;
  606.     lowpowers[1] = z1;
  607.     ans = _one_;
  608.  
  609.     hp = &z2.v[z2.len - 1];
  610.     curhalf = *hp;
  611.     curshift = BASEB - POWBITS;
  612.     while (curshift && ((curhalf >> curshift) == 0))
  613.         curshift -= POWBITS;
  614.  
  615.     /*
  616.      * Calculate the result by examining the power POWBITS bits at a time,
  617.      * and use the table of low powers at each iteration.
  618.      */
  619.     for (;;) {
  620.         curpow = (curhalf >> curshift) & (POWNUMS - 1);
  621.         pp = &lowpowers[curpow];
  622.  
  623.         /*
  624.          * If the small power is not yet saved in the table, then
  625.          * calculate it and remember it in the table for future use.
  626.          */
  627.         if (pp->len == 0) {
  628.             if (curpow & 0x1)
  629.                 zcopy(z1, &modpow);
  630.             else
  631.                 modpow = _one_;
  632.  
  633.             for (curbit = 0x2; curbit <= curpow; curbit *= 2) {
  634.                 pp = &lowpowers[curbit];
  635.                 if (pp->len == 0) {
  636.                     zsquare(lowpowers[curbit/2], &temp);
  637.                     zmod(temp, z3, pp);
  638.                     freeh(temp.v);
  639.                 }
  640.                 if (curbit & curpow) {
  641.                     zmul(*pp, modpow, &temp);
  642.                     freeh(modpow.v);
  643.                     zmod(temp, z3, &modpow);
  644.                     freeh(temp.v);
  645.                 }
  646.             }
  647.             pp = &lowpowers[curpow];
  648.             *pp = modpow;
  649.         }
  650.  
  651.         /*
  652.          * If the power is nonzero, then accumulate the small power
  653.          * into the result.
  654.          */
  655.         if (curpow) {
  656.             zmul(ans, *pp, &temp);
  657.             freeh(ans.v);
  658.             zmod(temp, z3, &ans);
  659.             freeh(temp.v);
  660.         }
  661.  
  662.         /*
  663.          * Select the next POWBITS bits of the power, if there is
  664.          * any more to generate.
  665.          */
  666.         curshift -= POWBITS;
  667.         if (curshift < 0) {
  668.             if (hp-- == z2.v)
  669.                 break;
  670.             curhalf = *hp;
  671.             curshift = BASEB - POWBITS;
  672.         }
  673.  
  674.         /*
  675.          * Square the result POWBITS times to make room for the next
  676.          * chunk of bits.
  677.          */
  678.         for (i = 0; i < POWBITS; i++) {
  679.             zsquare(ans, &temp);
  680.             freeh(ans.v);
  681.             zmod(temp, z3, &ans);
  682.             freeh(temp.v);
  683.         }
  684.     }
  685.  
  686.     for (pp = &lowpowers[2]; pp < &lowpowers[POWNUMS]; pp++) {
  687.         if (pp->len)
  688.             freeh(pp->v);
  689.     }
  690.     *res = ans;
  691. }
  692.  
  693.  
  694. /*
  695.  * Initialize the REDC algorithm for a particular modulus,
  696.  * returning a pointer to a structure that is used for other
  697.  * REDC calls.  An error is generated if the structure cannot
  698.  * be allocated.  The modulus must be odd and positive.
  699.  */
  700. REDC *
  701. zredcalloc(z1)
  702.     ZVALUE z1;        /* modulus to initialize for */
  703. {
  704.     REDC *rp;        /* REDC information */
  705.     ZVALUE tmp;
  706.     long bit;
  707.  
  708.     if (iseven(z1) || isneg(z1))
  709.         error("REDC requires positive odd modulus");
  710.  
  711.     rp = (REDC *) malloc(sizeof(REDC));
  712.     if (rp == NULL)
  713.         error("Cannot allocate REDC structure");
  714.  
  715.     /*
  716.      * Round up the binary modulus to the next power of two
  717.      * which is at a word boundary.  Then the shift and modulo
  718.      * operations mod the binary modulus can be done very cheaply.
  719.      * Calculate the REDC format for the number 1 for future use.
  720.      */
  721.     bit = zhighbit(z1) + 1;
  722.     if (bit % BASEB)
  723.         bit += (BASEB - (bit % BASEB));
  724.     zcopy(z1, &rp->mod);
  725.     zbitvalue(bit, &tmp);
  726.     z1.sign = 1;
  727.     (void) zmodinv(z1, tmp, &rp->inv);
  728.     zmod(tmp, rp->mod, &rp->one);
  729.     freeh(tmp.v);
  730.     rp->len = bit / BASEB;
  731.     return rp;
  732. }
  733.  
  734.  
  735. /*
  736.  * Free any numbers associated with the specified REDC structure,
  737.  * and then the REDC structure itself.
  738.  */
  739. void
  740. zredcfree(rp)
  741.     REDC *rp;        /* REDC information to be cleared */
  742. {
  743.     freeh(rp->mod.v);
  744.     freeh(rp->inv.v);
  745.     freeh(rp->one.v);
  746.     free(rp);
  747. }
  748.  
  749.  
  750. /*
  751.  * Convert a normal number into the specified REDC format.
  752.  * The number to be converted can be negative or out of modulo range.
  753.  * The resulting number can be used for multiplying, adding, subtracting,
  754.  * or comparing with any other such converted numbers, as if the numbers
  755.  * were being calculated modulo the number which initialized the REDC
  756.  * information.  When the final value is unconverted, the result is the
  757.  * same as if the usual operations were done with the original numbers.
  758.  */
  759. void
  760. zredcencode(rp, z1, res)
  761.     REDC *rp;        /* REDC information */
  762.     ZVALUE z1;        /* number to be converted */
  763.     ZVALUE *res;        /* returned converted number */
  764. {
  765.     ZVALUE tmp1, tmp2;
  766.  
  767.     /*
  768.      * Handle the cases 0, 1, -1, and 2 specially since these are
  769.      * easy to calculate.  Zero transforms to zero, and the others
  770.      * can be obtained from the precomputed REDC format for 1 since
  771.      * addition and subtraction act normally for REDC format numbers.
  772.      */
  773.     if (iszero(z1)) {
  774.         *res = _zero_;
  775.         return;
  776.     }
  777.     if (isone(z1)) {
  778.         zcopy(rp->one, res);
  779.         return;
  780.     }
  781.     if (isunit(z1)) {
  782.         zsub(rp->mod, rp->one, res);
  783.         return;
  784.     }
  785.     if (istwo(z1)) {
  786.         zadd(rp->one, rp->one, &tmp1);
  787.         if (zrel(tmp1, rp->mod) < 0) {
  788.             *res = tmp1;
  789.             return;
  790.         }
  791.         zsub(tmp1, rp->mod, res);
  792.         freeh(tmp1.v);
  793.         return;
  794.     }
  795.  
  796.     /*
  797.      * Not a trivial number to convert, so do the full transformation.
  798.      * Convert negative numbers to positive numbers before converting.
  799.      */
  800.     tmp1.len = 0;
  801.     if (isneg(z1)) {
  802.         zmod(z1, rp->mod, &tmp1);
  803.         z1 = tmp1;
  804.     }
  805.     zshift(z1, rp->len * BASEB, &tmp2);
  806.     if (tmp1.len)
  807.         freeh(tmp1.v);
  808.     zmod(tmp2, rp->mod, res);
  809.     freeh(tmp2.v);
  810. }
  811.  
  812.  
  813. /*
  814.  * The REDC algorithm used to convert numbers out of REDC format and also
  815.  * used after multiplication of two REDC numbers.  Using this routine
  816.  * avoids any divides, replacing the divide by two multiplications.
  817.  * If the numbers are very large, then these two multiplies will be
  818.  * quicker than the divide, since dividing is harder than multiplying.
  819.  */
  820. void
  821. zredcdecode(rp, z1, res)
  822.     REDC *rp;        /* REDC information */
  823.     ZVALUE z1;        /* number to be transformed */
  824.     ZVALUE *res;        /* returned transformed number */
  825. {
  826.     ZVALUE tmp1, tmp2;    /* temporaries */
  827.     HALF *hp;        /* saved pointer to tmp2 value */
  828.  
  829.     if (isneg(z1))
  830.         error("Negative number for zredc");
  831.  
  832.     /*
  833.      * Check first for the special values for 0 and 1 that are easy.
  834.      */
  835.     if (iszero(z1)) {
  836.         *res = _zero_;
  837.         return;
  838.     }
  839.     if ((z1.len == rp->one.len) && (z1.v[0] == rp->one.v[0]) &&
  840.         (zcmp(z1, rp->one) == 0)) {
  841.             *res = _one_;
  842.             return;
  843.     }
  844.  
  845.     /*
  846.      * First calculate the following:
  847.      *     tmp2 = ((z1 % 2^bitnum) * inv) % 2^bitnum.
  848.      * The mod operations can be done with no work since the bit
  849.      * number was selected as a multiple of the word size.  Just
  850.      * reduce the sizes of the numbers as required.
  851.      */
  852.     tmp1 = z1;
  853.     if (tmp1.len > rp->len)
  854.         tmp1.len = rp->len;
  855.     zmul(tmp1, rp->inv, &tmp2);
  856.     if (tmp2.len > rp->len)
  857.         tmp2.len = rp->len;
  858.  
  859.     /*
  860.      * Next calculate the following:
  861.      *    res = (z1 + tmp2 * modulus) / 2^bitnum
  862.      * The division by a power of 2 is always exact, and requires no
  863.      * work.  Just adjust the address and length of the number to do
  864.      * the divide, but save the original pointer for freeing later.
  865.      */
  866.     zmul(tmp2, rp->mod, &tmp1);
  867.     freeh(tmp2.v);
  868.     zadd(z1, tmp1, &tmp2);
  869.     freeh(tmp1.v);
  870.     hp = tmp2.v;
  871.     if (tmp2.len <= rp->len) {
  872.         freeh(hp);
  873.         *res = _zero_;
  874.         return;
  875.     }
  876.     tmp2.v += rp->len;
  877.     tmp2.len -= rp->len;
  878.  
  879.     /*
  880.      * Finally do a final modulo by a simple subtraction if necessary.
  881.      * This is all thanterations is equal to the size of the modulus.
  882.          * This acts as if the innermost loop was repeated for
  883.          * high digits of z2 that are zero.
  884.          */
  885.         len2 = modlen - z2.len;
  886.         while (len2--) {
  887.             sival2.ivalue = muln * ((FULL) *h3++);
  888.             sival3.ivalue = ((FULL) sival2.silow)
  889.                 + ((FULL) *hd)
  890.                 + ((FULL) carry.silow);
  891.             carry.ivalue = ((FULL) sival2.sihigh)
  892.                 + ((FULL) sival3.sihigh)
  893.                 + ((FULL) carry.sihigh);
  894.  
  895.             hd[-1] = sival3.silow;
  896.             hd++;
  897.         }
  898.  
  899.         res->v[modlen - 1] = carry.silow;
  900.         topdigit = carry.sihigh;
  901.     }
  902.  
  903.     /*
  904.      * Now continue the loop as necessary so the total number
  905.      * of interations is equal to the size of the modulus.
  906.      * This acts as if the outermost loop was repeated for high
  907.      * digits of z1 that are zero.
  908.      */
  909.     len = modlen - z1.len;
  910.     while (len--) {
  911.         /*
  912.          * Start off with the first digit of the modulus.
  913.          */
  914.         h3 = rp->mod.v;
  915.         hd = res->v;
  916.         muln = ((HALF) (*hd * Ninv));
  917.         sival2.ivalue = muln * ((FULL) *h3++);
  918.         sival3.ivalue = ((FULL) *hd++) + ((FULL) sival2.silow);
  919.         carry.ivalue = ((FULL) sival2.sihigh) + ((FULL) sival3.sihigh);
  920.  
  921.         /*
  922.          * Do this innermost loop for each digit of the modulus,
  923.          * except for the first digit which was just done above.
  924.          */
  925.         len2 = modlen;
  926.         while (--len2 > 0) {
  927.             sival2.ivalue = muln * ((FULL) *h3++);
  928.             sival3.ivalue = ((FULL) sival2.silow)
  929.                 + ((FULL) *hd)
  930.                 + ((FULL) carry.silow);
  931.             carry.ivalue = ((FULL) sival2.sihigh)
  932.                 + ((FULL) sival3.sihigh)
  933.                 + ((FULL) carry.sihigh);
  934.  
  935.             hd[-1] = sival3.silow;
  936.             hd++;
  937.         }
  938.         res->v[modlen - 1] = carry.silow;
  939.         topdigit = carry.sihigh;
  940.     }
  941.  
  942.     /*
  943.      * Determine the true size of the result, taking the top digit of
  944.      * the current result into account.  The top digit is not stored in
  945.      * the number because it is temporary and would become zero anyway
  946.      * after the final subtraction is done.
  947.      */
  948.     if (topdigit == 0) {
  949.         len = modlen;
  950.         hd = &res->v[len - 1];
  951.         while ((*hd == 0) && (len > 1)) {
  952.             hd--;
  953.             len--;
  954.         }
  955.         res->len = len;
  956.     }
  957.  
  958.     /*
  959.      * Compare the result with the modulus.
  960.      * If it is less than the modulus, then the calculation is complete.
  961.      */
  962.     if ((topdigit == 0) && ((len < modlen)
  963.         || (res->v[len - 1] < rp->mod.v[len - 1])
  964.         || (zrel(*res, rp->mod) < 0)))
  965.             return;
  966.  
  967.     /*
  968.      * Do a subtraction to reduce the result to a value less than
  969.      * the modulus.  The REDC algorithm guarantees that a single subtract
  970.      * is all that is needed.  Ignore any borrowing from the possible
  971.      * highest word of the current result because that would affect
  972.      * only the top digit value that was not stored and would become
  973.      * zero anyway.
  974.      */
  975.     carry.ivalue = 0;
  976.     h1 = rp->mod.v;
  977.     hd = res->v;
  978.     len = modlen;
  979.     while (len--) {
  980.         carry.ivalue = BASE1 - ((FULL) *hd) + ((FULL) *h1++)
  981.             + ((FULL) carry.silow);
  982.         *hd++ = BASE1 - carry.silow;
  983.         carry.silow = carry.sihigh;
  984.     }
  985.  
  986.     /*
  987.      * Now finally recompute the size of the result.
  988.      */
  989.     len = modlen;
  990.     hd = &res->v[len - 1];
  991.     while ((*hd == 0) && (len > 1)) {
  992.         hd--;
  993.         len--;
  994.     }
  995.     res->len = len;
  996. }
  997.  
  998.  
  999. /*
  1000.  * Square a number in REDC format producing a result also in REDC format.
  1001.  */
  1002. void
  1003. zredcsquare(rp, z1, res)
  1004.     REDC *rp;        /* REDC information */
  1005.     ZVALUE z1;        /* REDC number to be squared */
  1006.     ZVALUE *res;        /* resulting REDC number */
  1007. {
  1008.     ZVALUE tmp;
  1009.  
  1010.     if (isneg(z1))
  1011.         error("Negative number in zredcsquare");
  1012.     if (iszero(z1)) {
  1013.         *res = _zero_;
  1014.         return;
  1015.     }
  1016.     if ((z1.len == rp->one.len) && (z1.v[0] == rp->one.v[0]) &&
  1017.         (zcmp(z1, rp->one) == 0)) {
  1018.             zcopy(z1, res);
  1019.             return;
  1020.     }
  1021.  
  1022.     /*
  1023.      * If the modulus is small enough, then call the multiply
  1024.      * routine to produce the result.  Otherwise call the O(N^1.585)
  1025.      * routines to get the answer.
  1026.      */
  1027.     if (rp->mod.len < _redc2_) {
  1028.         zredcmul(rp, z1, z1, res);
  1029.         return;
  1030.     }
  1031.     zsquare(z1, &tmp);
  1032.     zredcdecode(rp, tmp, res);
  1033.     freeh(tmp.v);
  1034. }
  1035.  
  1036.  
  1037. /*
  1038.  * Compute the result of raising a REDC format number to a power.
  1039.  * The result is within the range 0 to the modulus - 1.
  1040.  * This calculates the result by examining the power POWBITS bits at a time,
  1041.  * using a small table of POWNUMS low powers to calculate powers for those bits,
  1042.  * and repeated squaring and multiplying by the partial powers to generate
  1043.  * the complete power.
  1044.  */
  1045. void
  1046. zredcpower(rp, z1, z2, res)
  1047.     REDC *rp;        /* REDC information */
  1048.     ZVALUE z1;        /* REDC number to be raised */
  1049.     ZVALUE z2;        /* normal number to raise number to */
  1050.     ZVALUE *res;        /* result */
  1051. {
  1052.     HALF *hp;        /* pointer to current word of the power */
  1053.     ZVALUE *pp;        /* pointer to low power table */
  1054.     ZVALUE ans, temp;    /* calculation values */
  1055.     ZVALUE modpow;        /* current small power */
  1056.     ZVALUE lowpowers[POWNUMS];    /* low powers */
  1057.     int curshift;        /* shift value for word of power */
  1058.     HALF curhalf;        /* current word of power */
  1059.     unsigned int curpow;    /* current low power */
  1060.     unsigned int curbit;    /* current bit of low power */
  1061.     int i;
  1062.  
  1063.     if (isneg(z1))
  1064.         error("Negative number in zredcpower");
  1065.     if (isneg(z2))
  1066.         error("Negative power in zredcpower");
  1067.  
  1068.     /*
  1069.      * Check for zero or the REDC format for one.
  1070.      */
  1071.     if (iszero(z1) || isunit(rp->mod)) {
  1072.         *res = _zero_;
  1073.         return;
  1074.     }
  1075.     if (zcmp(z1, rp->one) == 0) {
  1076.         zcopy(rp->one, res);
  1077.         return;
  1078.     }
  1079.  
  1080.     /*
  1081.      * See if the number being raised is the REDC format for -1.
  1082.      * If so, then the answer is the REDC format for one or minus one.
  1083.      * To do this check, calculate the REDC format for -1.
  1084.      */
  1085.     if (((HALF)(z1.v[0] + rp->one.v[0])) == rp->mod.v[0]) {
  1086.         zsub(rp->mod, rp->one, &temp);
  1087.         if (zcmp(z1, temp) == 0) {
  1088.             if (isodd(z2)) {
  1089.                 *res = temp;
  1090.                 return;
  1091.             }
  1092.             freeh(temp.v);
  1093.             zcopy(rp->one, res);
  1094.             return;
  1095.         }
  1096.         freeh(temp.v);
  1097.     }
  1098.  
  1099.     for (pp = &lowpowers[2]; pp < &lowpowers[POWNUMS]; pp++)
  1100.         pp->len = 0;
  1101.     zcopy(rp->one, &lowpowers[0]);
  1102.     zcopy(z1, &lowpowers[1]);
  1103.     zcopy(rp->one, &ans);
  1104.  
  1105.     hp = &z2.v[z2.len - 1];
  1106.     curhalf = *hp;
  1107.     curshift = BASEB - POWBITS;
  1108.     while (curshift && ((curhalf >> curshift) == 0))
  1109.         curshift -= POWBITS;
  1110.  
  1111.     /*
  1112.      * Calculate the result by examining the power POWBITS bits at a time,
  1113.      * and use the table of low powers at each iteration.
  1114.      */
  1115.     for (;;) {
  1116.         curpow = (curhalf >> curshift) & (POWNUMS - 1);
  1117.         pp = &lowpowers[curpow];
  1118.  
  1119.         /*
  1120.          * If the small power is not yet saved in the table, then
  1121.          * calculate it and remember it in the table for future use.
  1122.          */
  1123.         if (pp->len == 0) {
  1124.             if (curpow & 0x1)
  1125.                 zcopy(z1, &modpow);
  1126.             else
  1127.                 zcopy(rp->one, &modpow);
  1128.  
  1129.             for (curbit = 0x2; curbit <= curpow; curbit *= 2) {
  1130.                 pp = &lowpowers[curbit];
  1131.                 if (pp->len == 0)
  1132.                     zredcsquare(rp, lowpowers[curbit/2],
  1133.                         pp);
  1134.                 if (curbit & curpow) {
  1135.                     zredcmul(rp, *pp, modpow, &temp);
  1136.                     freeh(modpow.v);
  1137.                     modpow = temp;
  1138.                 }
  1139.             }
  1140.             pp = &lowpowers[curpow];
  1141.             *pp = modpow;
  1142.         }
  1143.  
  1144.         /*
  1145.          * If the power is nonzero, then accumulate the small power
  1146.          * into the result.
  1147.          */
  1148.         if (curpow) {
  1149.             zredcmul(rp, ans, *pp, &temp);
  1150.             freeh(ans.v);
  1151.             ans = temp;
  1152.         }
  1153.  
  1154.         /*
  1155.          * Select the next POWBITS bits of the power, if there is
  1156.          * any more to generate.
  1157.          */
  1158.         curshift -= POWBITS;
  1159.         if (curshift < 0) {
  1160.             if (hp-- == z2.v)
  1161.                 break;
  1162.             curhalf = *hp;
  1163.             curshift = BASEB - POWBITS;
  1164.         }
  1165.  
  1166.         /*
  1167.          * Square the result POWBITS times to make room for the next
  1168.          * chunk of bits.
  1169.          */
  1170.         for (i = 0; i < POWBITS; i++) {
  1171.             zredcsquare(rp, ans, &temp);
  1172.             freeh(ans.v);
  1173.             ans = temp;
  1174.         }
  1175.     }
  1176.  
  1177.     for (pp = lowpowers; pp < &lowpowers[POWNUMS]; pp++) {
  1178.         if (pp->len)
  1179.             freeh(pp->v);
  1180.     }
  1181.     *res = ans;
  1182. }
  1183.  
  1184. /* END CODE */
  1185.