home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1994 November / macformat-018.iso / Utility Spectacular / Utilities / Calc / zfunc.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-10-10  |  33.5 KB  |  1,676 lines  |  [TEXT/????]
  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.  * Extended precision integral arithmetic non-primitive routines
  7.  */
  8.  
  9. #include "xmath.h"
  10.  
  11. static ZVALUE primeprod;        /* product of primes under 100 */
  12. ZVALUE _tenpowers_[32];            /* table of 10^2^n */
  13.  
  14. #if 0
  15. static char *abortmsg = "Calculation aborted";
  16. static char *memmsg = "Not enough memory";
  17. #endif
  18.  
  19.  
  20. /*
  21.  * Compute the factorial of a number.
  22.  */
  23. void
  24. zfact(z, dest)
  25.     ZVALUE z, *dest;
  26. {
  27.     long ptwo;        /* count of powers of two */
  28.     long n;            /* current multiplication value */
  29.     long m;            /* reduced multiplication value */
  30.     long mul;        /* collected value to multiply by */
  31.     ZVALUE res, temp;
  32.  
  33.     if (isneg(z))
  34.         error("Negative argument for factorial");
  35.     if (isbig(z))
  36.         error("Very large factorial");
  37.     n = (istiny(z) ? z1tol(z) : z2tol(z));
  38.     ptwo = 0;
  39.     mul = 1;
  40.     res = _one_;
  41.     /*
  42.      * Multiply numbers together, but squeeze out all powers of two.
  43.      * We will put them back in at the end.  Also collect multiple
  44.      * numbers together until there is a risk of overflow.
  45.      */
  46.     for (; n > 1; n--) {
  47.         for (m = n; ((m & 0x1) == 0); m >>= 1)
  48.             ptwo++;
  49.         mul *= m;
  50.         if (mul < BASE1/2)
  51.             continue;
  52.         zmuli(res, mul, &temp);
  53.         freeh(res.v);
  54.         res = temp;
  55.         mul = 1;
  56.     }
  57.     /*
  58.      * Multiply by the remaining value, then scale result by
  59.      * the proper power of two.
  60.      */
  61.     if (mul > 1) {
  62.         zmuli(res, mul, &temp);
  63.         freeh(res.v);
  64.         res = temp;
  65.     }
  66.     zshift(res, ptwo, &temp);
  67.     freeh(res.v);
  68.     *dest = temp;
  69. }
  70.  
  71.  
  72. /*
  73.  * Compute the product of the primes up to the specified number.
  74.  */
  75. void
  76. zpfact(z, dest)
  77.     ZVALUE z, *dest;
  78. {
  79.     long n;            /* limiting number to multiply by */
  80.     long p;            /* current prime */
  81.     long i;            /* test value */
  82.     long mul;        /* collected value to multiply by */
  83.     ZVALUE res, temp;
  84.  
  85.     if (isneg(z))
  86.         error("Negative argument for factorial");
  87.     if (isbig(z))
  88.         error("Very large factorial");
  89.     n = (istiny(z) ? z1tol(z) : z2tol(z));
  90.     /*
  91.      * Multiply by the primes in order, collecting multiple numbers
  92.      * together until there is a change of overflow.
  93.      */
  94.     mul = 1 + (n > 1);
  95.     res = _one_;
  96.     for (p = 3; p <= n; p += 2) {
  97.         for (i = 3; (i * i) <= p; i += 2) {
  98.             if ((p % i) == 0)
  99.                 goto next;
  100.         }
  101.         mul *= p;
  102.         if (mul < BASE1/2)
  103.             continue;
  104.         zmuli(res, mul, &temp);
  105.         freeh(res.v);
  106.         res = temp;
  107.         mul = 1;
  108. next: ;
  109.     }
  110.     /*
  111.      * Multiply by the final value if any.
  112.      */
  113.     if (mul > 1) {
  114.         zmuli(res, mul, &temp);
  115.         freeh(res.v);
  116.         res = temp;
  117.     }
  118.     *dest = res;
  119. }
  120.  
  121.  
  122. /*
  123.  * Compute the least common multiple of all the numbers up to the
  124.  * specified number.
  125.  */
  126. void
  127. zlcmfact(z, dest)
  128.     ZVALUE z, *dest;
  129. {
  130.     long n;            /* limiting number to multiply by */
  131.     long p;            /* current prime */
  132.     long pp;        /* power of prime */
  133.     long i;            /* test value */
  134.     ZVALUE res, temp;
  135.  
  136.     if (isneg(z) || iszero(z))
  137.         error("Non-positive argument for lcmfact");
  138.     if (isbig(z))
  139.         error("Very large lcmfact");
  140.     n = (istiny(z) ? z1tol(z) : z2tol(z));
  141.     /*
  142.      * Multiply by powers of the necessary odd primes in order.
  143.      * The power for each prime is the highest one which is not
  144.      * more than the specified number.
  145.      */
  146.     res = _one_;
  147.     for (p = 3; p <= n; p += 2) {
  148.         for (i = 3; (i * i) <= p; i += 2) {
  149.             if ((p % i) == 0)
  150.                 goto next;
  151.         }
  152.         i = p;
  153.         while (i <= n) {
  154.             pp = i;
  155.             i *= p;
  156.         }
  157.         zmuli(res, pp, &temp);
  158.         freeh(res.v);
  159.         res = temp;
  160. next: ;
  161.     }
  162.     /*
  163.      * Finish by scaling by the necessary power of two.
  164.      */
  165.     zshift(res, zhighbit(z), dest);
  166.     freeh(res.v);
  167. }
  168.  
  169.  
  170. /*
  171.  * Compute the permuation function  M! / (M - N)!.
  172.  */
  173. void
  174. zperm(z1, z2, res)
  175.     ZVALUE z1, z2, *res;
  176. {
  177.     long count;
  178.     ZVALUE cur, tmp, ans;
  179.  
  180.     if (isneg(z1) || isneg(z2))
  181.         error("Negative argument for permutation");
  182.     if (zrel(z1, z2) < 0)
  183.         error("Second arg larger than first in permutation");
  184.     if (isbig(z2))
  185.         error("Very large permutation");
  186.     count = (istiny(z2) ? z1tol(z2) : z2tol(z2));
  187.     zcopy(z1, &ans);
  188.     zsub(z1, _one_, &cur);
  189.     while (--count > 0) {
  190.         zmul(ans, cur, &tmp);
  191.         freeh(ans.v);
  192.         ans = tmp;
  193.         zsub(cur, _one_, &tmp);
  194.         freeh(cur.v);
  195.         cur = tmp;
  196.     }
  197.     freeh(cur.v);
  198.     *res = ans;
  199. }
  200.  
  201.  
  202. /*
  203.  * Compute the combinatorial function  M! / ( N! * (M - N)! ).
  204.  */
  205. void
  206. zcomb(z1, z2, res)
  207.     ZVALUE z1, z2, *res;
  208. {
  209.     ZVALUE ans;
  210.     ZVALUE mul, div, temp;
  211.     FULL count, i;
  212.     HALF dh[2];
  213.  
  214.     if (isneg(z1) || isneg(z2))
  215.         error("Negative argument for combinatorial");
  216.     zsub(z1, z2, &temp);
  217.     if (isneg(temp)) {
  218.         freeh(temp.v);
  219.         error("Second arg larger than first for combinatorial");
  220.     }
  221.     if (isbig(z2) && isbig(temp)) {
  222.         freeh(temp.v);
  223.         error("Very large combinatorial");
  224.     }
  225.     count = (istiny(z2) ? z1tol(z2) : z2tol(z2));
  226.     i = (istiny(temp) ? z1tol(temp) : z2tol(temp));
  227.     if (isbig(z2) || (!isbig(temp) && (i < count)))
  228.         count = i;
  229.     freeh(temp.v);
  230.     mul = z1;
  231.     div.sign = 0;
  232.     div.v = dh;
  233.     ans = _one_;
  234.     for (i = 1; i <= count; i++) {
  235.         dh[0] = i & BASE1;
  236.         dh[1] = i / BASE;
  237.         div.len = 1 + (dh[1] != 0);
  238.         zmul(ans, mul, &temp);
  239.         freeh(ans.v);
  240.         zquo(temp, div, &ans);
  241.         freeh(temp.v);
  242.         zsub(mul, _one_, &temp);
  243.         if (mul.v != z1.v)
  244.             freeh(mul.v);
  245.         mul = temp;
  246.     }
  247.     if (mul.v != z1.v)
  248.         freeh(mul.v);
  249.     *res = ans;
  250. }
  251.  
  252.  
  253. /*
  254.  * Perform a probabilistic primality test (algorithm P in Knuth).
  255.  * Returns FALSE if definitely not prime, or TRUE if probably prime.
  256.  * Count determines how many times to check for primality.
  257.  * The chance of a non-prime passing this test is less than (1/4)^count.
  258.  * For example, a count of 100 fails for only 1 in 10^60 numbers.
  259.  */
  260. BOOL
  261. zprimetest(z, count)
  262.     ZVALUE z;        /* number to test for primeness */
  263.     long count;
  264. {
  265.     long ij, ik, ix;
  266.     ZVALUE zm1, z1, z2, z3, ztmp;
  267.     HALF val[2];
  268.  
  269.     z.sign = 0;
  270.     if (iseven(z))        /* if even, not prime if not 2 */
  271.         return (istwo(z) != 0);
  272.     /*
  273.      * See if the number is small, and is either a small prime,
  274.      * or is divisible by a small prime.
  275.      */
  276.     if (istiny(z) && (*z.v <= (HALF)(101*101-1))) {
  277.         ix = *z.v;
  278.         for (ik = 3; (ik <= 97) && ((ik * ik) <= ix); ik += 2)
  279.             if ((ix % ik) == 0)
  280.                 return FALSE;
  281.         return TRUE;
  282.     }
  283.     /*
  284.      * See if the number is divisible by one of the primes 3, 5,
  285.      * 7, 11, or 13.  This is a very easy check.
  286.      */
  287.     ij = zmodi(z, 15015L);
  288.     if (!(ij % 3) || !(ij % 5) || !(ij % 7) || !(ij % 11) || !(ij % 13))
  289.         return FALSE;
  290.     /*
  291.      * Check the gcd of the number and the product of more of the first
  292.      * few odd primes.  We must build the prime product on the first call.
  293.      */
  294.     ztmp.sign = 0;
  295.     ztmp.len = 1;
  296.     ztmp.v = val;
  297.     if (primeprod.len == 0) {
  298.         val[0] = 101;
  299.         zpfact(ztmp, &primeprod);
  300.     }
  301.     zgcd(z, primeprod, &z1);
  302.     if (!isunit(z1)) {
  303.         freeh(z1.v);
  304.         return FALSE;
  305.     }
  306.     freeh(z1.v);
  307.     /*
  308.      * Not divisible by a small prime, so onward with the real test.
  309.      * Make sure the count is limited by the number of odd numbers between
  310.      * three and the number being tested.
  311.      */
  312.     ix = ((istiny(z) ? z1tol(z) : z2tol(z) - 3) / 2);
  313.     if (count > ix) count = ix;
  314.     zsub(z, _one_, &zm1);
  315.     ik = zlowbit(zm1);
  316.     zshift(zm1, -ik, &z1);
  317.     /*
  318.      * Loop over various "random" numbers, testing each one.
  319.      * These numbers are the odd numbers starting from three.
  320.      */
  321.     for (ix = 0; ix < count; ix++) {
  322.         val[0] = (ix * 2) + 3;
  323.         ij = 0;
  324.         zpowermod(ztmp, z1, z, &z3);
  325.         for (;;) {
  326.             if (isone(z3)) {
  327.                 if (ij)    /* number is definitely not prime */
  328.                     goto notprime;
  329.                 break;
  330.             }
  331.             if (zcmp(z3, zm1) == 0)
  332.                 break;
  333.             if (++ij >= ik)
  334.                 goto notprime;    /* number is definitely not prime */
  335.             zsquare(z3, &z2);
  336.             freeh(z3.v);
  337.             zmod(z2, z, &z3);
  338.             freeh(z2.v);
  339.         }
  340.         freeh(z3.v);
  341.     }
  342.     freeh(zm1.v);
  343.     freeh(z1.v);
  344.     return TRUE;    /* number might be prime */
  345.  
  346. notprime:
  347.     freeh(z3.v);
  348.     freeh(zm1.v);
  349.     freeh(z1.v);
  350.     return FALSE;
  351. }
  352.  
  353.  
  354. /*
  355.  * Compute the Jacobi function (p / q) for odd q.
  356.  * If q is prime then the result is:
  357.  *    1 if p == x^2 (mod q) for some x.
  358.  *    -1 otherwise.
  359.  * If q is not prime, then the result is not meaningful if it is 1.
  360.  * This function returns 0 if q is even or q < 0.
  361.  */
  362. FLAG
  363. zjacobi(z1, z2)
  364.     ZVALUE z1, z2;
  365. {
  366.     ZVALUE p, q, tmp;
  367.     long lowbit;
  368.     int val;
  369.  
  370.     if (iseven(z2) || isneg(z2))
  371.         return 0;
  372.     val = 1;
  373.     if (iszero(z1) || isone(z1))
  374.         return val;
  375.     if (isunit(z1)) {
  376.         if ((*z2.v - 1) & 0x2)
  377.             val = -val;
  378.         return val;
  379.     }
  380.     zcopy(z1, &p);
  381.     zcopy(z2, &q);
  382.     for (;;) {
  383.         zmod(p, q, &tmp);
  384.         freeh(p.v);
  385.         p = tmp;
  386.         if (iszero(p)) {
  387.             freeh(p.v);
  388.             p = _one_;
  389.         }
  390.         if (iseven(p)) {
  391.             lowbit = zlowbit(p);
  392.             zshift(p, -lowbit, &tmp);
  393.             freeh(p.v);
  394.             p = tmp;
  395.             if ((lowbit & 1) && (((*q.v & 0x7) == 3) || ((*q.v & 0x7) == 5)))
  396.                 val = -val;
  397.         }
  398.         if (isunit(p)) {
  399.             freeh(p.v);
  400.             freeh(q.v);
  401.             return val;
  402.         }
  403.         if ((*p.v & *q.v & 0x3) == 3)
  404.             val = -val;
  405.         tmp = q;
  406.         q = p;
  407.         p = tmp;
  408.     }
  409. }
  410.  
  411.  
  412. /*
  413.  * Return the Fibonacci number F(n).
  414.  * This is evaluated by recursively using the formulas:
  415.  *    F(2N+1) = F(N+1)^2 + F(N)^2
  416.  * and
  417.  *    F(2N) = F(N+1)^2 - F(N-1)^2
  418.  */
  419. void
  420. zfib(z, res)
  421.     ZVALUE z, *res;
  422. {
  423.     unsigned long i;
  424.     long n;
  425.     int sign;
  426.     ZVALUE fnm1, fn, fnp1;        /* consecutive fibonacci values */
  427.     ZVALUE t1, t2, t3;
  428.  
  429.     if (isbig(z))
  430.         error("Very large Fibonacci number");
  431.     n = (istiny(z) ? z1tol(z) : z2tol(z));
  432.     if (n == 0) {
  433.         *res = _zero_;
  434.         return;
  435.     }
  436.     sign = z.sign && ((n & 0x1) == 0);
  437.     if (n <= 2) {
  438.         *res = _one_;
  439.         res->sign = (BOOL)sign;
  440.         return;
  441.     }
  442.     i = TOPFULL;
  443.     while ((i & n) == 0)
  444.         i >>= 1L;
  445.     i >>= 1L;
  446.     fnm1 = _zero_;
  447.     fn = _one_;
  448.     fnp1 = _one_;
  449.     while (i) {
  450.         zsquare(fnm1, &t1);
  451.         zsquare(fn, &t2);
  452.         zsquare(fnp1, &t3);
  453.         freeh(fnm1.v);
  454.         freeh(fn.v);
  455.         freeh(fnp1.v);
  456.         zadd(t2, t3, &fnp1);
  457.         zsub(t3, t1, &fn);
  458.         freeh(t1.v);
  459.         freeh(t2.v);
  460.         freeh(t3.v);
  461.         if (i & n) {
  462.             fnm1 = fn;
  463.             fn = fnp1;
  464.             zadd(fnm1, fn, &fnp1);
  465.         } else
  466.             zsub(fnp1, fn, &fnm1);
  467.         i >>= 1L;
  468.     }
  469.     freeh(fnm1.v);
  470.     freeh(fnp1.v);
  471.     *res = fn;
  472.     res->sign = (BOOL)sign;
  473. }
  474.  
  475.  
  476. /*
  477.  * Compute the result of raising one number to the power of another
  478.  * The second number is assumed to be non-negative.
  479.  * It cannot be too large except for trivial cases.
  480.  */
  481. void
  482. zpowi(z1, z2, res)
  483.     ZVALUE z1, z2, *res;
  484. {
  485.     int sign;        /* final sign of number */
  486.     unsigned long power;    /* power to raise to */
  487.     unsigned long bit;    /* current bit value */
  488.     long twos;        /* count of times 2 is in result */
  489.     ZVALUE ans, temp;
  490.  
  491.     sign = (z1.sign && isodd(z2));
  492.     z1.sign = 0;
  493.     z2.sign = 0;
  494.     if (iszero(z2)) {    /* number raised to power 0 */
  495.         if (iszero(z1))
  496.             error("Zero raised to zero power");
  497.         *res = _one_;
  498.         return;
  499.     }
  500.     if (isleone(z1)) {    /* 0, 1, or -1 raised to a power */
  501.         ans = _one_;
  502.         ans.sign = (BOOL)sign;
  503.         if (*z1.v == 0)
  504.             ans = _zero_;
  505.         *res = ans;
  506.         return;
  507.     }
  508.     if (isbig(z2))
  509.         error("Raising to very large power");
  510.     power = (istiny(z2) ? z1tol(z2) : z2tol(z2));
  511.     if (istwo(z1)) {    /* two raised to a power */
  512.         zbitvalue((long) power, res);
  513.         return;
  514.     }
  515.     /*
  516.      * See if this is a power of ten
  517.      */
  518.     if (istiny(z1) && (*z1.v == 10)) {
  519.         ztenpow((long) power, res);
  520.         res->sign = (BOOL)sign;
  521.         return;
  522.     }
  523.     /*
  524.      * Handle low powers specially
  525.      */
  526.     if (power <= 4) {
  527.         switch ((int) power) {
  528.             case 1:
  529.                 ans.len = z1.len;
  530.                 ans.v = alloc(ans.len);
  531.                 copyval(z1, ans);
  532.                 ans.sign = (BOOL)sign;
  533.                 *res = ans;
  534.                 return;
  535.             case 2:
  536.                 zsquare(z1, res);
  537.                 return;
  538.             case 3:
  539.                 zsquare(z1, &temp);
  540.                 zmul(z1, temp, res);
  541.                 freeh(temp.v);
  542.                 res->sign = (BOOL)sign;
  543.                 return;
  544.             case 4:
  545.                 zsquare(z1, &temp);
  546.                 zsquare(temp, res);
  547.                 freeh(temp.v);
  548.                 return;
  549.         }
  550.     }
  551.     /*
  552.      * Shift out all powers of twos so the multiplies are smaller.
  553.      * We will shift back the right amount when done.
  554.      */
  555.     twos = 0;
  556.     if (iseven(z1)) {
  557.         twos = zlowbit(z1);
  558.         ans.v = alloc(z1.len);
  559.         ans.len = z1.len;
  560.         copyval(z1, ans);
  561.         shiftr(ans, twos);
  562.         trim(&ans);
  563.         z1 = ans;
  564.         twos *= power;
  565.     }
  566.     /*
  567.      * Compute the power by squaring and multiplying.
  568.      * This uses the left to right method of power raising.
  569.      */
  570.     bit = TOPFULL;
  571.     while ((bit & power) == 0)
  572.         bit >>= 1L;
  573.     bit >>= 1L;
  574.     zsquare(z1, &ans);
  575.     if (bit & power) {
  576.         zmul(ans, z1, &temp);
  577.         freeh(ans.v);
  578.         ans = temp;
  579.     }
  580.     bit >>= 1L;
  581.     while (bit) {
  582.         zsquare(ans, &temp);
  583.         freeh(ans.v);
  584.         ans = temp;
  585.         if (bit & power) {
  586.             zmul(ans, z1, &temp);
  587.             freeh(ans.v);
  588.             ans = temp;
  589.         }
  590.         bit >>= 1L;
  591.     }
  592.     /*
  593.      * Scale back up by proper power of two
  594.      */
  595.     if (twos) {
  596.         zshift(ans, twos, &temp);
  597.         freeh(ans.v);
  598.         ans = temp;
  599.         freeh(z1.v);
  600.     }
  601.     ans.sign = (BOOL)sign;
  602.     *res = ans;
  603. }
  604.  
  605.  
  606. /*
  607.  * Compute ten to the specified power
  608.  * This saves some work since the squares of ten are saved.
  609.  */
  610. void
  611. ztenpow(power, res)
  612.     long power;
  613.     ZVALUE *res;
  614. {
  615.     long i;
  616.     ZVALUE ans;
  617.     ZVALUE temp;
  618.  
  619.     if (power <= 0) {
  620.         *res = _one_;
  621.         return;
  622.     }
  623.     ans = _one_;
  624.     _tenpowers_[0] = _ten_;
  625.     for (i = 0; power; i++) {
  626.         if (_tenpowers_[i].len == 0)
  627.             zsquare(_tenpowers_[i-1], &_tenpowers_[i]);
  628.         if (power & 0x1) {
  629.             zmul(ans, _tenpowers_[i], &temp);
  630.             freeh(ans.v);
  631.             ans = temp;
  632.         }
  633.         power /= 2;
  634.     }
  635.     *res = ans;
  636. }
  637.  
  638.  
  639. /*
  640.  * Calculate modular inverse suppressing unnecessary divisions.
  641.  * This is based on the Euclidian algorithm for large numbers.
  642.  * (Algorithm X from Knuth Vol 2, section 4.5.2. and exercise 17)
  643.  * Returns TRUE if there is no solution because the numbers
  644.  * are not relatively prime.
  645.  */
  646. BOOL
  647. zmodinv(u, v, res)
  648.     ZVALUE u, v;
  649.     ZVALUE *res;
  650. {
  651.     FULL    q1, q2, ui3, vi3, uh, vh, A, B, C, D, T;
  652.     ZVALUE    u2, u3, v2, v3, qz, tmp1, tmp2, tmp3;
  653.  
  654.     if (isneg(u) || isneg(v) || (zrel(u, v) >= 0))
  655.         zmod(u, v, &v3);
  656.     else
  657.         zcopy(u, &v3);
  658.     zcopy(v, &u3);
  659.     u2 = _zero_;
  660.     v2 = _one_;
  661.  
  662.     /*
  663.      * Loop here while the size of the numbers remain above
  664.      * the size of a FULL.  Throughout this loop u3 >= v3.
  665.      */
  666.     while ((u3.len > 1) && !iszero(v3)) {
  667.         uh = (((FULL) u3.v[u3.len - 1]) << BASEB) + u3.v[u3.len - 2];
  668.         vh = 0;
  669.         if ((v3.len + 1) >= u3.len)
  670.             vh = v3.v[v3.len - 1];
  671.         if (v3.len == u3.len)
  672.             vh = (vh << BASEB) + v3.v[v3.len - 2];
  673.         A = 1;
  674.         B = 0;
  675.         C = 0;
  676.         D = 1;
  677.  
  678.         /*
  679.          * Calculate successive quotients of the continued fraction
  680.          * expansion using only single precision arithmetic until
  681.          * greater precision is required.
  682.          */
  683.         while ((vh + C) && (vh + D)) {
  684.             q1 = (uh + A) / (vh + C);
  685.             q2 = (uh + B) / (vh + D);
  686.             if (q1 != q2)
  687.                 break;
  688.             T = A - q1 * C;
  689.             A = C;
  690.             C = T;
  691.             T = B - q1 * D;
  692.             B = D;
  693.             D = T;
  694.             T = uh - q1 * vh;
  695.             uh = vh;
  696.             vh = T;
  697.         }
  698.     
  699.         /*
  700.          * If B is zero, then we made no progress because
  701.          * the calculation requires a very large quotient.
  702.          * So we must do this step of the calculation in
  703.          * full precision
  704.          */
  705.         if (B == 0) {
  706.             zquo(u3, v3, &qz);
  707.             zmul(qz, v2, &tmp1);
  708.             zsub(u2, tmp1, &tmp2);
  709.             freeh(tmp1.v);
  710.             freeh(u2.v);
  711.             u2 = v2;
  712.             v2 = tmp2;
  713.             zmul(qz, v3, &tmp1);
  714.             zsub(u3, tmp1, &tmp2);
  715.             freeh(tmp1.v);
  716.             freeh(u3.v);
  717.             u3 = v3;
  718.             v3 = tmp2;
  719.             freeh(qz.v);
  720.             continue;
  721.         }
  722.         /*
  723.          * Apply the calculated A,B,C,D numbers to the current
  724.          * values to update them as if the full precision
  725.          * calculations had been carried out.
  726.          */
  727.         zmuli(u2, (long) A, &tmp1);
  728.         zmuli(v2, (long) B, &tmp2);
  729.         zadd(tmp1, tmp2, &tmp3);
  730.         freeh(tmp1.v);
  731.         freeh(tmp2.v);
  732.         zmuli(u2, (long) C, &tmp1);
  733.         zmuli(v2, (long) D, &tmp2);
  734.         freeh(u2.v);
  735.         freeh(v2.v);
  736.         u2 = tmp3;
  737.         zadd(tmp1, tmp2, &v2);
  738.         freeh(tmp1.v);
  739.         freeh(tmp2.v);
  740.         zmuli(u3, (long) A, &tmp1);
  741.         zmuli(v3, (long) B, &tmp2);
  742.         zadd(tmp1, tmp2, &tmp3);
  743.         freeh(tmp1.v);
  744.         freeh(tmp2.v);
  745.         zmuli(u3, (long) C, &tmp1);
  746.         zmuli(v3, (long) D, &tmp2);
  747.         freeh(u3.v);
  748.         freeh(v3.v);
  749.         u3 = tmp3;
  750.         zadd(tmp1, tmp2, &v3);
  751.         freeh(tmp1.v);
  752.         freeh(tmp2.v);
  753.     }
  754.  
  755.     /*
  756.      * Here when the remaining numbers become single precision in size.
  757.      * Finish the procedure using single precision calculations.
  758.      */
  759.     if (iszero(v3) && !isone(u3)) {
  760.         freeh(u3.v);
  761.         freeh(v3.v);
  762.         freeh(u2.v);
  763.         freeh(v2.v);
  764.         return TRUE;
  765.     }
  766.     ui3 = (istiny(u3) ? z1tol(u3) : z2tol(u3));
  767.     vi3 = (istiny(v3) ? z1tol(v3) : z2tol(v3));
  768.     freeh(u3.v);
  769.     freeh(v3.v);
  770.     while (vi3) {
  771.         q1 = ui3 / vi3;
  772.         zmuli(v2, (long) q1, &tmp1);
  773.         zsub(u2, tmp1, &tmp2);
  774.         freeh(tmp1.v);
  775.         freeh(u2.v);
  776.         u2 = v2;
  777.         v2 = tmp2;
  778.         q2 = ui3 - q1 * vi3;
  779.         ui3 = vi3;
  780.         vi3 = q2;
  781.     }
  782.     freeh(v2.v);
  783.     if (ui3 != 1) {
  784.         freeh(u2.v);
  785.         return TRUE;
  786.     }
  787.     if (isneg(u2)) {
  788.         zadd(v, u2, res);
  789.         freeh(u2.v);
  790.         return FALSE;
  791.     }
  792.     *res = u2;
  793.     return FALSE;
  794. }
  795.  
  796.  
  797. #if 0
  798. /*
  799.  * Approximate the quotient of two integers by another set of smaller
  800.  * integers.  This uses continued fractions to determine the smaller set.
  801.  */
  802. void
  803. zapprox(z1, z2, res1, res2)
  804.     ZVALUE z1, z2, *res1, *res2;
  805. {
  806.     int sign;
  807.     ZVALUE u1, v1, u3, v3, q, t1, t2, t3;
  808.  
  809.     sign = ((z1.sign != 0) ^ (z2.sign != 0));
  810.     z1.sign = 0;
  811.     z2.sign = 0;
  812.     v3 = z2;
  813.     u3 = z1;
  814.     u1 = _one_;
  815.     v1 = _zero_;
  816.     while (!iszero(v3)) {
  817.         zdiv(u3, v3, &q, &t1);
  818.         zmul(v1, q, &t2);
  819.         zsub(u1, t2, &t3);
  820.         freeh(q.v);
  821.         freeh(t2.v);
  822.         freeh(u1.v);
  823.         if ((u3.v != z1.v) && (u3.v != z2.v))
  824.             freeh(u3.v);
  825.         u1 = v1;
  826.         u3 = v3;
  827.         v1 = t3;
  828.         v3 = t1;
  829.     }
  830.     if (!isunit(u3))
  831.         error("Non-relativly prime numbers for approx");
  832.     if ((u3.v != z1.v) && (u3.v != z2.v))
  833.         freeh(u3.v);
  834.     if ((v3.v != z1.v) && (v3.v != z2.v))
  835.         freeh(v3.v);
  836.     freeh(v1.v);
  837.     zmul(u1, z1, &t1);
  838.     zsub(t1, _one_, &t2);
  839.     freeh(t1.v);
  840.     zquo(t2, z2, &t1);
  841.     freeh(t2.v);
  842.     u1.sign = (BOOL)sign;
  843.     t1.sign = 0;
  844.     *res1 = t1;
  845.     *res2 = u1;
  846. }
  847. #endif
  848.  
  849.  
  850. /*
  851.  * Binary gcd algorithm
  852.  * This algorithm taken from Knuth
  853.  */
  854. void
  855. zgcd(z1, z2, res)
  856.     ZVALUE z1, z2, *res;
  857. {
  858.     ZVALUE u, v, t;
  859.     register long j, k, olen, mask;
  860.     register HALF h;
  861.     HALF *oldv1, *oldv2;
  862.  
  863.     /*
  864.      * First see if one number is very much larger than the other.
  865.      * If so, then divide as necessary to get the numbers near each other.
  866.      */
  867.     z1.sign = 0;
  868.     z2.sign = 0;
  869.     oldv1 = z1.v;
  870.     oldv2 = z2.v;
  871.     if (z1.len < z2.len) {
  872.         t = z1;
  873.         z1 = z2;
  874.         z2 = t;
  875.     }
  876.     while ((z1.len > (z2.len + 5)) && !iszero(z2)) {
  877.         zmod(z1, z2, &t);
  878.         if ((z1.v != oldv1) && (z1.v != oldv2))
  879.             freeh(z1.v);
  880.         z1 = z2;
  881.         z2 = t;
  882.     }
  883.     /*
  884.      * Ok, now do the binary method proper
  885.      */
  886.     u.len = z1.len;
  887.     v.len = z2.len;
  888.     u.sign = 0;
  889.     v.sign = 0;
  890.     if (!ztest(z1)) {
  891.         v.v = alloc(v.len);
  892.         copyval(z2, v);
  893.         *res = v;
  894.         goto done;
  895.     }
  896.     if (!ztest(z2)) {
  897.         u.v = alloc(u.len);
  898.         copyval(z1, u);
  899.         *res = u;
  900.         goto done;
  901.     }
  902.     u.v = alloc(u.len);
  903.     v.v = alloc(v.len);
  904.     copyval(z1, u);
  905.     copyval(z2, v);
  906.     k = 0;
  907.     while (u.v[k] == 0 && v.v[k] == 0)
  908.         ++k;
  909.     mask = 01;
  910.     h = u.v[k] | v.v[k];
  911.     k *= BASEB;
  912.     while (!(h & mask)) {
  913.         mask <<= 1;
  914.         k++;
  915.     }
  916.     shiftr(u, k);
  917.     shiftr(v, k);
  918.     trim(&u);
  919.     trim(&v);
  920.     if (isodd(u)) {
  921.         t.v = alloc(v.len);
  922.         t.len = v.len;
  923.         copyval(v, t);
  924.         t.sign = !v.sign;
  925.     } else {
  926.         t.v = alloc(u.len);
  927.         t.len = u.len;
  928.         copyval(u, t);
  929.         t.sign = u.sign;
  930.     }
  931.     while (ztest(t)) {
  932.         j = 0;
  933.         while (!t.v[j])
  934.             ++j;
  935.         mask = 01;
  936.         h = t.v[j];
  937.         j *= BASEB;
  938.         while (!(h & mask)) {
  939.             mask <<= 1;
  940.             j++;
  941.         }
  942.         shiftr(t, j);
  943.         trim(&t);
  944.         if (ztest(t) > 0) {
  945.             freeh(u.v);
  946.             u = t;
  947.         } else {
  948.             freeh(v.v);
  949.             v = t;
  950.             v.sign = !t.sign;
  951.         }
  952.         zsub(u, v, &t);
  953.     }
  954.     freeh(t.v);
  955.     freeh(v.v);
  956.     if (k) {
  957.         olen = u.len;
  958.         u.len += k / BASEB + 1;
  959.         u.v = (HALF *) realloc(u.v, u.len * sizeof(HALF));
  960.         while (olen != u.len)
  961.             u.v[olen++] = 0;
  962.         shiftl(u, k);
  963.     }
  964.     trim(&u);
  965.     *res = u;
  966.  
  967. done:
  968.     if ((z1.v != oldv1) && (z1.v != oldv2))
  969.         freeh(z1.v);
  970.     if ((z2.v != oldv1) && (z2.v != oldv2))
  971.         freeh(z2.v);
  972. }
  973.  
  974.  
  975. /*
  976.  * Compute the lcm of two integers (least common multiple).
  977.  * This is done using the formula:  gcd(a,b) * lcm(a,b) = a * b.
  978.  */
  979. void
  980. zlcm(z1, z2, res)
  981.     ZVALUE z1, z2, *res;
  982. {
  983.     ZVALUE temp1, temp2;
  984.  
  985.     zgcd(z1, z2, &temp1);
  986.     zquo(z1, temp1, &temp2);
  987.     freeh(temp1.v);
  988.     zmul(temp2, z2, res);
  989.     freeh(temp2.v);
  990. }
  991.  
  992.  
  993. /*
  994.  * Return whether or not two numbers are relatively prime to each other.
  995.  */
  996. BOOL
  997. zrelprime(z1, z2)
  998.     ZVALUE z1, z2;            /* numbers to be tested */
  999. {
  1000.     FULL rem1, rem2;        /* remainders */
  1001.     ZVALUE rem;
  1002.     BOOL result;
  1003.  
  1004.     z1.sign = 0;
  1005.     z2.sign = 0;
  1006.     if (iseven(z1) && iseven(z2))    /* false if both even */
  1007.         return FALSE;
  1008.     if (isunit(z1) || isunit(z2))    /* true if either is a unit */
  1009.         return TRUE;
  1010.     if (iszero(z1) || iszero(z2))    /* false if either is zero */
  1011.         return FALSE;
  1012.     if (istwo(z1) || istwo(z2))    /* true if either is two */
  1013.         return TRUE;
  1014.     /*
  1015.      * Try reducing each number by the product of the first few odd primes
  1016.      * to see if any of them are a common factor.
  1017.      */
  1018.     rem1 = zmodi(z1, 3L * 5 * 7 * 11 * 13);
  1019.     rem2 = zmodi(z2, 3L * 5 * 7 * 11 * 13);
  1020.     if (((rem1 % 3) == 0) && ((rem2 % 3) == 0))
  1021.         return FALSE;
  1022.     if (((rem1 % 5) == 0) && ((rem2 % 5) == 0))
  1023.         return FALSE;
  1024.     if (((rem1 % 7) == 0) && ((rem2 % 7) == 0))
  1025.         return FALSE;
  1026.     if (((rem1 % 11) == 0) && ((rem2 % 11) == 0))
  1027.         return FALSE;
  1028.     if (((rem1 % 13) == 0) && ((rem2 % 13) == 0))
  1029.         return FALSE;
  1030.     /*
  1031.      * Try a new batch of primes now
  1032.      */
  1033.     rem1 = zmodi(z1, 17L * 19 * 23);
  1034.     rem2 = zmodi(z2, 17L * 19 * 23);
  1035.     if (((rem1 % 17) == 0) && ((rem2 % 17) == 0))
  1036.         return FALSE;
  1037.     if (((rem1 % 19) == 0) && ((rem2 % 19) == 0))
  1038.         return FALSE;
  1039.     if (((rem1 % 23) == 0) && ((rem2 % 23) == 0))
  1040.         return FALSE;
  1041.     /*
  1042.      * Yuk, we must actually compute the gcd to know the answer
  1043.      */
  1044.     zgcd(z1, z2, &rem);
  1045.     result = isunit(rem);
  1046.     freeh(rem.v);
  1047.     return result;
  1048. }
  1049.  
  1050.  
  1051. /*
  1052.  * Compute the log of one number base another, to the closest integer.
  1053.  * This is the largest integer which when the second number is raised to it,
  1054.  * the resulting value is less than or equal to the first number.
  1055.  * Example:  zlog(123456, 10) = 5.
  1056.  */
  1057. long
  1058. zlog(z1, z2)
  1059.     ZVALUE z1, z2;
  1060. {
  1061.     register ZVALUE *zp;        /* current square */
  1062.     long power;            /* current power */
  1063.     long worth;            /* worth of current square */
  1064.     ZVALUE val;            /* current value of power */
  1065.     ZVALUE temp;            /* temporary */
  1066.     ZVALUE squares[32];        /* table of squares of base */
  1067.  
  1068.     /*
  1069.      * Make sure that the numbers are nonzero and the base is greater than one.
  1070.      */
  1071.     if (isneg(z1) || iszero(z1) || isneg(z2) || isleone(z2))
  1072.         error("Bad arguments for log");
  1073.     /*
  1074.      * Reject trivial cases.
  1075.      */
  1076.     if (z1.len < z2.len)
  1077.         return 0;
  1078.     if ((z1.len == z2.len) && (z1.v[z1.len-1] < z2.v[z2.len-1]))
  1079.         return 0;
  1080.     power = zrel(z1, z2);
  1081.     if (power <= 0)
  1082.         return (power + 1);
  1083.     /*
  1084.      * Handle any power of two special.
  1085.      */
  1086.     if (zisonebit(z2))
  1087.         return (zhighbit(z1) / zlowbit(z2));
  1088.     /*
  1089.      * Handle base 10 special
  1090.      */
  1091.     if ((z2.len == 1) && (*z2.v == 10))
  1092.         return zlog10(z1);
  1093.     /*
  1094.      * Now loop by squaring the base each time, and see whether or
  1095.      * not each successive square is still smaller than the number.
  1096.      */
  1097.     worth = 1;
  1098.     zp = &squares[0];
  1099.     *zp = z2;
  1100.     while (((zp->len * 2) - 1) <= z1.len) {    /* while square not too large */
  1101.         zsquare(*zp, zp + 1);
  1102.         zp++;
  1103.         worth *= 2;
  1104.     }
  1105.     /*
  1106.      * Now back down the squares, and multiply them together to see
  1107.      * exactly how many times the base can be raised by.
  1108.      */
  1109.     val = _one_;
  1110.     power = 0;
  1111.     for (; zp >= squares; zp--, worth /= 2) {
  1112.         if ((val.len + zp->len - 1) <= z1.len) {
  1113.             zmul(val, *zp, &temp);
  1114.             if (zrel(z1, temp) >= 0) {
  1115.                 freeh(val.v);
  1116.                 val = temp;
  1117.                 power += worth;
  1118.             } else
  1119.                 freeh(temp.v);
  1120.         }
  1121.         if (zp != squares)
  1122.             freeh(zp->v);
  1123.     }
  1124.     return power;
  1125. }
  1126.  
  1127.  
  1128. /*
  1129.  * Return the integral log base 10 of a number.
  1130.  */
  1131. long
  1132. zlog10(z)
  1133.     ZVALUE z;
  1134. {
  1135.     register ZVALUE *zp;        /* current square */
  1136.     long power;            /* current power */
  1137.     long worth;            /* worth of current square */
  1138.     ZVALUE val;            /* current value of power */
  1139.     ZVALUE temp;            /* temporary */
  1140.  
  1141.     if (!ispos(z))
  1142.         error("Non-positive number for log10");
  1143.     /*
  1144.      * Loop by squaring the base each time, and see whether or
  1145.      * not each successive square is still smaller than the number.
  1146.      */
  1147.     worth = 1;
  1148.     zp = &_tenpowers_[0];
  1149.     *zp = _ten_;
  1150.     while (((zp->len * 2) - 1) <= z.len) {    /* while square not too large */
  1151.         if (zp[1].len == 0)
  1152.             zsquare(*zp, zp + 1);
  1153.         zp++;
  1154.         worth *= 2;
  1155.     }
  1156.     /*
  1157.      * Now back down the squares, and multiply them together to see
  1158.      * exactly how many times the base can be raised by.
  1159.      */
  1160.     val = _one_;
  1161.     power = 0;
  1162.     for (; zp >= _tenpowers_; zp--, worth /= 2) {
  1163.         if ((val.len + zp->len - 1) <= z.len) {
  1164.             zmul(val, *zp, &temp);
  1165.             if (zrel(z, temp) >= 0) {
  1166.                 freeh(val.v);
  1167.                 val = temp;
  1168.                 power += worth;
  1169.             } else
  1170.                 freeh(temp.v);
  1171.         }
  1172.     }
  1173.     return power;
  1174. }
  1175.  
  1176.  
  1177. /*
  1178.  * Return the number of times that one number will divide another.
  1179.  * This works similarly to zlog, except that divisions must be exact.
  1180.  * For example, zdivcount(540, 3) = 3, since 3^3 divides 540 but 3^4 won't.
  1181.  */
  1182. long
  1183. zdivcount(z1, z2)
  1184.     ZVALUE z1, z2;
  1185. {
  1186.     long count;        /* number of factors removed */
  1187.     ZVALUE tmp;        /* ignored return value */
  1188.  
  1189.     count = zfacrem(z1, z2, &tmp);
  1190.     freeh(tmp.v);
  1191.     return count;
  1192. }
  1193.  
  1194.  
  1195. /*
  1196.  * Remove all occurances of the specified factor from a number.
  1197.  * Also returns the number of factors removed as a function return value.
  1198.  * Example:  zfacrem(540, 3, &x) returns 3 and sets x to 20.
  1199.  */
  1200. long
  1201. zfacrem(z1, z2, rem)
  1202.     ZVALUE z1, z2, *rem;
  1203. {
  1204.     register ZVALUE *zp;        /* current square */
  1205.     long count;            /* total count of divisions */
  1206.     long worth;            /* worth of current square */
  1207.     ZVALUE temp1, temp2, temp3;    /* temporaries */
  1208.     ZVALUE squares[32];        /* table of squares of factor */
  1209.  
  1210.     z1.sign = 0;
  1211.     z2.sign = 0;
  1212.     /*
  1213.      * Make sure factor isn't 0 or 1.
  1214.      */
  1215.     if (isleone(z2))
  1216.         error("Bad argument for facrem");
  1217.     /*
  1218.      * Reject trivial cases.
  1219.      */
  1220.     if ((z1.len < z2.len) || (isodd(z1) && iseven(z2)) ||
  1221.         ((z1.len == z2.len) && (z1.v[z1.len-1] < z2.v[z2.len-1]))) {
  1222.         rem->v = alloc(z1.len);
  1223.         rem->len = z1.len;
  1224.         rem->sign = 0;
  1225.         copyval(z1, *rem);
  1226.         return 0;
  1227.     }
  1228.     /*
  1229.      * Handle any power of two special.
  1230.      */
  1231.     if (zisonebit(z2)) {
  1232.         count = zlowbit(z1) / zlowbit(z2);
  1233.         rem->v = alloc(z1.len);
  1234.         rem->len = z1.len;
  1235.         rem->sign = 0;
  1236.         copyval(z1, *rem);
  1237.         shiftr(*rem, count);
  1238.         trim(rem);
  1239.         return count;
  1240.     }
  1241.     /*
  1242.      * See if the factor goes in even once.
  1243.      */
  1244.     zdiv(z1, z2, &temp1, &temp2);
  1245.     if (!iszero(temp2)) {
  1246.         freeh(temp1.v);
  1247.         freeh(temp2.v);
  1248.         rem->v = alloc(z1.len);
  1249.         rem->len = z1.len;
  1250.         rem->sign = 0;
  1251.         copyval(z1, *rem);
  1252.         return 0;
  1253.     }
  1254.     freeh(temp2.v);
  1255.     z1 = temp1;
  1256.     /*
  1257.      * Now loop by squaring the factor each time, and see whether
  1258.      * or not each successive square will still divide the number.
  1259.      */
  1260.     count = 1;
  1261.     worth = 1;
  1262.     zp = &squares[0];
  1263.     *zp = z2;
  1264.     while (((zp->len * 2) - 1) <= z1.len) {    /* while square not too large */
  1265.         zsquare(*zp, &temp1);
  1266.         zdiv(z1, temp1, &temp2, &temp3);
  1267.         if (!iszero(temp3)) {
  1268.             freeh(temp1.v);
  1269.             freeh(temp2.v);
  1270.             freeh(temp3.v);
  1271.             break;
  1272.         }
  1273.         freeh(temp3.v);
  1274.         freeh(z1.v);
  1275.         z1 = temp2;
  1276.         *++zp = temp1;
  1277.         worth *= 2;
  1278.         count += worth;
  1279.     }
  1280.     /*
  1281.      * Now back down the list of squares, and see if the lower powers
  1282.      * will divide any more times.
  1283.      */
  1284.     for (; zp >= squares; zp--, worth /= 2) {
  1285.         if (zp->len <= z1.len) {
  1286.             zdiv(z1, *zp, &temp1, &temp2);
  1287.             if (iszero(temp2)) {
  1288.                 temp3 = z1;
  1289.                 z1 = temp1;
  1290.                 temp1 = temp3;
  1291.                 count += worth;
  1292.             }
  1293.             freeh(temp1.v);
  1294.             freeh(temp2.v);
  1295.         }
  1296.         if (zp != squares)
  1297.             freeh(zp->v);
  1298.     }
  1299.     *rem = z1;
  1300.     return count;
  1301. }
  1302.  
  1303.  
  1304. /*
  1305.  * Keep dividing a number by the gcd of it with another number until the
  1306.  * result is relatively prime to the second number.
  1307.  */
  1308. void
  1309. zgcdrem(z1, z2, res)
  1310.     ZVALUE z1, z2, *res;
  1311. {
  1312.     ZVALUE tmp1, tmp2;
  1313.  
  1314.     /*
  1315.      * Begin by taking the gcd for the first time.
  1316.      * If the number is already relatively prime, then we are done.
  1317.      */
  1318.     zgcd(z1, z2, &tmp1);
  1319.     if (isunit(tmp1) || iszero(tmp1)) {
  1320.         res->len = z1.len;
  1321.         res->v = alloc(z1.len);
  1322.         res->sign = z1.sign;
  1323.         copyval(z1, *res);
  1324.         return;
  1325.     }
  1326.     zquo(z1, tmp1, &tmp2);
  1327.     z1 = tmp2;
  1328.     z2 = tmp1;
  1329.     /*
  1330.      * Now keep alternately taking the gcd and removing factors until
  1331.      * the gcd becomes one.
  1332.      */
  1333.     while (!isunit(z2)) {
  1334.         (void) zfacrem(z1, z2, &tmp1);
  1335.         freeh(z1.v);
  1336.         z1 = tmp1;
  1337.         zgcd(z1, z2, &tmp1);
  1338.         freeh(z2.v);
  1339.         z2 = tmp1;
  1340.     }
  1341.     *res = z1;
  1342. }
  1343.  
  1344.  
  1345. /*
  1346.  * Find the lowest prime factor of a number if one can be found.
  1347.  * Search is conducted for the first count primes.
  1348.  * Returns 1 if no factor was found.
  1349.  */
  1350. long
  1351. zlowfactor(z, count)
  1352.     ZVALUE z;
  1353.     long count;
  1354. {
  1355.     FULL p, d;
  1356.     ZVALUE div, tmp;
  1357.     HALF divval[2];
  1358.  
  1359.     if ((--count < 0) || iszero(z))
  1360.         return 1;
  1361.     if (iseven(z))
  1362.         return 2;
  1363.     div.sign = 0;
  1364.     div.v = divval;
  1365.     for (p = 3; (count > 0); p += 2) {
  1366.         for (d = 3; (d * d) <= p; d += 2)
  1367.             if ((p % d) == 0)
  1368.                 goto next;
  1369.         divval[0] = (p & BASE1);
  1370.         divval[1] = (p / BASE);
  1371.         div.len = 1 + (p >= BASE);
  1372.         zmod(z, div, &tmp);
  1373.         if (iszero(tmp))
  1374.             return p;
  1375.         freeh(tmp.v);
  1376.         count--;
  1377. next:;
  1378.     }
  1379.     return 1;
  1380. }
  1381.  
  1382.  
  1383. /*
  1384.  * Return the number of digits (base 10) in a number, ignoring the sign.
  1385.  */
  1386. long
  1387. zdigits(z1)
  1388.     ZVALUE z1;
  1389. {
  1390.     long count, val;
  1391.  
  1392.     z1.sign = 0;
  1393.     if (istiny(z1)) {    /* do small numbers ourself */
  1394.         count = 1;
  1395.         val = 10;
  1396.         while (*z1.v >= (HALF)val) {
  1397.             count++;
  1398.             val *= 10;
  1399.         }
  1400.         return count;
  1401.     }
  1402.     return (zlog10(z1) + 1);
  1403. }
  1404.  
  1405.  
  1406. /*
  1407.  * Return the single digit at the specified decimal place of a number,
  1408.  * where 0 means the rightmost digit.  Example:  zdigit(1234, 1) = 3.
  1409.  */
  1410. FLAG
  1411. zdigit(z1, n)
  1412.     ZVALUE z1;
  1413.     long n;
  1414. {
  1415.     ZVALUE tmp1, tmp2;
  1416.     FLAG res;
  1417.  
  1418.     z1.sign = 0;
  1419.     if (iszero(z1) || (n < 0) || (n / BASEDIG >= z1.len))
  1420.         return 0;
  1421.     if (n == 0)
  1422.         return zmodi(z1, 10L);
  1423.     if (n == 1)
  1424.         return zmodi(z1, 100L) / 10;
  1425.     if (n == 2)
  1426.         return zmodi(z1, 1000L) / 100;
  1427.     if (n == 3)
  1428.         return zmodi(z1, 10000L) / 1000;
  1429.     ztenpow(n, &tmp1);
  1430.     zquo(z1, tmp1, &tmp2);
  1431.     res = zmodi(tmp2, 10L);
  1432.     freeh(tmp1.v);
  1433.     freeh(tmp2.v);
  1434.     return res;
  1435. }
  1436.  
  1437.  
  1438. /*
  1439.  * Find the square root of a number.  This is the greatest integer whose
  1440.  * square is less than or equal to the number. Returns TRUE if the
  1441.  * square root is exact.
  1442.  */
  1443. BOOL
  1444. zsqrt(z1, dest)
  1445.     ZVALUE z1, *dest;
  1446. {
  1447.     ZVALUE try, quo, rem, old, temp;
  1448.     FULL iquo, val;
  1449.     long i,j;
  1450.  
  1451.     if (z1.sign)
  1452.         error("Square root of negative number");
  1453.     if (iszero(z1)) {
  1454.         *dest = _zero_;
  1455.         return TRUE;
  1456.     }
  1457.     if ((*z1.v < 4) && istiny(z1)) {
  1458.         *dest = _one_;
  1459.         return (*z1.v == 1);
  1460.     }
  1461.     /*
  1462.      * Pick the square root of the leading one or two digits as a first guess.
  1463.      */
  1464.     val = z1.v[z1.len-1];
  1465.     if ((z1.len & 0x1) == 0)
  1466.         val = (val * BASE) + z1.v[z1.len-2];
  1467.  
  1468.     /*
  1469.      * Find the largest power of 2 that when squared
  1470.      * is <= val > 0.  Avoid multiply overflow by doing 
  1471.      * a careful check at the BASE boundary.
  1472.      */
  1473.     j = 1<<(BASEB+BASEB-2);
  1474.     if (val > j) {
  1475.         iquo = BASE;
  1476.     } else {
  1477.         i = BASEB-1;
  1478.         while (j > val) {
  1479.             --i;
  1480.             j >>= 2;
  1481.         }
  1482.         iquo = bitmask[i];
  1483.     }
  1484.  
  1485.     for (i = 8; i > 0; i--)
  1486.         iquo = (iquo + (val / iquo)) / 2;
  1487.     if (iquo > BASE1)
  1488.         iquo = BASE1;
  1489.     /*
  1490.      * Allocate the numbers to use for the main loop.
  1491.      * The size and high bits of the final result are correctly set here.
  1492.      * Notice that the remainder of the test value is rubbish, but this
  1493.      * is unimportant.
  1494.      */
  1495.     try.sign = 0;
  1496.     try.len = (z1.len + 1) / 2;
  1497.     try.v = alloc(try.len);
  1498.     try.v[try.len-1] = (HALF)iquo;
  1499.     old.sign = 0;
  1500.     old.v = alloc(try.len);
  1501.     old.len = 1;
  1502.     *old.v = 0;
  1503.     /*
  1504.      * Main divide and average loop
  1505.      */
  1506.     for (;;) {
  1507.         zdiv(z1, try, &quo, &rem);
  1508.         i = zrel(try, quo);
  1509.         if ((i == 0) && iszero(rem)) {    /* exact square root */
  1510.             freeh(rem.v);
  1511.             freeh(quo.v);
  1512.             freeh(old.v);
  1513.             *dest = try;
  1514.             return TRUE;
  1515.         }
  1516.         freeh(rem.v);
  1517.         if (i <= 0) {
  1518.             /*
  1519.             * Current try is less than or equal to the square root since
  1520.             * it is less than the quotient.  If the quotient is equal to
  1521.             * the try, we are all done.  Also, if the try is equal to the
  1522.             * old try value, we are done since no improvement occurred.
  1523.             * If not, save the improved value and loop some more.
  1524.             */
  1525.             if ((i == 0) || (zcmp(old, try) == 0)) {
  1526.                 freeh(quo.v);
  1527.                 freeh(old.v);
  1528.                 *dest = try;
  1529.                 return FALSE;
  1530.             }
  1531.             old.len = try.len;
  1532.             copyval(try, old);
  1533.         }
  1534.         /* average current try and quotent for the new try */
  1535.         zadd(try, quo, &temp);
  1536.         freeh(quo.v);
  1537.         freeh(try.v);
  1538.         try = temp;
  1539.         shiftr(try, 1L);
  1540.         quicktrim(try);
  1541.     }
  1542. }
  1543.  
  1544.  
  1545. /*
  1546.  * Take an arbitrary root of a number (to the greatest integer).
  1547.  * This uses the following iteration to get the Kth root of N:
  1548.  *    x = ((K-1) * x + N / x^(K-1)) / K
  1549.  */
  1550. void
  1551. zroot(z1, z2, dest)
  1552.     ZVALUE z1, z2, *dest;
  1553. {
  1554.     ZVALUE try, quo, old, temp, temp2;
  1555.     ZVALUE k1;            /* holds k - 1 */
  1556.     int sign;
  1557.     long i, k, highbit;
  1558.     SIUNION sival;
  1559.  
  1560.     sign = z1.sign;
  1561.     if (sign && iseven(z2))
  1562.         error("Even root of negative number");
  1563.     if (iszero(z2) || isneg(z2))
  1564.         error("Non-positive root");
  1565.     if (iszero(z1)) {    /* root of zero */
  1566.         *dest = _zero_;
  1567.         return;
  1568.     }
  1569.     if (isunit(z2)) {    /* first root */
  1570.         zcopy(z1, dest);
  1571.         return;
  1572.     }
  1573.     if (isbig(z2)) {    /* humongous root */
  1574.         *dest = _one_;
  1575.         dest->sign = (HALF)sign;
  1576.         return;
  1577.     }
  1578.     k = (istiny(z2) ? z1tol(z2) : z2tol(z2));
  1579.     highbit = zhighbit(z1);
  1580.     if (highbit < k) {    /* too high a root */
  1581.         *dest = _one_;
  1582.         dest->sign = (HALF)sign;
  1583.         return;
  1584.     }
  1585.     sival.ivalue = k - 1;
  1586.     k1.v = &sival.silow;
  1587.     k1.len = 1 + (sival.sihigh != 0);
  1588.     k1.sign = 0;
  1589.     z1.sign = 0;
  1590.     /*
  1591.      * Allocate the numbers to use for the main loop.
  1592.      * The size and high bits of the final result are correctly set here.
  1593.      * Notice that the remainder of the test value is rubbish, but this
  1594.      * is unimportant.
  1595.      */
  1596.     highbit = (highbit + k - 1) / k;
  1597.     try.len = (highbit / BASEB) + 1;
  1598.     try.v = alloc(try.len);
  1599.     try.v[try.len-1] = ((HALF)1 << (highbit % BASEB));
  1600.     try.sign = 0;
  1601.     old.v = alloc(try.len);
  1602.     old.len = 1;
  1603.     *old.v = 0;
  1604.     old.sign = 0;
  1605.     /*
  1606.      * Main divide and average loop
  1607.      */
  1608.     for (;;) {
  1609.         zpowi(try, k1, &temp);
  1610.         zquo(z1, temp, &quo);
  1611.         freeh(temp.v);
  1612.         i = zrel(try, quo);
  1613.         if (i <= 0) {
  1614.             /*
  1615.              * Current try is less than or equal to the root since it is
  1616.              * less than the quotient. If the quotient is equal to the try,
  1617.              * we are all done.  Also, if the try is equal to the old value,
  1618.              * we are done since no improvement occurred.
  1619.              * If not, save the improved value and loop some more.
  1620.              */
  1621.             if ((i == 0) || (zcmp(old, try) == 0)) {
  1622.                 freeh(quo.v);
  1623.                 freeh(old.v);
  1624.                 try.sign = (HALF)sign;
  1625.                 quicktrim(try);
  1626.                 *dest = try;
  1627.                 return;
  1628.             }
  1629.             old.len = try.len;
  1630.             copyval(try, old);
  1631.         }
  1632.         /* average current try and quotent for the new try */
  1633.         zmul(try, k1, &temp);
  1634.         freeh(try.v);
  1635.         zadd(quo, temp, &temp2);
  1636.         freeh(temp.v);
  1637.         freeh(quo.v);
  1638.         zquo(temp2, z2, &try);
  1639.         freeh(temp2.v);
  1640.     }
  1641. }
  1642.  
  1643.  
  1644. /*
  1645.  * Test to see if a number is an exact square or not.
  1646.  */
  1647. BOOL
  1648. zissquare(z)
  1649.     ZVALUE z;        /* number to be tested */
  1650. {
  1651.     long n, i;
  1652.     ZVALUE tmp;
  1653.  
  1654.     if (z.sign)        /* negative */
  1655.         return FALSE;
  1656.     while ((z.len > 1) && (*z.v == 0)) {    /* ignore trailing zero words */
  1657.         z.len--;
  1658.         z.v++;
  1659.     }
  1660.     if (isleone(z))        /* zero or one */
  1661.         return TRUE;
  1662.     n = *z.v & 0xf;        /* check mod 16 values */
  1663.     if ((n != 0) && (n != 1) && (n != 4) && (n != 9))
  1664.         return FALSE;
  1665.     n = *z.v & 0xff;    /* check mod 256 values */
  1666.     i = 0x80;
  1667.     while (((i * i) & 0xff) != n)
  1668.         if (--i <= 0)
  1669.             return FALSE;
  1670.     n = zsqrt(z, &tmp);    /* must do full square root test now */
  1671.     freeh(tmp.v);
  1672.     return n;
  1673. }
  1674.  
  1675. /* END CODE */
  1676.