home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1994 November / macformat-018.iso / Utility Spectacular / Utilities / Calc / qfunc.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-10-10  |  23.7 KB  |  1,192 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 rational arithmetic non-primitive functions
  7.  */
  8.  
  9. #include "xmath.h"
  10.  
  11.  
  12. NUMBER *_epsilon_;    /* default precision for real functions */
  13. long _epsilonprec_;    /* bits of precision for epsilon */
  14.  
  15. #if 0
  16. static char *abortmsg = "Calculation aborted";
  17. static char *memmsg = "Not enough memory";
  18. #endif
  19.  
  20.  
  21. /*
  22.  * Set the default precision for real calculations.
  23.  * The precision must be between zero and one.
  24.  */
  25. void
  26. setepsilon(q)
  27.     NUMBER *q;        /* number to be set as the new epsilon */
  28. {
  29.     NUMBER *old;
  30.  
  31.     if (qisneg(q) || qiszero(q) || (qreli(q, 1L) >= 0))
  32.         error("Epsilon value must be between zero and one");
  33.     old = _epsilon_;
  34.     _epsilonprec_ = qprecision(q);
  35.     _epsilon_ = qlink(q);
  36.     if (old)
  37.         qfree(old);
  38. }
  39.  
  40.  
  41. /*
  42.  * Return the inverse of one number modulo another.
  43.  * That is, find x such that:
  44.  *    Ax = 1 (mod B)
  45.  * Returns zero if the numbers are not relatively prime (temporary hack).
  46.  */
  47. NUMBER *
  48. qminv(q1, q2)
  49.     NUMBER *q1, *q2;
  50. {
  51.     NUMBER *r;
  52.  
  53.     if (qisfrac(q1) || qisfrac(q2))
  54.         error("Non-integers for minv");
  55.     r = qalloc();
  56.     if (zmodinv(q1->num, q2->num, &r->num)) {
  57.         qfree(r);
  58.         return qlink(&_qzero_);
  59.     }
  60.     return r;
  61. }
  62.  
  63.  
  64. /*
  65.  * Return the modulo of one number raised to another.
  66.  * Here q1 is the number to be raised, q2 is the power to raise it to,
  67.  * and q3 is the number to take the modulo with the result.
  68.  * The second and third numbers are assumed nonnegative.
  69.  * Returned answer is non-negative.
  70.  *    q4 = qpowermod(q1, q2, q3);
  71.  */
  72. NUMBER *
  73. qpowermod(q1, q2, q3)
  74.     NUMBER *q1, *q2, *q3;
  75. {
  76.     NUMBER *r;
  77.  
  78.     if (qisfrac(q1) || qisfrac(q2) || qisfrac(q3))
  79.         error("Non-integers for powermod");
  80.     r = qalloc();
  81.     zpowermod(q1->num, q2->num, q3->num, &r->num);
  82.     return r;
  83. }
  84.  
  85.  
  86. /*
  87.  * Return the power function of one number with another.
  88.  * The power must be integral.
  89.  *    q3 = qpowi(q1, q2);
  90.  */
  91. NUMBER *
  92. qpowi(q1, q2)
  93.     NUMBER *q1, *q2;
  94. {
  95.     register NUMBER *r;
  96.     BOOL invert, sign;
  97.     ZVALUE num, den, z2;
  98.  
  99.     if (qisfrac(q2))
  100.         error("Raising number to fractional power");
  101.     num = q1->num;
  102.     den = q1->den;
  103.     z2 = q2->num;
  104.     sign = (num.sign && isodd(z2));
  105.     invert = z2.sign;
  106.     num.sign = 0;
  107.     z2.sign = 0;
  108.     /*
  109.     * Check for trivial cases first.
  110.     */
  111.     if (iszero(num)) {    /* zero raised to a power */
  112.         if (invert || iszero(z2))
  113.             error("Zero raised to non-positive power");
  114.         return qlink(&_qzero_);
  115.     }
  116.     if (isunit(num) && isunit(den)) {    /* 1 or -1 raised to a power */
  117.         r = (sign ? q1 : &_qone_);
  118.         r->links++;
  119.         return r;
  120.     }
  121.     if (iszero(z2))    /* raising to zeroth power */
  122.         return qlink(&_qone_);
  123.     if (isunit(z2)) {    /* raising to power 1 or -1 */
  124.         if (invert)
  125.             return qinv(q1);
  126.         return qlink(q1);
  127.     }
  128.     /*
  129.      * Not a trivial case.  Do the real work.
  130.      */
  131.     r = qalloc();
  132.     if (!isunit(num))
  133.         zpowi(num, z2, &r->num);
  134.     if (!isunit(den))
  135.         zpowi(den, z2, &r->den);
  136.     if (invert) {
  137.         z2 = r->num;
  138.         r->num = r->den;
  139.         r->den = z2;
  140.     }
  141.     r->num.sign = sign;
  142.     return r;
  143. }
  144.  
  145.  
  146. /*
  147.  * Given the legs of a right triangle, compute its hypothenuse within
  148.  * the specified error.  This is sqrt(a^2 + b^2).
  149.  */
  150. NUMBER *
  151. qhypot(q1, q2, epsilon)
  152.     NUMBER *q1, *q2, *epsilon;
  153. {
  154.     NUMBER *tmp1, *tmp2, *tmp3;
  155.  
  156.     if (qisneg(epsilon) || qiszero(epsilon))
  157.         error("Bad epsilon value for hypot");
  158.     if (qiszero(q1))
  159.         return qabs(q2);
  160.     if (qiszero(q2))
  161.         return qabs(q1);
  162.     tmp1 = qsquare(q1);
  163.     tmp2 = qsquare(q2);
  164.     tmp3 = qadd(tmp1, tmp2);
  165.     qfree(tmp1);
  166.     qfree(tmp2);
  167.     tmp1 = qsqrt(tmp3, epsilon);
  168.     qfree(tmp3);
  169.     return tmp1;
  170. }
  171.  
  172.  
  173. /*
  174.  * Given one leg of a right triangle with unit hypothenuse, calculate
  175.  * the other leg within the specified error.  This is sqrt(1 - a^2).
  176.  * If the wantneg flag is nonzero, then negative square root is returned.
  177.  */
  178. NUMBER *
  179. qlegtoleg(q, epsilon, wantneg)
  180.     NUMBER *q, *epsilon;
  181.     BOOL wantneg;
  182. {
  183.     NUMBER *qt, *res, qtmp;
  184.     ZVALUE num, ztmp1, ztmp2;
  185.  
  186.     if (qisneg(epsilon) || qiszero(epsilon))
  187.         error("Bad epsilon value for legtoleg");
  188.     if (qisunit(q))
  189.         return qlink(&_qzero_);
  190.     if (qiszero(q)) {
  191.         if (wantneg)
  192.             return qlink(&_qnegone_);
  193.         return qlink(&_qone_);
  194.     }
  195.     num = q->num;
  196.     num.sign = 0;
  197.     if (zrel(num, q->den) >= 0)
  198.         error("Leg too large in legtoleg");
  199.     zsquare(q->den, &ztmp1);
  200.     zsquare(num, &ztmp2);
  201.     zsub(ztmp1, ztmp2, &qtmp.num);
  202.     freeh(ztmp1.v);
  203.     freeh(ztmp2.v);
  204.     qtmp.den = _one_;
  205.     qt = qsqrt(&qtmp, epsilon);
  206.     freeh(qtmp.num.v);
  207.     qtmp.num = q->den;
  208.     res = qdiv(qt, &qtmp);
  209.     qfree(qt);
  210.     qt = qbappr(res, epsilon);
  211.     qfree(res);
  212.     if (wantneg) {
  213.         res = qneg(qt);
  214.         qfree(qt);
  215.         qt = res;
  216.     }
  217.     return qt;
  218. }
  219.  
  220.  
  221. /*
  222.  * Compute the square root of a number to within the specified error.
  223.  * If the number is an exact square, the exact result is returned.
  224.  *    q3 = qsqrt(q1, q2);
  225.  */
  226. NUMBER *
  227. qsqrt(q1, epsilon)
  228.     register NUMBER *q1, *epsilon;
  229. {
  230.     long bits, bits2;    /* number of bits of precision */
  231.     int exact;
  232.     NUMBER *r;
  233.     ZVALUE t1, t2;
  234.  
  235.     if (qisneg(q1))
  236.         error("Square root of negative number");
  237.     if (qisneg(epsilon) || qiszero(epsilon))
  238.         error("Bad epsilon value for sqrt");
  239.     if (qiszero(q1))
  240.         return qlink(&_qzero_);
  241.     if (qisunit(q1))
  242.         return qlink(&_qone_);
  243.     if (qiszero(epsilon) && qisint(q1) && istiny(q1->num) && (*q1->num.v <= 3))
  244.         return qlink(&_qone_);
  245.     bits = zhighbit(epsilon->den) - zhighbit(epsilon->num) + 1;
  246.     if (bits < 0)
  247.         bits = 0;
  248.     bits2 = zhighbit(q1->num) - zhighbit(q1->den) + 1;
  249.     if (bits2 > 0)
  250.         bits += bits2;
  251.     r = qalloc();
  252.     zshift(q1->num, bits * 2, &t2);
  253.     zmul(q1->den, t2, &t1);
  254.     freeh(t2.v);
  255.     exact = zsqrt(t1, &t2);
  256.     freeh(t1.v);
  257.     if (exact) {
  258.         zshift(q1->den, bits, &t1);
  259.         zreduce(t2, t1, &r->num, &r->den);
  260.     } else {
  261.         zquo(t2, q1->den, &t1);
  262.         freeh(t2.v);
  263.         zbitvalue(bits, &t2);
  264.         zreduce(t1, t2, &r->num, &r->den);
  265.     }
  266.     freeh(t1.v);
  267.     freeh(t2.v);
  268.     if (qiszero(r)) {
  269.         qfree(r);
  270.         r = qlink(&_qzero_);
  271.     }
  272.     return r;
  273. }
  274.  
  275.  
  276. /*
  277.  * Calculate the integral part of the square root of a number.
  278.  * Example:  qisqrt(13) = 3.
  279.  */
  280. NUMBER *
  281. qisqrt(q)
  282.     NUMBER *q;
  283. {
  284.     NUMBER *r;
  285.     ZVALUE tmp;
  286.  
  287.     if (qisneg(q))
  288.         error("Square root of negative number");
  289.     if (qiszero(q))
  290.         return qlink(&_qzero_);
  291.     if (qisint(q) && istiny(q->num) && (z1tol(q->num) < 4))
  292.         return qlink(&_qone_);
  293.     r = qalloc();
  294.     if (qisint(q)) {
  295.         (void) zsqrt(q->num, &r->num);
  296.         return r;
  297.     }
  298.     zquo(q->num, q->den, &tmp);
  299.     (void) zsqrt(tmp, &r->num);
  300.     freeh(tmp.v);
  301.     return r;
  302. }
  303.  
  304.  
  305. /*
  306.  * Return whether or not a number is an exact square.
  307.  */
  308. BOOL
  309. qissquare(q)
  310.     NUMBER *q;
  311. {
  312.     BOOL flag;
  313.  
  314.     flag = zissquare(q->num);
  315.     if (qisint(q) || !flag)
  316.         return flag;
  317.     return zissquare(q->den);
  318. }
  319.  
  320.  
  321. /*
  322.  * Compute the greatest integer of the Kth root of a number.
  323.  * Example:  qiroot(85, 3) = 4.
  324.  */
  325. NUMBER *
  326. qiroot(q1, q2)
  327.     register NUMBER *q1, *q2;
  328. {
  329.     NUMBER *r;
  330.     ZVALUE tmp;
  331.  
  332.     if (qisneg(q2) || qiszero(q2) || qisfrac(q2))
  333.         error("Taking number to bad root value");
  334.     if (qiszero(q1))
  335.         return qlink(&_qzero_);
  336.     if (qisone(q1) || qisone(q2))
  337.         return qlink(q1);
  338.     if (qistwo(q2))
  339.         return qisqrt(q1);
  340.     r = qalloc();
  341.     if (qisint(q1)) {
  342.         zroot(q1->num, q2->num, &r->num);
  343.         return r;
  344.     }
  345.     zquo(q1->num, q1->den, &tmp);
  346.     zroot(tmp, q2->num, &r->num);
  347.     freeh(tmp.v);
  348.     return r;
  349. }
  350.  
  351.  
  352. /*
  353.  * Return the greatest integer of the base 2 log of a number.
  354.  * This is the number such that  1 <= x / log2(x) < 2.
  355.  * Examples:  qilog2(8) = 3, qilog2(1.3) = 1, qilog2(1/7) = -3.
  356.  */
  357. long
  358. qilog2(q)
  359.     NUMBER *q;        /* number to take log of */
  360. {
  361.     long n;            /* power of two */
  362.     int c;            /* result of comparison */
  363.     ZVALUE tmp;        /* temporary value */
  364.  
  365.     if (qisneg(q) || qiszero(q))
  366.         error("Non-positive number for log2");
  367.     if (qisint(q))
  368.         return zhighbit(q->num);
  369.     n = zhighbit(q->num) - zhighbit(q->den);
  370.     if (n == 0)
  371.         c = zrel(q->num, q->den);
  372.     else if (n > 0) {
  373.         zshift(q->den, n, &tmp);
  374.         c = zrel(q->num, tmp);
  375.     } else {
  376.         zshift(q->num, n, &tmp);
  377.         c = zrel(tmp, q->den);
  378.     }
  379.     if (n)
  380.         freeh(tmp.v);
  381.     if (c < 0)
  382.         n--;
  383.     return n;
  384. }
  385.  
  386.  
  387. /*
  388.  * Return the greatest integer of the base 10 log of a number.
  389.  * This is the number such that  1 <= x / log10(x) < 10.
  390.  * Examples:  qilog10(100) = 2, qilog10(12.3) = 1, qilog10(.023) = -2.
  391.  */
  392. long
  393. qilog10(q)
  394.     NUMBER *q;        /* number to take log of */
  395. {
  396.     long n;            /* log value */
  397.     ZVALUE temp;        /* temporary value */
  398.  
  399.     if (qisneg(q) || qiszero(q))
  400.         error("Non-positive number for log10");
  401.     if (qisint(q))
  402.         return zlog10(q->num);
  403.     /*
  404.      * The number is not an integer.
  405.      * Compute the result if the number is greater than one.
  406.      */
  407.     if ((q->num.len > q->den.len) ||
  408.         ((q->num.len == q->den.len) && (zrel(q->num, q->den) > 0))) {
  409.             zquo(q->num, q->den, &temp);
  410.             n = zlog10(temp);
  411.             freeh(temp.v);
  412.             return n;
  413.     }
  414.     /*
  415.      * Here if the number is less than one.
  416.      * If the number is the inverse of a power of ten, then the obvious answer
  417.      * will be off by one.  Subtracting one if the number is the inverse of an
  418.      * integer will fix it.
  419.      */
  420.     if (isunit(q->num))
  421.         zsub(q->den, _one_, &temp);
  422.     else
  423.         zquo(q->den, q->num, &temp);
  424.     n = -zlog10(temp) - 1;
  425.     freeh(temp.v);
  426.     return n;
  427. }
  428.  
  429.  
  430. /*
  431.  * Return the number of digits in a number, ignoring the sign.
  432.  * For fractions, this is the number of digits of its greatest integer.
  433.  * Examples: qdigits(3456) = 4, qdigits(-23.45) = 2, qdigits(.0120) = 1.
  434.  */
  435. long
  436. qdigits(q)
  437.     NUMBER *q;        /* number to count digits of */
  438. {
  439.     long n;            /* number of digits */
  440.     ZVALUE temp;        /* temporary value */
  441.  
  442.     if (qisint(q))
  443.         return zdigits(q->num);
  444.     zquo(q->num, q->den, &temp);
  445.     n = zdigits(temp);
  446.     freeh(temp.v);
  447.     return n;
  448. }
  449.  
  450.  
  451. /*
  452.  * Return the digit at the specified decimal place of a number represented
  453.  * in floating point.  The lowest digit of the integral part of a number
  454.  * is the zeroth decimal place.  Negative decimal places indicate digits
  455.  * to the right of the decimal point.  Examples: qdigit(1234.5678, 1) = 3,
  456.  * qdigit(1234.5678, -3) = 7.
  457.  */
  458. FLAG
  459. qdigit(q, n)
  460.     NUMBER *q;
  461.     long n;
  462. {
  463.     ZVALUE tenpow, tmp1, tmp2;
  464.     FLAG res;
  465.  
  466.     /*
  467.      * Zero number or negative decimal place of integer is trivial.
  468.      */
  469.     if (qiszero(q) || (qisint(q) && (n < 0)))
  470.         return 0;
  471.     /*
  472.      * For non-negative decimal places, answer is easy.
  473.      */
  474.     if (n >= 0) {
  475.         if (qisint(q))
  476.             return zdigit(q->num, n);
  477.         zquo(q->num, q->den, &tmp1);
  478.         res = zdigit(tmp1, n);
  479.         freeh(tmp1.v);
  480.         return res;
  481.     }
  482.     /*
  483.      * Fractional value and want negative digit, must work harder.
  484.      */
  485.     ztenpow(-n, &tenpow);
  486.     zmul(q->num, tenpow, &tmp1);
  487.     freeh(tenpow.v);
  488.     zquo(tmp1, q->den, &tmp2);
  489.     res = zmodi(tmp2, 10L);
  490.     freeh(tmp1.v);
  491.     freeh(tmp2.v);
  492.     return res;
  493. }
  494.  
  495.  
  496. /*
  497.  * Return whether or not a bit is set at the specified bit position in a
  498.  * number.  The lowest bit of the integral part of a number is the zeroth
  499.  * bit position.  Negative bit positions indicate bits to the right of the
  500.  * binary decimal point.  Examples: qdigit(17.1, 0) = 1, qdigit(17.1, -1) = 0.
  501.  */
  502. BOOL
  503. qisset(q, n)
  504.     NUMBER *q;
  505.     long n;
  506. {
  507.     NUMBER *qtmp1, *qtmp2;
  508.     ZVALUE ztmp;
  509.     BOOL res;
  510.  
  511.     /*
  512.      * Zero number or negative bit position place of integer is trivial.
  513.      */
  514.     if (qiszero(q) || (qisint(q) && (n < 0)))
  515.         return FALSE;
  516.     /*
  517.      * For non-negative bit positions, answer is easy.
  518.      */
  519.     if (n >= 0) {
  520.         if (qisint(q))
  521.             return zisset(q->num, n);
  522.         zquo(q->num, q->den, &ztmp);
  523.         res = zisset(ztmp, n);
  524.         freeh(ztmp.v);
  525.         return res;
  526.     }
  527.     /*
  528.      * Fractional value and want negative bit position, must work harder.
  529.      */
  530.     qtmp1 = qscale(q, -n);
  531.     qtmp2 = qint(qtmp1);
  532.     qfree(qtmp1);
  533.     res = ((qtmp2->num.v[0] & 0x01) != 0);
  534.     qfree(qtmp2);
  535.     return res;
  536. }
  537.  
  538.  
  539. /*
  540.  * Compute the factorial of an integer.
  541.  *    q2 = qfact(q1);
  542.  */
  543. NUMBER *
  544. qfact(q)
  545.     register NUMBER *q;
  546. {
  547.     register NUMBER *r;
  548.  
  549.     if (qisfrac(q))
  550.         error("Non-integral factorial");
  551.     if (qiszero(q) || isone(q->num))
  552.         return qlink(&_qone_);
  553.     r = qalloc();
  554.     zfact(q->num, &r->num);
  555.     return r;
  556. }
  557.  
  558.  
  559. /*
  560.  * Compute the product of the primes less than or equal to a number.
  561.  *    q2 = qpfact(q1);
  562.  */
  563. NUMBER *
  564. qpfact(q)
  565.     register NUMBER *q;
  566. {
  567.     NUMBER *r;
  568.  
  569.     if (qisfrac(q))
  570.         error("Non-integral factorial");
  571.     r = qalloc();
  572.     zpfact(q->num, &r->num);
  573.     return r;
  574. }
  575.  
  576.  
  577. /*
  578.  * Compute the lcm of all the numbers less than or equal to a number.
  579.  *    q2 = qlcmfact(q1);
  580.  */
  581. NUMBER *
  582. qlcmfact(q)
  583.     register NUMBER *q;
  584. {
  585.     NUMBER *r;
  586.  
  587.     if (qisfrac(q))
  588.         error("Non-integral lcmfact");
  589.     r = qalloc();
  590.     zlcmfact(q->num, &r->num);
  591.     return r;
  592. }
  593.  
  594.  
  595. /*
  596.  * Compute the permutation function  M! / (M - N)!.
  597.  */
  598. NUMBER *
  599. qperm(q1, q2)
  600.     register NUMBER *q1, *q2;
  601. {
  602.     register NUMBER *r;
  603.  
  604.     if (qisfrac(q1) || qisfrac(q2))
  605.         error("Non-integral arguments for permutation");
  606.     r = qalloc();
  607.     zperm(q1->num, q2->num, &r->num);
  608.     return r;
  609. }
  610.  
  611.  
  612. /*
  613.  * Compute the combinatorial function  M! / (N! * (M - N)!).
  614.  */
  615. NUMBER *
  616. qcomb(q1, q2)
  617.     register NUMBER *q1, *q2;
  618. {
  619.     register NUMBER *r;
  620.  
  621.     if (qisfrac(q1) || qisfrac(q2))
  622.         error("Non-integral arguments for combinatorial");
  623.     r = qalloc();
  624.     zcomb(q1->num, q2->num, &r->num);
  625.     return r;
  626. }
  627.  
  628.  
  629. /*
  630.  * Compute the Jacobi function (a / b).
  631.  * -1 => a is not quadratic residue mod b
  632.  * 1 => b is composite, or a is quad residue of b
  633.  * 0 => b is even or b < 0
  634.  */
  635. NUMBER *
  636. qjacobi(q1, q2)
  637.     register NUMBER *q1, *q2;
  638. {
  639.     if (qisfrac(q1) || qisfrac(q2))
  640.         error("Non-integral arguments for jacobi");
  641.     return itoq((long) zjacobi(q1->num, q2->num));
  642. }
  643.  
  644.  
  645. /*
  646.  * Compute the Fibonacci number F(n).
  647.  */
  648. NUMBER *
  649. qfib(q)
  650.     register NUMBER *q;
  651. {
  652.     register NUMBER *r;
  653.  
  654.     if (qisfrac(q))
  655.         error("Non-integral Fibonacci number");
  656.     r = qalloc();
  657.     zfib(q->num, &r->num);
  658.     return r;
  659. }
  660.  
  661.  
  662. /*
  663.  * Truncate a number to the specified number of decimal places.
  664.  * Specifying zero places makes the result identical to qint.
  665.  * Example: qtrunc(2/3, 3) = .666
  666.  */
  667. NUMBER *
  668. qtrunc(q1, q2)
  669.     NUMBER *q1, *q2;
  670. {
  671.     long places;
  672.     NUMBER *r;
  673.     ZVALUE tenpow, tmp1, tmp2;
  674.  
  675.     if (qisfrac(q2) || qisneg(q2) || !istiny(q2->num))
  676.         error("Bad number of places for qtrunc");
  677.     if (qisint(q1))
  678.         return qlink(q1);
  679.     r = qalloc();
  680.     places = z1tol(q2->num);
  681.     /*
  682.      * Ok, produce the result.
  683.      * First see if we want no places, in which case just take integer part.
  684.      */
  685.     if (places == 0) {
  686.         zquo(q1->num, q1->den, &r->num);
  687.         return r;
  688.     }
  689.     ztenpow(places, &tenpow);
  690.     zmul(q1->num, tenpow, &tmp1);
  691.     zquo(tmp1, q1->den, &tmp2);
  692.     freeh(tmp1.v);
  693.     if (iszero(tmp2)) {
  694.         freeh(tmp2.v);
  695.         return qlink(&_qzero_);
  696.     }
  697.     /*
  698.      * Now reduce the result to the lowest common denominator.
  699.      */
  700.     zgcd(tmp2, tenpow, &tmp1);
  701.     if (isunit(tmp1)) {
  702.         freeh(tmp1.v);
  703.         r->num = tmp2;
  704.         r->den = tenpow;
  705.         return r;
  706.     }
  707.     zquo(tmp2, tmp1, &r->num);
  708.     zquo(tenpow, tmp1, &r->den);
  709.     freeh(tmp1.v);
  710.     freeh(tmp2.v);
  711.     freeh(tenpow.v);
  712.     return r;
  713. }
  714.  
  715.  
  716. /*
  717.  * Round a number to the specified number of decimal places.
  718.  * Zero decimal places means round to the nearest integer.
  719.  * Example: qround(2/3, 3) = .667
  720.  */
  721. NUMBER *
  722. qround(q, places)
  723.     NUMBER *q;        /* number to be rounded */
  724.     long places;        /* number of decimal places to round to */
  725. {
  726.     NUMBER *r;
  727.     ZVALUE tenpow, roundval, tmp1, tmp2;
  728.  
  729.     if (places < 0)
  730.         error("Negative places for qround");
  731.     if (qisint(q))
  732.         return qlink(q);
  733.     /*
  734.      * Calculate one half of the denominator, ignoring fractional results.
  735.      * This is the value we will add in order to cause rounding.
  736.      */
  737.     zshift(q->den, -1L, &roundval);
  738.     roundval.sign = q->num.sign;
  739.     /*
  740.      * Ok, now do the actual work to produce the result.
  741.      */
  742.     r = qalloc();
  743.     ztenpow(places, &tenpow);
  744.     zmul(q->num, tenpow, &tmp2);
  745.     zadd(tmp2, roundval, &tmp1);
  746.     freeh(tmp2.v);
  747.     freeh(roundval.v);
  748.     zquo(tmp1, q->den, &tmp2);
  749.     freeh(tmp1.v);
  750.     if (iszero(tmp2)) {
  751.         freeh(tmp2.v);
  752.         return qlink(&_qzero_);
  753.     }
  754.     /*
  755.      * Now reduce the result to the lowest common denominator.
  756.      */
  757.     zgcd(tmp2, tenpow, &tmp1);
  758.     if (isunit(tmp1)) {
  759.         freeh(tmp1.v);
  760.         r->num = tmp2;
  761.         r->den = tenpow;
  762.         return r;
  763.     }
  764.     zquo(tmp2, tmp1, &r->num);
  765.     zquo(tenpow, tmp1, &r->den);
  766.     freeh(tmp1.v);
  767.     freeh(tmp2.v);
  768.     freeh(tenpow.v);
  769.     return r;
  770. }
  771.  
  772.  
  773. /*
  774.  * Truncate a number to the specified number of binary places.
  775.  * Specifying zero places makes the result identical to qint.
  776.  */
  777. NUMBER *
  778. qbtrunc(q1, q2)
  779.     NUMBER *q1, *q2;
  780. {
  781.     long places, twopow;
  782.     NUMBER *r;
  783.     ZVALUE tmp1, tmp2;
  784.  
  785.     if (qisfrac(q2) || qisneg(q2) || !istiny(q2->num))
  786.         error("Bad number of places for qtrunc");
  787.     if (qisint(q1))
  788.         return qlink(q1);
  789.     r = qalloc();
  790.     places = z1tol(q2->num);
  791.     /*
  792.      * Ok, produce the result.
  793.      * First see if we want no places, in which case just take integer part.
  794.      */
  795.     if (places == 0) {
  796.         zquo(q1->num, q1->den, &r->num);
  797.         return r;
  798.     }
  799.     zshift(q1->num, places, &tmp1);
  800.     zquo(tmp1, q1->den, &tmp2);
  801.     freeh(tmp1.v);
  802.     if (iszero(tmp2)) {
  803.         freeh(tmp2.v);
  804.         return qlink(&_qzero_);
  805.     }
  806.     /*
  807.      * Now reduce the result to the lowest common denominator.
  808.      */
  809.     if (isodd(tmp2)) {
  810.         r->num = tmp2;
  811.         zbitvalue(places, &r->den);
  812.         return r;
  813.     }
  814.     twopow = zlowbit(tmp2);
  815.     if (twopow > places)
  816.         twopow = places;
  817.     places -= twopow;
  818.     zshift(tmp2, -twopow, &r->num);
  819.     freeh(tmp2.v);
  820.     zbitvalue(places, &r->den);
  821.     return r;
  822. }
  823.  
  824.  
  825. /*
  826.  * Round a number to the specified number of binary places.
  827.  * Zero binary places means round to the nearest integer.
  828.  */
  829. NUMBER *
  830. qbround(q, places)
  831.     NUMBER *q;        /* number to be rounded */
  832.     long places;        /* number of binary places to round to */
  833. {
  834.     long twopow;
  835.     NUMBER *r;
  836.     ZVALUE roundval, tmp1, tmp2;
  837.  
  838.     if (places < 0)
  839.         error("Negative places for qbround");
  840.     if (qisint(q))
  841.         return qlink(q);
  842.     r = qalloc();
  843.     /*
  844.      * Calculate one half of the denominator, ignoring fractional results.
  845.      * This is the value we will add in order to cause rounding.
  846.      */
  847.     zshift(q->den, -1L, &roundval);
  848.     roundval.sign = q->num.sign;
  849.     /*
  850.      * Ok, produce the result.
  851.      */
  852.     zshift(q->num, places, &tmp1);
  853.     zadd(tmp1, roundval, &tmp2);
  854.     freeh(roundval.v);
  855.     freeh(tmp1.v);
  856.     zquo(tmp2, q->den, &tmp1);
  857.     freeh(tmp2.v);
  858.     if (iszero(tmp1)) {
  859.         freeh(tmp1.v);
  860.         return qlink(&_qzero_);
  861.     }
  862.     /*
  863.      * Now reduce the result to the lowest common denominator.
  864.      */
  865.     if (isodd(tmp1)) {
  866.         r->num = tmp1;
  867.         zbitvalue(places, &r->den);
  868.         return r;
  869.     }
  870.     twopow = zlowbit(tmp1);
  871.     if (twopow > places)
  872.         twopow = places;
  873.     places -= twopow;
  874.     zshift(tmp1, -twopow, &r->num);
  875.     freeh(tmp1.v);
  876.     zbitvalue(places, &r->den);
  877.     return r;
  878. }
  879.  
  880.  
  881. /*
  882.  * Approximate a number by using binary rounding with the minimum number
  883.  * of binary places so that the resulting number is within the specified
  884.  * epsilon of the original number.
  885.  */
  886. NUMBER *
  887. qbappr(q, e)
  888.     NUMBER *q, *e;
  889. {
  890.     long bits;
  891.  
  892.     if (qisneg(e) || qiszero(e))
  893.         error("Bad epsilon value for qbappr");
  894.     if (e == _epsilon_)
  895.         bits = _epsilonprec_ + 1;
  896.     else
  897.         bits = qprecision(e) + 1;
  898.     return qbround(q, bits);
  899. }
  900.  
  901.  
  902. /*
  903.  * Approximate a number using continued fractions until the approximation
  904.  * error is less than the specified value.  If a NULL pointer is given
  905.  * for the error value, then the closest simpler fraction is returned.
  906.  */
  907. NUMBER *
  908. qcfappr(q, e)
  909.     NUMBER *q, *e;
  910. {
  911.     NUMBER qtest, *qtmp;
  912.     ZVALUE u1, u2, u3, v1, v2, v3, t1, t2, t3, qq, tt;
  913.     int i;
  914.     BOOL haveeps;
  915.  
  916.     haveeps = TRUE;
  917.     if (e == NULL) {
  918.         haveeps = FALSE;
  919.         e = &_qzero_;
  920.     }
  921.     if (qisneg(e))
  922.         error("Negative epsilon for cfappr");
  923.     if (qisint(q) || isunit(q->num) || (haveeps && qiszero(e)))
  924.         return qlink(q);
  925.     u1 = _one_;
  926.     u2 = _zero_;
  927.     u3 = q->num;
  928.     u3.sign = 0;
  929.     v1 = _zero_;
  930.     v2 = _one_;
  931.     v3 = q->den;
  932.     while (!iszero(v3)) {
  933.         if (!iszero(u2) && !iszero(u1)) {
  934.             qtest.num = u2;
  935.             qtest.den = u1;
  936.             qtest.den.sign = 0;
  937.             qtest.num.sign = q->num.sign;
  938.             qtmp = qsub(q, &qtest);
  939.             qtest = *qtmp;
  940.             qtest.num.sign = 0;
  941.             i = qrel(&qtest, e);
  942.             qfree(qtmp);
  943.             if (i <= 0)
  944.                 break;
  945.         }
  946.         zquo(u3, v3, &qq);
  947.         zmul(qq, v1, &tt); zsub(u1, tt, &t1); freeh(tt.v);
  948.         zmul(qq, v2, &tt); zsub(u2, tt, &t2); freeh(tt.v);
  949.         zmul(qq, v3, &tt); zsub(u3, tt, &t3); freeh(tt.v);
  950.         freeh(qq.v); freeh(u1.v); freeh(u2.v);
  951.         if ((u3.v != q->num.v) && (u3.v != q->den.v))
  952.             freeh(u3.v);
  953.         u1 = v1; u2 = v2; u3 = v3;
  954.         v1 = t1; v2 = t2; v3 = t3;
  955.     }
  956.     if (u3.v != q->den.v)
  957.         freeh(u3.v);
  958.     freeh(v1.v);
  959.     freeh(v2.v);
  960.     i = iszero(v3);
  961.     freeh(v3.v);
  962.     if (i && haveeps) {
  963.         freeh(u1.v);
  964.         freeh(u2.v);
  965.         return qlink(q);
  966.     }
  967.     qtest.num = u2;
  968.     qtest.den = u1;
  969.     qtest.den.sign = 0;
  970.     qtest.num.sign = q->num.sign;
  971.     qtmp = qcopy(&qtest);
  972.     freeh(u1.v);
  973.     freeh(u2.v);
  974.     return qtmp;
  975. }
  976.  
  977.  
  978. /*
  979.  * Return an indication on whether or not two fractions are approximately
  980.  * equal within the specified epsilon. Returns negative if the absolute value
  981.  * of the difference between the two numbers is less than epsilon, zero if
  982.  * the difference is equal to epsilon, and positive if the difference is
  983.  * greater than epsilon.
  984.  */
  985. FLAG
  986. qnear(q1, q2, epsilon)
  987.     NUMBER *q1, *q2, *epsilon;
  988. {
  989.     int res;
  990.     NUMBER qtemp, *qq;
  991.  
  992.     if (qisneg(epsilon))
  993.         error("Negative epsilon for near");
  994.     if (q1 == q2) {
  995.         if (qiszero(epsilon))
  996.             return 0;
  997.         return -1;
  998.     }
  999.     if (qiszero(epsilon))
  1000.         return qcmp(q1, q2);
  1001.     if (qiszero(q2)) {
  1002.         qtemp = *q1;
  1003.         qtemp.num.sign = 0;
  1004.         return qrel(&qtemp, epsilon);
  1005.     }
  1006.     if (qiszero(q1)) {
  1007.         qtemp = *q2;
  1008.         qtemp.num.sign = 0;
  1009.         return qrel(&qtemp, epsilon);
  1010.     }
  1011.     qq = qsub(q1, q2);
  1012.     qtemp = *qq;
  1013.     qtemp.num.sign = 0;
  1014.     res = qrel(&qtemp, epsilon);
  1015.     qfree(qq);
  1016.     return res;
  1017. }
  1018.  
  1019.  
  1020. /*
  1021.  * Compute the gcd (greatest common divisor) of two numbers.
  1022.  *    q3 = qgcd(q1, q2);
  1023.  */
  1024. NUMBER *
  1025. qgcd(q1, q2)
  1026.     register NUMBER *q1, *q2;
  1027. {
  1028.     ZVALUE z;
  1029.     NUMBER *q;
  1030.  
  1031.     if (qisfrac(q1) || qisfrac(q2))
  1032.         error("Non-integers for gcd");
  1033.     zgcd(q1->num, q2->num, &z);
  1034.     if (isunit(z)) {
  1035.         freeh(z.v);
  1036.         return qlink(&_qone_);
  1037.     }
  1038.     q = qalloc();
  1039.     q->num = z;
  1040.     return q;
  1041. }
  1042.  
  1043.  
  1044. /*
  1045.  * Compute the lcm (least common denominator) of two numbers.
  1046.  *    q3 = qlcm(q1, q2);
  1047.  */
  1048. NUMBER *
  1049. qlcm(q1, q2)
  1050.     register NUMBER *q1, *q2;
  1051. {
  1052.     NUMBER *q;
  1053.  
  1054.     if (qisfrac(q1) || qisfrac(q2))
  1055.         error("Non-integers for lcm");
  1056.     if (qisunit(q1))
  1057.         return qlink(q2);
  1058.     if (qisunit(q2))
  1059.         return qlink(q1);
  1060.     q = qalloc();
  1061.     zlcm(q1->num, q2->num, &q->num);
  1062.     return q;
  1063. }
  1064.  
  1065.  
  1066. /*
  1067.  * Remove all occurances of the specified factor from a number.
  1068.  * Returned number is always positive.
  1069.  */
  1070. NUMBER *
  1071. qfacrem(q1, q2)
  1072.     NUMBER *q1, *q2;
  1073. {
  1074.     long count;
  1075.     ZVALUE tmp;
  1076.     NUMBER *r;
  1077.  
  1078.     if (qisfrac(q1) || qisfrac(q2))
  1079.         error("Non-integers for factor removal");
  1080.     count = zfacrem(q1->num, q2->num, &tmp);
  1081.     if (isunit(tmp)) {
  1082.         freeh(tmp.v);
  1083.         return qlink(&_qone_);
  1084.     }
  1085.     if (count == 0) {
  1086.         freeh(tmp.v);
  1087.         return qlink(q1);
  1088.     }
  1089.     r = qalloc();
  1090.     r->num = tmp;
  1091.     return r;
  1092. }
  1093.  
  1094.  
  1095. /*
  1096.  * Divide one number by the gcd of it with another number repeatedly until
  1097.  * the number is relatively prime.
  1098.  */
  1099. NUMBER *
  1100. qgcdrem(q1, q2)
  1101.     NUMBER *q1, *q2;
  1102. {
  1103.     ZVALUE tmp;
  1104.     NUMBER *r;
  1105.  
  1106.     if (qisfrac(q1) || qisfrac(q2))
  1107.         error("Non-integers for gcdrem");
  1108.     zgcdrem(q1->num, q2->num, &tmp);
  1109.     if (isunit(tmp)) {
  1110.         freeh(tmp.v);
  1111.         return qlink(&_qone_);
  1112.     }
  1113.     if (zcmp(q1->num, tmp) == 0) {
  1114.         freeh(tmp.v);
  1115.         return qlink(q1);
  1116.     }
  1117.     r = qalloc();
  1118.     r->num = tmp;
  1119.     return r;
  1120. }
  1121.  
  1122.  
  1123. /*
  1124.  * Return the lowest prime factor of a number.
  1125.  * Search is conducted for the specified number of primes.
  1126.  * Returns one if no factor was found.
  1127.  */
  1128. NUMBER *
  1129. qlowfactor(q1, q2)
  1130.     NUMBER *q1, *q2;
  1131. {
  1132.     if (qisfrac(q1) || qisfrac(q2))
  1133.         error("Non-integers for lowfactor");
  1134.     return itoq(zlowfactor(q1->num, ztoi(q2->num)));
  1135. }
  1136.  
  1137.  
  1138. /*
  1139.  * Return the number of places after the decimal point needed to exactly
  1140.  * represent the specified number as a real number.  Integers return zero,
  1141.  * and non-terminating decimals return minus one.  Examples:
  1142.  *    qplaces(1/7)=-1, qplaces(3/10)= 1, qplaces(1/8)=3, qplaces(4)=0.
  1143.  */
  1144. long
  1145. qplaces(q)
  1146.     NUMBER *q;
  1147. {
  1148.     long twopow, fivepow;
  1149.     HALF fiveval[2];
  1150.     ZVALUE five;
  1151.     ZVALUE tmp;
  1152.  
  1153.     if (qisint(q))    /* no decimal places if number is integer */
  1154.         return 0;
  1155.     /*
  1156.      * The number of decimal places of a fraction in lowest terms is finite
  1157.      * if an only if the denominator is of the form 2^A * 5^B, and then the
  1158.      * number of decimal places is equal to MAX(A, B).
  1159.      */
  1160.     five.sign = 0;
  1161.     five.len = 1;
  1162.     five.v = fiveval;
  1163.     fiveval[0] = 5;
  1164.     fivepow = zfacrem(q->den, five, &tmp);
  1165.     if (!zisonebit(tmp)) {
  1166.         freeh(tmp.v);
  1167.         return -1;
  1168.     }
  1169.     twopow = zlowbit(tmp);
  1170.     freeh(tmp.v);
  1171.     if (twopow < fivepow)
  1172.         twopow = fivepow;
  1173.     return twopow;
  1174. }
  1175.  
  1176.  
  1177. /*
  1178.  * Perform a probabilistic primality test (algorithm P in Knuth).
  1179.  * Returns FALSE if definitely not prime, or TRUE if probably prime.
  1180.  * Second arg determines how many times to check for primality.
  1181.  */
  1182. BOOL
  1183. qprimetest(q1, q2)
  1184.     NUMBER *q1, *q2;
  1185. {
  1186.     if (qisfrac(q1) || qisfrac(q2) || qisneg(q2))
  1187.         error("Bad arguments for qprimetest");
  1188.     return zprimetest(q1->num, qtoi(q2));
  1189. }
  1190.  
  1191. /* END CODE */
  1192.