home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1994 November / macformat-018.iso / Utility Spectacular / Utilities / Calc / qtrans.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-10-10  |  20.1 KB  |  938 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.  * Transcendental functions for real numbers.
  7.  * These are sin, cos, exp, ln, power, cosh, sinh.
  8.  */
  9.  
  10. #include "xmath.h"
  11.  
  12. BOOL _sinisneg_;    /* whether sin(x) < 0 (set by cos(x)) */
  13.  
  14.  
  15. /*
  16.  * Calculate the cosine of a number with an accuracy within epsilon.
  17.  * This also saves the sign of the corresponding sin function.
  18.  */
  19. NUMBER *
  20. qcos(q, epsilon)
  21.     NUMBER *q, *epsilon;
  22. {
  23.     NUMBER *term, *sum, *qsq, *epsilon2, *tmp;
  24.     FULL n, i;
  25.     long scale, bits, bits2;
  26.  
  27.     _sinisneg_ = qisneg(q);
  28.     if (qisneg(epsilon) || qiszero(epsilon))
  29.         error("Illegal epsilon value for cosine");
  30.     if (qiszero(q))
  31.         return qlink(&_qone_);
  32.     bits = qprecision(epsilon) + 1;
  33.     epsilon = qscale(epsilon, -4L);
  34.     /*
  35.      * If the argument is larger than one, then divide it by a power of two
  36.      * so that it is one or less.  This will make the series converge quickly.
  37.      * We will extrapolate the result for the original argument afterwards.
  38.      */
  39.     scale = zhighbit(q->num) - zhighbit(q->den) + 1;
  40.     if (scale < 0)
  41.         scale = 0;
  42.     if (scale > 0) {
  43.         q = qscale(q, -scale);
  44.         tmp = qscale(epsilon, -scale);
  45.         qfree(epsilon);
  46.         epsilon = tmp;
  47.     }
  48.     epsilon2 = qscale(epsilon, -4L);
  49.     qfree(epsilon);
  50.     bits2 = qprecision(epsilon2) + 10;
  51.     /*
  52.      * Now use the Taylor series expansion to calculate the cosine.
  53.      * Keep using approximations so that the fractions don't get too large.
  54.      */
  55.     qsq = qsquare(q);
  56.     if (scale > 0)
  57.         qfree(q);
  58.     term = qlink(&_qone_);
  59.     sum = qlink(&_qone_);
  60.     n = 0;
  61.     while (qrel(term, epsilon2) > 0) {
  62.         i = ++n;
  63.         i *= ++n;
  64.         tmp = qmul(term, qsq);
  65.         qfree(term);
  66.         term = qdivi(tmp, (long) i);
  67.         qfree(tmp);
  68.         tmp = qbround(term, bits2);
  69.         qfree(term);
  70.         term = tmp;
  71.         if (n & 2)
  72.             tmp = qsub(sum, term);
  73.         else
  74.             tmp = qadd(sum, term);
  75.         qfree(sum);
  76.         sum = qbround(tmp, bits2);
  77.         qfree(tmp);
  78.     }
  79.     qfree(term);
  80.     qfree(qsq);
  81.     qfree(epsilon2);
  82.     /*
  83.      * Now scale back up to the original value of x by using the formula:
  84.      *    cos(2 * x) = 2 * (cos(x) ^ 2) - 1.
  85.      */
  86.     while (--scale >= 0) {
  87.         if (qisneg(sum))
  88.             _sinisneg_ = !_sinisneg_;
  89.         tmp = qsquare(sum);
  90.         qfree(sum);
  91.         sum = qscale(tmp, 1L);
  92.         qfree(tmp);
  93.         tmp = qdec(sum);
  94.         qfree(sum);
  95.         sum = qbround(tmp, bits2);
  96.         qfree(tmp);
  97.     }
  98.     tmp = qbround(sum, bits);
  99.     qfree(sum);
  100.     return tmp;
  101. }
  102.  
  103.  
  104. /*
  105.  * Calculate the sine of a number with an accuracy within epsilon.
  106.  * This is calculated using the formula:
  107.  *    sin(x)^2 + cos(x)^2 = 1.
  108.  * The only tricky bit is resolving the sign of the result.
  109.  * Future: Use sin(3*x) = 3*sin(x) - 4*sin(x)^3.
  110.  */
  111. NUMBER *
  112. qsin(q, epsilon)
  113.     NUMBER *q, *epsilon;
  114. {
  115.     NUMBER *tmp1, *tmp2, *epsilon2;
  116.  
  117.     if (qisneg(epsilon) || qiszero(epsilon))
  118.         error("Illegal epsilon value for sine");
  119.     if (qiszero(q))
  120.         return qlink(q);
  121.     epsilon2 = qscale(epsilon, -4L);
  122.     tmp1 = qcos(q, epsilon2);
  123.     qfree(epsilon2);
  124.     tmp2 = qlegtoleg(tmp1, epsilon, _sinisneg_);
  125.     qfree(tmp1);
  126.     return tmp2;
  127. }
  128.  
  129.  
  130. /*
  131.  * Calculate the tangent function.
  132.  */
  133. NUMBER *
  134. qtan(q, epsilon)
  135.     NUMBER *q, *epsilon;
  136. {
  137.     NUMBER *cosval, *sinval, *epsilon2, *tmp, *res;
  138.  
  139.     if (qisneg(epsilon) || qiszero(epsilon))
  140.         error("Illegal epsilon value for tangent");
  141.     if (qiszero(q))
  142.         return qlink(q);
  143.     epsilon2 = qscale(epsilon, -4L);
  144.     cosval = qcos(q, epsilon2);
  145.     sinval = qlegtoleg(cosval, epsilon2, _sinisneg_);
  146.     qfree(epsilon2);
  147.     tmp = qdiv(sinval, cosval);
  148.     qfree(cosval);
  149.     qfree(sinval);
  150.     res = qbround(tmp, qprecision(epsilon) + 1);
  151.     qfree(tmp);
  152.     return res;
  153. }
  154.  
  155.  
  156. /*
  157.  * Calculate the arcsine function.
  158.  * The result is in the range -pi/2 to pi/2.
  159.  */
  160. NUMBER *
  161. qasin(q, epsilon)
  162.     NUMBER *q, *epsilon;
  163. {
  164.     NUMBER *sum, *term, *epsilon2, *qsq, *tmp;
  165.     FULL n, i;
  166.     long bits, bits2;
  167.     int neg;
  168.     NUMBER mulnum;
  169.     HALF numval[2];
  170.     HALF denval[2];
  171.  
  172.     if (qisneg(epsilon) || qiszero(epsilon))
  173.         error("Illegal epsilon value for arcsine");
  174.     if (qiszero(q))
  175.         return qlink(&_qzero_);
  176.     if ((qrel(q, &_qone_) > 0) || (qrel(q, &_qnegone_) < 0))
  177.         error("Argument too large for asin");
  178.     neg = qisneg(q);
  179.     q = qabs(q);
  180.     epsilon = qscale(epsilon, -4L);
  181.     epsilon2 = qscale(epsilon, -4L);
  182.     mulnum.num.sign = 0;
  183.     mulnum.num.len = 1;
  184.     mulnum.num.v = numval;
  185.     mulnum.den.sign = 0;
  186.     mulnum.den.len = 1;
  187.     mulnum.den.v = denval;
  188.     /*
  189.      * If the argument is too near one (we use .5) then reduce the
  190.      * argument to a more accurate range using the formula:
  191.      *    asin(x) = 2 * asin(sqrt((1 - sqrt(1 - x^2)) / 2)).
  192.      */
  193.     if (qrel(q, &_qonehalf_) > 0) {
  194.         sum = qlegtoleg(q, epsilon2, FALSE);
  195.         qfree(q);
  196.         tmp = qsub(&_qone_, sum);
  197.         qfree(sum);
  198.         sum = qscale(tmp, -1L);
  199.         qfree(tmp);
  200.         tmp = qsqrt(sum, epsilon2);
  201.         qfree(sum);
  202.         qfree(epsilon2);
  203.         sum = qasin(tmp, epsilon);
  204.         qfree(tmp);
  205.         qfree(epsilon);
  206.         tmp = qscale(sum, 1L);
  207.         qfree(sum);
  208.         if (neg) {
  209.             sum = qneg(tmp);
  210.             qfree(tmp);
  211.             tmp = sum;
  212.         }
  213.         return tmp;
  214.     }
  215.     /*
  216.      * Argument is between zero and .5, so use the series.
  217.      */
  218.     epsilon = qscale(epsilon, -4L);
  219.     epsilon2 = qscale(epsilon, -4L);
  220.     bits = qprecision(epsilon) + 1;
  221.     bits2 = bits + 10;
  222.     sum = qlink(q);
  223.     term = qlink(q);
  224.     qsq = qsquare(q);
  225.     qfree(q);
  226.     n = 1;
  227.     while (qrel(term, epsilon2) > 0) {
  228.         i = n * n;
  229.         numval[0] = i & BASE1;
  230.         if (i >= BASE) {
  231.             numval[1] = i / BASE;
  232.             mulnum.den.len = 2;
  233.         }
  234.         i = (n + 1) * (n + 2);
  235.         denval[0] = i & BASE1;
  236.         if (i >= BASE) {
  237.             denval[1] = i / BASE;
  238.             mulnum.den.len = 2;
  239.         }
  240.         tmp = qmul(term, qsq);
  241.         qfree(term);
  242.         term = qmul(tmp, &mulnum);
  243.         qfree(tmp);
  244.         tmp = qbround(term, bits2);
  245.         qfree(term);
  246.         term = tmp;
  247.         tmp = qadd(sum, term);
  248.         qfree(sum);
  249.         sum = qbround(tmp, bits2);
  250.         qfree(tmp);
  251.         n += 2;
  252.     }
  253.     qfree(epsilon);
  254.     qfree(epsilon2);
  255.     qfree(term);
  256.     qfree(qsq);
  257.     tmp = qbround(sum, bits);
  258.     qfree(sum);
  259.     if (neg) {
  260.         term = qneg(tmp);
  261.         qfree(tmp);
  262.         tmp = term;
  263.     }
  264.     return tmp;
  265. }
  266.  
  267.  
  268. /*
  269.  * Calculate the acos function.
  270.  * The result is in the range 0 to pi.
  271.  */
  272. NUMBER *
  273. qacos(q, epsilon)
  274.     NUMBER *q, *epsilon;
  275. {
  276.     NUMBER *tmp1, *tmp2, *epsilon2;
  277.  
  278.     if (qisneg(epsilon) || qiszero(epsilon))
  279.         error("Illegal epsilon value for arccosine");
  280.     if (qisone(q))
  281.         return qlink(&_qzero_);
  282.     if ((qrel(q, &_qone_) > 0) || (qrel(q, &_qnegone_) < 0))
  283.         error("Argument too large for acos");
  284.     /*
  285.      * Calculate the result using the formula:
  286.      *    acos(x) = asin(sqrt(1 - x^2)).
  287.      */
  288.     epsilon2 = qscale(epsilon, -8L);
  289.     tmp1 = qlegtoleg(q, epsilon2, FALSE);
  290.     qfree(epsilon2);
  291.     tmp2 = qasin(tmp1, epsilon);
  292.     qfree(tmp1);
  293.     return tmp2;
  294. }
  295.  
  296.  
  297. /*
  298.  * Calculate the arctangent function with a accuracy less than epsilon.
  299.  * This uses the formula:
  300.  *    atan(x) = asin(sqrt(x^2 / (x^2 + 1))).
  301.  */
  302. NUMBER *
  303. qatan(q, epsilon)
  304.     NUMBER *q, *epsilon;
  305. {
  306.     NUMBER *tmp1, *tmp2, *tmp3, *epsilon2;
  307.  
  308.     if (qisneg(epsilon) || qiszero(epsilon))
  309.         error("Illegal epsilon value for arctangent");
  310.     if (qiszero(q))
  311.         return qlink(&_qzero_);
  312.     tmp1 = qsquare(q);
  313.     tmp2 = qinc(tmp1);
  314.     tmp3 = qdiv(tmp1, tmp2);
  315.     qfree(tmp1);
  316.     qfree(tmp2);
  317.     epsilon2 = qscale(epsilon, -8L);
  318.     tmp1 = qsqrt(tmp3, epsilon2);
  319.     qfree(epsilon2);
  320.     qfree(tmp3);
  321.     tmp2 = qasin(tmp1, epsilon);
  322.     qfree(tmp1);
  323.     if (qisneg(q)) {
  324.         tmp1 = qneg(tmp2);
  325.         qfree(tmp2);
  326.         tmp2 = tmp1;
  327.     }
  328.     return tmp2;
  329. }
  330.  
  331.  
  332. /*
  333.  * Calculate the angle which is determined by the point (x,y).
  334.  * This is the same as arctan for non-negative x, but gives the correct
  335.  * value for negative x.  By convention, y is the first argument.
  336.  * For example, qatan2(1, -1) = 3/4 * pi.
  337.  */
  338. NUMBER *
  339. qatan2(qy, qx, epsilon)
  340.     NUMBER *qy, *qx, *epsilon;
  341. {
  342.     NUMBER *tmp1, *tmp2, *epsilon2;
  343.  
  344.     if (qisneg(epsilon) || qiszero(epsilon))
  345.         error("Illegal epsilon value for atan2");
  346.     if (qiszero(qy) && qiszero(qx))
  347.         error("Zero coordinates for atan2");
  348.     /*
  349.      * If the point is on the negative real axis, then the answer is pi.
  350.      */
  351.     if (qiszero(qy) && qisneg(qx))
  352.         return qpi(epsilon);
  353.     /*
  354.      * If the point is in the right half plane, then use the normal atan.
  355.      */
  356.     if (!qisneg(qx) && !qiszero(qx)) {
  357.         if (qiszero(qy))
  358.             return qlink(&_qzero_);
  359.         tmp1 = qdiv(qy, qx);
  360.         tmp2 = qatan(tmp1, epsilon);
  361.         qfree(tmp1);
  362.         return tmp2;
  363.     }
  364.     /*
  365.      * The point is in the left half plane.  Calculate the angle by finding
  366.      * the atan of half the angle using the formula:
  367.      *    atan2(y,x) = 2 * atan((sqrt(x^2 + y^2) - x) / y).
  368.      */
  369.     epsilon2 = qscale(epsilon, -4L);
  370.     tmp1 = qhypot(qx, qy, epsilon2);
  371.     tmp2 = qsub(tmp1, qx);
  372.     qfree(tmp1);
  373.     tmp1 = qdiv(tmp2, qy);
  374.     qfree(tmp2);
  375.     tmp2 = qatan(tmp1, epsilon2);
  376.     qfree(tmp1);
  377.     qfree(epsilon2);
  378.     tmp1 = qscale(tmp2, 1L);
  379.     qfree(tmp2);
  380.     return tmp1;
  381. }
  382.  
  383.  
  384. /*
  385.  * Calculate the value of pi to within the required epsilon.
  386.  * This uses the following formula which only needs integer calculations
  387.  * except for the final operation:
  388.  *    pi = 1 / SUMOF(comb(2 * N, N) ^ 3 * (42 * N + 5) / 2 ^ (12 * N + 4)),
  389.  * where the summation runs from N=0.  This formula gives about 6 bits of
  390.  * accuracy per term.  Since the denominator for each term is a power of two,
  391.  * we can simply use shifts to sum the terms.  The combinatorial numbers
  392.  * in the formula are calculated recursively using the formula:
  393.  *    comb(2*(N+1), N+1) = 2 * comb(2 * N, N) * (2 * N + 1) / N.
  394.  */
  395. NUMBER *
  396. qpi(epsilon)
  397.     NUMBER *epsilon;
  398. {
  399.     ZVALUE comb;            /* current combinatorial value */
  400.     ZVALUE sum;            /* current sum */
  401.     ZVALUE tmp1, tmp2;
  402.     NUMBER *r, *t1, qtmp;
  403.     long shift;            /* current shift of result */
  404.     long N;                /* current term number */
  405.     long bits;            /* needed number of bits of precision */
  406.     long t;
  407.  
  408.     if (qiszero(epsilon) || qisneg(epsilon))
  409.         error("Bad epsilon value for pi");
  410.     bits = qprecision(epsilon) + 4;
  411.     comb = _one_;
  412.     itoz(5L, &sum);
  413.     N = 0;
  414.     shift = 4;
  415.     do {
  416.         t = 1 + (++N & 0x1);
  417.         (void) zdivi(comb, N / (3 - t), &tmp1);
  418.         freeh(comb.v);
  419.         zmuli(tmp1, t * (2 * N - 1), &comb);
  420.         freeh(tmp1.v);
  421.         zsquare(comb, &tmp1);
  422.         zmul(comb, tmp1, &tmp2);
  423.         freeh(tmp1.v);
  424.         zmuli(tmp2, 42 * N + 5, &tmp1);
  425.         freeh(tmp2.v);
  426.         zshift(sum, 12L, &tmp2);
  427.         freeh(sum.v);
  428.         zadd(tmp1, tmp2, &sum);
  429.         t = zhighbit(tmp1);
  430.         freeh(tmp1.v);
  431.         freeh(tmp2.v);
  432.         shift += 12;
  433.     } while ((shift - t) < bits);
  434.     qtmp.num = _one_;
  435.     qtmp.den = sum;
  436.     t1 = qscale(&qtmp, shift);
  437.     freeh(sum.v);
  438.     r = qbround(t1, bits);
  439.     qfree(t1);
  440.     return r;
  441. }
  442.  
  443.  
  444. /*
  445.  * Calculate the exponential function with a relative accuracy less than
  446.  * epsilon.
  447.  */
  448. NUMBER *
  449. qexp(q, epsilon)
  450.     NUMBER *q, *epsilon;
  451. {
  452.     long scale;
  453.     FULL n;
  454.     long bits, bits2;
  455.     NUMBER *sum, *term, *qs, *epsilon2, *tmp;
  456.  
  457.     if (qisneg(epsilon) || qiszero(epsilon))
  458.         error("Illegal epsilon value for exp");
  459.     if (qiszero(q))
  460.         return qlink(&_qone_);
  461.     epsilon = qscale(epsilon, -4L);
  462.     /*
  463.      * If the argument is larger than one, then divide it by a power of two
  464.      * so that it is one or less.  This will make the series converge quickly.
  465.      * We will extrapolate the result for the original argument afterwards.
  466.      * Also make the argument non-negative.
  467.      */
  468.     qs = qabs(q);
  469.     scale = zhighbit(q->num) - zhighbit(q->den) + 1;
  470.     if (scale < 0)
  471.         scale = 0;
  472.     if (scale > 0) {
  473.         if (scale > 100000)
  474.             error("Very large argument for exp");
  475.         tmp = qscale(qs, -scale);
  476.         qfree(qs);
  477.         qs = tmp;
  478.         tmp = qscale(epsilon, -scale);
  479.         qfree(epsilon);
  480.         epsilon = tmp;
  481.     }
  482.     epsilon2 = qscale(epsilon, -4L);
  483.     bits = qprecision(epsilon) + 1;
  484.     bits2 = bits + 10;
  485.     qfree(epsilon);
  486.     /*
  487.      * Now use the Taylor series expansion to calculate the exponential.
  488.      * Keep using approximations so that the fractions don't get too large.
  489.      */
  490.     sum = qlink(&_qone_);
  491.     term = qlink(&_qone_);
  492.     n = 0;
  493.     while (qrel(term, epsilon2) > 0) {
  494.         n++;
  495.         tmp = qmul(term, qs);
  496.         qfree(term);
  497.         term = qdivi(tmp, (long) n);
  498.         qfree(tmp);
  499.         tmp = qbround(term, bits2);
  500.         qfree(term);
  501.         term = tmp;
  502.         tmp = qadd(sum, term);
  503.         qfree(sum);
  504.         sum = qbround(tmp, bits2);
  505.         qfree(tmp);
  506.     }
  507.     qfree(term);
  508.     qfree(qs);
  509.     qfree(epsilon2);
  510.     /*
  511.      * Now repeatedly square the answer to get the final result.
  512.      * Then invert it if the original argument was negative.
  513.      */
  514.     while (--scale >= 0) {
  515.         tmp = qsquare(sum);
  516.         qfree(sum);
  517.         sum = qbround(tmp, bits2);
  518.         qfree(tmp);
  519.     }
  520.     tmp = qbround(sum, bits);
  521.     qfree(sum);
  522.     if (qisneg(q)) {
  523.         sum = qinv(tmp);
  524.         qfree(tmp);
  525.         tmp = sum;
  526.     }
  527.     return tmp;
  528. }
  529.  
  530.  
  531. /*
  532.  * Calculate the natural logarithm of a number accurate to the specified
  533.  * epsilon.
  534.  */
  535. NUMBER *
  536. qln(q, epsilon)
  537.     NUMBER *q, *epsilon;
  538. {
  539.     NUMBER *term, *term2, *sum, *epsilon2, *tmp1, *tmp2;
  540.     NUMBER *minr, *maxr;
  541.     long shift, bits, bits2;
  542.     int neg;
  543.     FULL n;
  544.  
  545.     if (qisneg(q) || qiszero(q))
  546.         error("log of non-positive number");
  547.     if (qisneg(epsilon) || qiszero(epsilon))
  548.         error("Illegal epsilon for ln");
  549.     epsilon = qscale(epsilon, -4L);
  550.     epsilon2 = qscale(epsilon, -4L);
  551.     bits = qprecision(epsilon) + 1;
  552.     bits2 = bits + 10;
  553.     qfree(epsilon);
  554.     /*
  555.      * Scale down the logarithm to the convergent range by taking square
  556.      * roots.  The result will be modified to get the final value later.
  557.      */
  558.     q = qlink(q);
  559.     minr = iitoq(7L, 10L);
  560.     maxr = iitoq(7L, 5L);
  561.     shift = 1;
  562.     while ((qrel(q, minr) < 0) || (qrel(q, maxr) > 0)) {
  563.         tmp1 = qsqrt(q, epsilon2);
  564.         qfree(q);
  565.         q = tmp1;
  566.         shift++;
  567.     }
  568.     qfree(minr);
  569.     qfree(maxr);
  570.     /*
  571.      * Calculate a value which will always converge using the formula:
  572.      *    ln((1+x)/(1-x)) = ln(1+x) - ln(1-x).
  573.      */
  574.     tmp1 = qdec(q);
  575.     tmp2 = qinc(q);
  576.     term = qdiv(tmp1, tmp2);
  577.     qfree(tmp1);
  578.     qfree(tmp2);
  579.     qfree(q);
  580.     neg = 0;
  581.     if (qisneg(term)) {
  582.         neg = 1;
  583.         tmp1 = qabs(term);
  584.         qfree(term);
  585.         term = tmp1;
  586.     }
  587.     /*
  588.      * Now use the Taylor series expansion to calculate the result.
  589.      */
  590.     n = 1;
  591.     term2 = qsquare(term);
  592.     sum = qlink(term);
  593.     while (qrel(term, epsilon2) > 0) {
  594.         n += 2;
  595.         tmp1 = qmul(term, term2);
  596.         qfree(term);
  597.         term = qbround(tmp1, bits2);
  598.         qfree(tmp1);
  599.         tmp1 = qdivi(term, (long) n);
  600.         tmp2 = qadd(sum, tmp1);
  601.         qfree(tmp1);
  602.         qfree(sum);
  603.         sum = qbround(tmp2, bits2);
  604.     }
  605.     qfree(epsilon2);
  606.     qfree(term);
  607.     qfree(term2);
  608.     if (neg) {
  609.         tmp1 = qneg(sum);
  610.         qfree(sum);
  611.         sum = tmp1;
  612.     }
  613.     /*
  614.      * Calculate the final result by multiplying by the proper power of two
  615.      * to undo the square roots done at the top.
  616.      */
  617.     tmp1 = qscale(sum, shift);
  618.     qfree(sum);
  619.     sum = qbround(tmp1, bits);
  620.     qfree(tmp1);
  621.     return sum;
  622. }
  623.  
  624.  
  625. /*
  626.  * Calculate the result of raising one number to the power of another.
  627.  * The result is calculated to within the specified relative error.
  628.  */
  629. NUMBER *
  630. qpower(q1, q2, epsilon)
  631.     NUMBER *q1, *q2, *epsilon;
  632. {
  633.     NUMBER *tmp1, *tmp2, *epsilon2;
  634.  
  635.     if (qisint(q2))
  636.         return qpowi(q1, q2);
  637.     epsilon2 = qscale(epsilon, -4L);
  638.     tmp1 = qln(q1, epsilon2);
  639.     tmp2 = qmul(tmp1, q2);
  640.     qfree(tmp1);
  641.     tmp1 = qexp(tmp2, epsilon);
  642.     qfree(tmp2);
  643.     qfree(epsilon2);
  644.     return tmp1;
  645. }
  646.  
  647.  
  648. /*
  649.  * Calculate the Kth root of a number to within the specified accuracy.
  650.  */
  651. NUMBER *
  652. qroot(q1, q2, epsilon)
  653.     NUMBER *q1, *q2, *epsilon;
  654. {
  655.     NUMBER *tmp1, *tmp2, *epsilon2;
  656.     int neg;
  657.  
  658.     if (qisneg(q2) || qiszero(q2) || qisfrac(q2))
  659.         error("Taking bad root of number");
  660.     if (qiszero(q1) || qisone(q1) || qisone(q2))
  661.         return qlink(q1);
  662.     if (qistwo(q2))
  663.         return qsqrt(q1, epsilon);
  664.     neg = qisneg(q1);
  665.     if (neg) {
  666.         if (iseven(q2->num))
  667.             error("Taking even root of negative number");
  668.         q1 = qabs(q1);
  669.     }
  670.     epsilon2 = qscale(epsilon, -4L);
  671.     tmp1 = qln(q1, epsilon2);
  672.     tmp2 = qdiv(tmp1, q2);
  673.     qfree(tmp1);
  674.     tmp1 = qexp(tmp2, epsilon);
  675.     qfree(tmp2);
  676.     qfree(epsilon2);
  677.     if (neg) {
  678.         tmp2 = qneg(tmp1);
  679.         qfree(tmp1);
  680.         tmp1 = tmp2;
  681.     }
  682.     return tmp1;
  683. }
  684.  
  685.  
  686. /*
  687.  * Calculate the hyperbolic cosine function with a relative accuracy less
  688.  * than epsilon.  This is defined by:
  689.  *    cosh(x) = (exp(x) + exp(-x)) / 2.
  690.  */
  691. NUMBER *
  692. qcosh(q, epsilon)
  693.     NUMBER *q, *epsilon;
  694. {
  695.     long scale;
  696.     FULL n;
  697.     FULL m;
  698.     long bits, bits2;
  699.     NUMBER *sum, *term, *qs, *epsilon2, *tmp;
  700.  
  701.     if (qisneg(epsilon) || qiszero(epsilon))
  702.         error("Illegal epsilon value for exp");
  703.     if (qiszero(q))
  704.         return qlink(&_qone_);
  705.     epsilon = qscale(epsilon, -4L);
  706.     /*
  707.      * If the argument is larger than one, then divide it by a power of two
  708.      * so that it is one or less.  This will make the series converge quickly.
  709.      * We will extrapolate the result for the original argument afterwards.
  710.      */
  711.     qs = qabs(q);
  712.     scale = zhighbit(q->num) - zhighbit(q->den) + 1;
  713.     if (scale < 0)
  714.         scale = 0;
  715.     if (scale > 0) {
  716.         if (scale > 100000)
  717.             error("Very large argument for exp");
  718.         tmp = qscale(qs, -scale);
  719.         qfree(qs);
  720.         qs = tmp;
  721.         tmp = qscale(epsilon, -scale);
  722.         qfree(epsilon);
  723.         epsilon = tmp;
  724.     }
  725.     epsilon2 = qscale(epsilon, -4L);
  726.     bits = qprecision(epsilon) + 1;
  727.     bits2 = bits + 10;
  728.     qfree(epsilon);
  729.     tmp = qsquare(qs);
  730.     qfree(qs);
  731.     qs = tmp;
  732.     /*
  733.      * Now use the Taylor series expansion to calculate the exponential.
  734.      * Keep using approximations so that the fractions don't get too large.
  735.      */
  736.     sum = qlink(&_qone_);
  737.     term = qlink(&_qone_);
  738.     n = 0;
  739.     while (qrel(term, epsilon2) > 0) {
  740.         m = ++n;
  741.         m *= ++n;
  742.         tmp = qmul(term, qs);
  743.         qfree(term);
  744.         term = qdivi(tmp, (long) m);
  745.         qfree(tmp);
  746.         tmp = qbround(term, bits2);
  747.         qfree(term);
  748.         term = tmp;
  749.         tmp = qadd(sum, term);
  750.         qfree(sum);
  751.         sum = qbround(tmp, bits2);
  752.         qfree(tmp);
  753.     }
  754.     qfree(term);
  755.     qfree(qs);
  756.     qfree(epsilon2);
  757.     /*
  758.      * Now bring the number back up into range to get the final result.
  759.      * This uses the formula:
  760.      *    cosh(2 * x) = 2 * cosh(x)^2 - 1.
  761.      */
  762.     while (--scale >= 0) {
  763.         tmp = qsquare(sum);
  764.         qfree(sum);
  765.         sum = qscale(tmp, 1L);
  766.         qfree(tmp);
  767.         tmp = qdec(sum);
  768.         qfree(sum);
  769.         sum = qbround(tmp, bits2);
  770.         qfree(tmp);
  771.     }
  772.     tmp = qbround(sum, bits);
  773.     qfree(sum);
  774.     return tmp;
  775. }
  776.  
  777.  
  778. /*
  779.  * Calculate the hyperbolic sine with an accurary less than epsilon.
  780.  * This is calculated using the formula:
  781.  *    cosh(x)^2 - sinh(x)^2 = 1.
  782.  */
  783. NUMBER *
  784. qsinh(q, epsilon)
  785.     NUMBER *q, *epsilon;
  786. {
  787.     NUMBER *tmp1, *tmp2;
  788.  
  789.     if (qisneg(epsilon) || qiszero(epsilon))
  790.         error("Illegal epsilon value for sinh");
  791.     if (qiszero(q))
  792.         return qlink(q);
  793.     epsilon = qscale(epsilon, -4L);
  794.     tmp1 = qcosh(q, epsilon);
  795.     tmp2 = qsquare(tmp1);
  796.     qfree(tmp1);
  797.     tmp1 = qdec(tmp2);
  798.     qfree(tmp2);
  799.     tmp2 = qsqrt(tmp1, epsilon);
  800.     qfree(tmp1);
  801.     if (qisneg(q)) {
  802.         tmp1 = qneg(tmp2);
  803.         qfree(tmp2);
  804.         tmp2 = tmp1;
  805.     }
  806.     qfree(epsilon);
  807.     return tmp2;
  808. }
  809.  
  810.  
  811. /*
  812.  * Calculate the hyperbolic tangent with an accurary less than epsilon.
  813.  * This is calculated using the formula:
  814.  *    tanh(x) = sinh(x) / cosh(x).
  815.  */
  816. NUMBER *
  817. qtanh(q, epsilon)
  818.     NUMBER *q, *epsilon;
  819. {
  820.     NUMBER *tmp1, *tmp2, *coshval;
  821.  
  822.     if (qisneg(epsilon) || qiszero(epsilon))
  823.         error("Illegal epsilon value for tanh");
  824.     if (qiszero(q))
  825.         return qlink(q);
  826.     epsilon = qscale(epsilon, -4L);
  827.     coshval = qcosh(q, epsilon);
  828.     tmp2 = qsquare(coshval);
  829.     tmp1 = qdec(tmp2);
  830.     qfree(tmp2);
  831.     tmp2 = qsqrt(tmp1, epsilon);
  832.     qfree(tmp1);
  833.     if (qisneg(q)) {
  834.         tmp1 = qneg(tmp2);
  835.         qfree(tmp2);
  836.         tmp2 = tmp1;
  837.     }
  838.     qfree(epsilon);
  839.     tmp1 = qdiv(tmp2, coshval);
  840.     qfree(tmp2);
  841.     qfree(coshval);
  842.     return tmp1;
  843. }
  844.  
  845.  
  846. /*
  847.  * Compute the hyperbolic arccosine within the specified accuracy.
  848.  * This is calculated using the formula:
  849.  *    acosh(x) = ln(x + sqrt(x^2 - 1)).
  850.  */
  851. NUMBER *
  852. qacosh(q, epsilon)
  853.     NUMBER *q, *epsilon;
  854. {
  855.     NUMBER *tmp1, *tmp2, *epsilon2;
  856.  
  857.     if (qisneg(epsilon) || qiszero(epsilon))
  858.         error("Illegal epsilon value for acosh");
  859.     if (qisone(q))
  860.         return qlink(&_qzero_);
  861.     if (qreli(q, 1L) < 0)
  862.         error("Argument less than one for acosh");
  863.     epsilon2 = qscale(epsilon, -8L);
  864.     tmp1 = qsquare(q);
  865.     tmp2 = qdec(tmp1);
  866.     qfree(tmp1);
  867.     tmp1 = qsqrt(tmp2, epsilon2);
  868.     qfree(tmp2);
  869.     tmp2 = qadd(tmp1, q);
  870.     qfree(tmp1);
  871.     tmp1 = qln(tmp2, epsilon);
  872.     qfree(tmp2);
  873.     qfree(epsilon2);
  874.     return tmp1;
  875. }
  876.  
  877.  
  878. /*
  879.  * Compute the hyperbolic arcsine within the specified accuracy.
  880.  * This is calculated using the formula:
  881.  *    asinh(x) = ln(x + sqrt(x^2 + 1)).
  882.  */
  883. NUMBER *
  884. qasinh(q, epsilon)
  885.     NUMBER *q, *epsilon;
  886. {
  887.     NUMBER *tmp1, *tmp2, *epsilon2;
  888.  
  889.     if (qisneg(epsilon) || qiszero(epsilon))
  890.         error("Illegal epsilon value for asinh");
  891.     if (qiszero(q))
  892.         return qlink(&_qzero_);
  893.     epsilon2 = qscale(epsilon, -8L);
  894.     tmp1 = qsquare(q);
  895.     tmp2 = qinc(tmp1);
  896.     qfree(tmp1);
  897.     tmp1 = qsqrt(tmp2, epsilon2);
  898.     qfree(tmp2);
  899.     tmp2 = qadd(tmp1, q);
  900.     qfree(tmp1);
  901.     tmp1 = qln(tmp2, epsilon);
  902.     qfree(tmp2);
  903.     qfree(epsilon2);
  904.     return tmp1;
  905. }
  906.  
  907.  
  908. /*
  909.  * Compute the hyperbolic arctangent within the specified accuracy.
  910.  * This is calculated using the formula:
  911.  *    atanh(x) = ln((1 + u) / (1 - u)) / 2.
  912.  */
  913. NUMBER *
  914. qatanh(q, epsilon)
  915.     NUMBER *q, *epsilon;
  916. {
  917.     NUMBER *tmp1, *tmp2, *tmp3;
  918.  
  919.     if (qisneg(epsilon) || qiszero(epsilon))
  920.         error("Illegal epsilon value for atanh");
  921.     if (qiszero(q))
  922.         return qlink(&_qzero_);
  923.     if ((qreli(q, 1L) > 0) || (qreli(q, -1L) < 0))
  924.         error("Argument not between -1 and 1 for atanh");
  925.     tmp1 = qinc(q);
  926.     tmp2 = qsub(&_qone_, q);
  927.     tmp3 = qdiv(tmp1, tmp2);
  928.     qfree(tmp1);
  929.     qfree(tmp2);
  930.     tmp1 = qln(tmp3, epsilon);
  931.     qfree(tmp3);
  932.     tmp2 = qscale(tmp1, -1L);
  933.     qfree(tmp1);
  934.     return tmp2;
  935. }
  936.  
  937. /* END CODE */
  938.