# 4.4BSD/usr/src/contrib/calc-1.26.4/qtrans.c

```/*
* Copyright (c) 1993 David I. Bell
* Permission is granted to use, distribute, or modify this source,
* provided that this copyright notice remains intact.
*
* Transcendental functions for real numbers.
* These are sin, cos, exp, ln, power, cosh, sinh.
*/

#include "math.h"

BOOL _sinisneg_;	/* whether sin(x) < 0 (set by cos(x)) */

/*
* Calculate the cosine of a number with an accuracy within epsilon.
* This also saves the sign of the corresponding sin function.
*/
NUMBER *
qcos(q, epsilon)
NUMBER *q, *epsilon;
{
NUMBER *term, *sum, *qsq, *epsilon2, *tmp;
FULL n, i;
long scale, bits, bits2;

_sinisneg_ = qisneg(q);
if (qisneg(epsilon) || qiszero(epsilon))
error("Illegal epsilon value for cosine");
if (qiszero(q))
bits = qprecision(epsilon) + 1;
epsilon = qscale(epsilon, -4L);
/*
* If the argument is larger than one, then divide it by a power of two
* so that it is one or less.  This will make the series converge quickly.
* We will extrapolate the result for the original argument afterwards.
*/
scale = zhighbit(q->num) - zhighbit(q->den) + 1;
if (scale < 0)
scale = 0;
if (scale > 0) {
q = qscale(q, -scale);
tmp = qscale(epsilon, -scale);
qfree(epsilon);
epsilon = tmp;
}
epsilon2 = qscale(epsilon, -4L);
qfree(epsilon);
bits2 = qprecision(epsilon2) + 10;
/*
* Now use the Taylor series expansion to calculate the cosine.
* Keep using approximations so that the fractions don't get too large.
*/
qsq = qsquare(q);
if (scale > 0)
qfree(q);
n = 0;
while (qrel(term, epsilon2) > 0) {
i = ++n;
i *= ++n;
tmp = qmul(term, qsq);
qfree(term);
term = qdivi(tmp, (long) i);
qfree(tmp);
tmp = qbround(term, bits2);
qfree(term);
term = tmp;
if (n & 2)
tmp = qsub(sum, term);
else
qfree(sum);
sum = qbround(tmp, bits2);
qfree(tmp);
}
qfree(term);
qfree(qsq);
qfree(epsilon2);
/*
* Now scale back up to the original value of x by using the formula:
*	cos(2 * x) = 2 * (cos(x) ^ 2) - 1.
*/
while (--scale >= 0) {
if (qisneg(sum))
_sinisneg_ = !_sinisneg_;
tmp = qsquare(sum);
qfree(sum);
sum = qscale(tmp, 1L);
qfree(tmp);
tmp = qdec(sum);
qfree(sum);
sum = qbround(tmp, bits2);
qfree(tmp);
}
tmp = qbround(sum, bits);
qfree(sum);
return tmp;
}

/*
* Calculate the sine of a number with an accuracy within epsilon.
* This is calculated using the formula:
*	sin(x)^2 + cos(x)^2 = 1.
* The only tricky bit is resolving the sign of the result.
* Future: Use sin(3*x) = 3*sin(x) - 4*sin(x)^3.
*/
NUMBER *
qsin(q, epsilon)
NUMBER *q, *epsilon;
{
NUMBER *tmp1, *tmp2, *epsilon2;

if (qisneg(epsilon) || qiszero(epsilon))
error("Illegal epsilon value for sine");
if (qiszero(q))
epsilon2 = qscale(epsilon, -4L);
tmp1 = qcos(q, epsilon2);
qfree(epsilon2);
tmp2 = qlegtoleg(tmp1, epsilon, _sinisneg_);
qfree(tmp1);
return tmp2;
}

/*
* Calculate the tangent function.
*/
NUMBER *
qtan(q, epsilon)
NUMBER *q, *epsilon;
{
NUMBER *cosval, *sinval, *epsilon2, *tmp, *res;

if (qisneg(epsilon) || qiszero(epsilon))
error("Illegal epsilon value for tangent");
if (qiszero(q))
epsilon2 = qscale(epsilon, -4L);
cosval = qcos(q, epsilon2);
sinval = qlegtoleg(cosval, epsilon2, _sinisneg_);
qfree(epsilon2);
tmp = qdiv(sinval, cosval);
qfree(cosval);
qfree(sinval);
res = qbround(tmp, qprecision(epsilon) + 1);
qfree(tmp);
return res;
}

/*
* Calculate the arcsine function.
* The result is in the range -pi/2 to pi/2.
*/
NUMBER *
qasin(q, epsilon)
NUMBER *q, *epsilon;
{
NUMBER *sum, *term, *epsilon2, *qsq, *tmp;
FULL n, i;
long bits, bits2;
int neg;
NUMBER mulnum;
HALF numval[2];
HALF denval[2];

if (qisneg(epsilon) || qiszero(epsilon))
error("Illegal epsilon value for arcsine");
if (qiszero(q))
if ((qrel(q, &_qone_) > 0) || (qrel(q, &_qnegone_) < 0))
error("Argument too large for asin");
neg = qisneg(q);
q = qabs(q);
epsilon = qscale(epsilon, -4L);
epsilon2 = qscale(epsilon, -4L);
mulnum.num.sign = 0;
mulnum.num.len = 1;
mulnum.num.v = numval;
mulnum.den.sign = 0;
mulnum.den.len = 1;
mulnum.den.v = denval;
/*
* If the argument is too near one (we use .5) then reduce the
* argument to a more accurate range using the formula:
*	asin(x) = 2 * asin(sqrt((1 - sqrt(1 - x^2)) / 2)).
*/
if (qrel(q, &_qonehalf_) > 0) {
sum = qlegtoleg(q, epsilon2, FALSE);
qfree(q);
tmp = qsub(&_qone_, sum);
qfree(sum);
sum = qscale(tmp, -1L);
qfree(tmp);
tmp = qsqrt(sum, epsilon2);
qfree(sum);
qfree(epsilon2);
sum = qasin(tmp, epsilon);
qfree(tmp);
qfree(epsilon);
tmp = qscale(sum, 1L);
qfree(sum);
if (neg) {
sum = qneg(tmp);
qfree(tmp);
tmp = sum;
}
return tmp;
}
/*
* Argument is between zero and .5, so use the series.
*/
epsilon = qscale(epsilon, -4L);
epsilon2 = qscale(epsilon, -4L);
bits = qprecision(epsilon) + 1;
bits2 = bits + 10;
qsq = qsquare(q);
qfree(q);
n = 1;
while (qrel(term, epsilon2) > 0) {
i = n * n;
numval[0] = i & BASE1;
if (i >= BASE) {
numval[1] = i / BASE;
mulnum.den.len = 2;
}
i = (n + 1) * (n + 2);
denval[0] = i & BASE1;
if (i >= BASE) {
denval[1] = i / BASE;
mulnum.den.len = 2;
}
tmp = qmul(term, qsq);
qfree(term);
term = qmul(tmp, &mulnum);
qfree(tmp);
tmp = qbround(term, bits2);
qfree(term);
term = tmp;
qfree(sum);
sum = qbround(tmp, bits2);
qfree(tmp);
n += 2;
}
qfree(epsilon);
qfree(epsilon2);
qfree(term);
qfree(qsq);
tmp = qbround(sum, bits);
qfree(sum);
if (neg) {
term = qneg(tmp);
qfree(tmp);
tmp = term;
}
return tmp;
}

/*
* Calculate the acos function.
* The result is in the range 0 to pi.
*/
NUMBER *
qacos(q, epsilon)
NUMBER *q, *epsilon;
{
NUMBER *tmp1, *tmp2, *tmp3, *epsilon2;

if (qisneg(epsilon) || qiszero(epsilon))
error("Illegal epsilon value for arccosine");
if (qisone(q))
if ((qrel(q, &_qone_) > 0) || (qrel(q, &_qnegone_) < 0))
error("Argument too large for acos");
/*
* Calculate the result using the formula:
*	acos(x) = asin(sqrt(1 - x^2)).
* The formula is only good for positive x, so we must fix up the
* result for negative values.
*/
epsilon2 = qscale(epsilon, -8L);
tmp1 = qlegtoleg(q, epsilon2, FALSE);
qfree(epsilon2);
tmp2 = qasin(tmp1, epsilon);
qfree(tmp1);
if (!qisneg(q))
return tmp2;
/*
* For negative values, we need to subtract the asin from pi.
*/
tmp1 = qpi(epsilon);
tmp3 = qsub(tmp1, tmp2);
qfree(tmp1);
qfree(tmp2);
tmp1 = qbround(tmp3, qprecision(epsilon) + 1);
qfree(tmp3);
return tmp1;
}

/*
* Calculate the arctangent function with a accuracy less than epsilon.
* This uses the formula:
*	atan(x) = asin(sqrt(x^2 / (x^2 + 1))).
*/
NUMBER *
qatan(q, epsilon)
NUMBER *q, *epsilon;
{
NUMBER *tmp1, *tmp2, *tmp3, *epsilon2;

if (qisneg(epsilon) || qiszero(epsilon))
error("Illegal epsilon value for arctangent");
if (qiszero(q))
tmp1 = qsquare(q);
tmp2 = qinc(tmp1);
tmp3 = qdiv(tmp1, tmp2);
qfree(tmp1);
qfree(tmp2);
epsilon2 = qscale(epsilon, -8L);
tmp1 = qsqrt(tmp3, epsilon2);
qfree(epsilon2);
qfree(tmp3);
tmp2 = qasin(tmp1, epsilon);
qfree(tmp1);
if (qisneg(q)) {
tmp1 = qneg(tmp2);
qfree(tmp2);
tmp2 = tmp1;
}
return tmp2;
}

/*
* Calculate the angle which is determined by the point (x,y).
* This is the same as arctan for non-negative x, but gives the correct
* value for negative x.  By convention, y is the first argument.
* For example, qatan2(1, -1) = 3/4 * pi.
*/
NUMBER *
qatan2(qy, qx, epsilon)
NUMBER *qy, *qx, *epsilon;
{
NUMBER *tmp1, *tmp2, *epsilon2;

if (qisneg(epsilon) || qiszero(epsilon))
error("Illegal epsilon value for atan2");
if (qiszero(qy) && qiszero(qx))
error("Zero coordinates for atan2");
/*
* If the point is on the negative real axis, then the answer is pi.
*/
if (qiszero(qy) && qisneg(qx))
return qpi(epsilon);
/*
* If the point is in the right half plane, then use the normal atan.
*/
if (!qisneg(qx) && !qiszero(qx)) {
if (qiszero(qy))
tmp1 = qdiv(qy, qx);
tmp2 = qatan(tmp1, epsilon);
qfree(tmp1);
return tmp2;
}
/*
* The point is in the left half plane.  Calculate the angle by finding
* the atan of half the angle using the formula:
*	atan2(y,x) = 2 * atan((sqrt(x^2 + y^2) - x) / y).
*/
epsilon2 = qscale(epsilon, -4L);
tmp1 = qhypot(qx, qy, epsilon2);
tmp2 = qsub(tmp1, qx);
qfree(tmp1);
tmp1 = qdiv(tmp2, qy);
qfree(tmp2);
tmp2 = qatan(tmp1, epsilon2);
qfree(tmp1);
qfree(epsilon2);
tmp1 = qscale(tmp2, 1L);
qfree(tmp2);
return tmp1;
}

/*
* Calculate the value of pi to within the required epsilon.
* This uses the following formula which only needs integer calculations
* except for the final operation:
*	pi = 1 / SUMOF(comb(2 * N, N) ^ 3 * (42 * N + 5) / 2 ^ (12 * N + 4)),
* where the summation runs from N=0.  This formula gives about 6 bits of
* accuracy per term.  Since the denominator for each term is a power of two,
* we can simply use shifts to sum the terms.  The combinatorial numbers
* in the formula are calculated recursively using the formula:
*	comb(2*(N+1), N+1) = 2 * comb(2 * N, N) * (2 * N + 1) / N.
*/
NUMBER *
qpi(epsilon)
NUMBER *epsilon;
{
ZVALUE comb;			/* current combinatorial value */
ZVALUE sum;			/* current sum */
ZVALUE tmp1, tmp2;
NUMBER *r, *t1, qtmp;
long shift;			/* current shift of result */
long N;				/* current term number */
long bits;			/* needed number of bits of precision */
long t;

if (qiszero(epsilon) || qisneg(epsilon))
bits = qprecision(epsilon) + 4;
comb = _one_;
itoz(5L, &sum);
N = 0;
shift = 4;
do {
t = 1 + (++N & 0x1);
(void) zdivi(comb, N / (3 - t), &tmp1);
freeh(comb.v);
zmuli(tmp1, t * (2 * N - 1), &comb);
freeh(tmp1.v);
zsquare(comb, &tmp1);
zmul(comb, tmp1, &tmp2);
freeh(tmp1.v);
zmuli(tmp2, 42 * N + 5, &tmp1);
freeh(tmp2.v);
zshift(sum, 12L, &tmp2);
freeh(sum.v);
t = zhighbit(tmp1);
freeh(tmp1.v);
freeh(tmp2.v);
shift += 12;
} while ((shift - t) < bits);
qtmp.num = _one_;
qtmp.den = sum;
t1 = qscale(&qtmp, shift);
freeh(sum.v);
r = qbround(t1, bits);
qfree(t1);
return r;
}

/*
* Calculate the exponential function with a relative accuracy less than
* epsilon.
*/
NUMBER *
qexp(q, epsilon)
NUMBER *q, *epsilon;
{
long scale;
FULL n;
long bits, bits2;
NUMBER *sum, *term, *qs, *epsilon2, *tmp;

if (qisneg(epsilon) || qiszero(epsilon))
error("Illegal epsilon value for exp");
if (qiszero(q))
epsilon = qscale(epsilon, -4L);
/*
* If the argument is larger than one, then divide it by a power of two
* so that it is one or less.  This will make the series converge quickly.
* We will extrapolate the result for the original argument afterwards.
* Also make the argument non-negative.
*/
qs = qabs(q);
scale = zhighbit(q->num) - zhighbit(q->den) + 1;
if (scale < 0)
scale = 0;
if (scale > 0) {
if (scale > 100000)
error("Very large argument for exp");
tmp = qscale(qs, -scale);
qfree(qs);
qs = tmp;
tmp = qscale(epsilon, -scale);
qfree(epsilon);
epsilon = tmp;
}
epsilon2 = qscale(epsilon, -4L);
bits = qprecision(epsilon) + 1;
bits2 = bits + 10;
qfree(epsilon);
/*
* Now use the Taylor series expansion to calculate the exponential.
* Keep using approximations so that the fractions don't get too large.
*/
n = 0;
while (qrel(term, epsilon2) > 0) {
n++;
tmp = qmul(term, qs);
qfree(term);
term = qdivi(tmp, (long) n);
qfree(tmp);
tmp = qbround(term, bits2);
qfree(term);
term = tmp;
qfree(sum);
sum = qbround(tmp, bits2);
qfree(tmp);
}
qfree(term);
qfree(qs);
qfree(epsilon2);
/*
* Now repeatedly square the answer to get the final result.
* Then invert it if the original argument was negative.
*/
while (--scale >= 0) {
tmp = qsquare(sum);
qfree(sum);
sum = qbround(tmp, bits2);
qfree(tmp);
}
tmp = qbround(sum, bits);
qfree(sum);
if (qisneg(q)) {
sum = qinv(tmp);
qfree(tmp);
tmp = sum;
}
return tmp;
}

/*
* Calculate the natural logarithm of a number accurate to the specified
* epsilon.
*/
NUMBER *
qln(q, epsilon)
NUMBER *q, *epsilon;
{
NUMBER *term, *term2, *sum, *epsilon2, *tmp1, *tmp2;
NUMBER *minr, *maxr;
long shift, bits, bits2;
int neg;
FULL n;

if (qisneg(q) || qiszero(q))
error("log of non-positive number");
if (qisneg(epsilon) || qiszero(epsilon))
error("Illegal epsilon for ln");
epsilon = qscale(epsilon, -4L);
epsilon2 = qscale(epsilon, -4L);
bits = qprecision(epsilon) + 1;
bits2 = bits + 10;
qfree(epsilon);
/*
* Scale down the logarithm to the convergent range by taking square
* roots.  The result will be modified to get the final value later.
*/
minr = iitoq(7L, 10L);
maxr = iitoq(7L, 5L);
shift = 1;
while ((qrel(q, minr) < 0) || (qrel(q, maxr) > 0)) {
tmp1 = qsqrt(q, epsilon2);
qfree(q);
q = tmp1;
shift++;
}
qfree(minr);
qfree(maxr);
/*
* Calculate a value which will always converge using the formula:
*	ln((1+x)/(1-x)) = ln(1+x) - ln(1-x).
*/
tmp1 = qdec(q);
tmp2 = qinc(q);
term = qdiv(tmp1, tmp2);
qfree(tmp1);
qfree(tmp2);
qfree(q);
neg = 0;
if (qisneg(term)) {
neg = 1;
tmp1 = qabs(term);
qfree(term);
term = tmp1;
}
/*
* Now use the Taylor series expansion to calculate the result.
*/
n = 1;
term2 = qsquare(term);
while (qrel(term, epsilon2) > 0) {
n += 2;
tmp1 = qmul(term, term2);
qfree(term);
term = qbround(tmp1, bits2);
qfree(tmp1);
tmp1 = qdivi(term, (long) n);
qfree(tmp1);
qfree(sum);
sum = qbround(tmp2, bits2);
}
qfree(epsilon2);
qfree(term);
qfree(term2);
if (neg) {
tmp1 = qneg(sum);
qfree(sum);
sum = tmp1;
}
/*
* Calculate the final result by multiplying by the proper power of two
* to undo the square roots done at the top.
*/
tmp1 = qscale(sum, shift);
qfree(sum);
sum = qbround(tmp1, bits);
qfree(tmp1);
return sum;
}

/*
* Calculate the result of raising one number to the power of another.
* The result is calculated to within the specified relative error.
*/
NUMBER *
qpower(q1, q2, epsilon)
NUMBER *q1, *q2, *epsilon;
{
NUMBER *tmp1, *tmp2, *epsilon2;

if (qisint(q2))
return qpowi(q1, q2);
epsilon2 = qscale(epsilon, -4L);
tmp1 = qln(q1, epsilon2);
tmp2 = qmul(tmp1, q2);
qfree(tmp1);
tmp1 = qexp(tmp2, epsilon);
qfree(tmp2);
qfree(epsilon2);
return tmp1;
}

/*
* Calculate the Kth root of a number to within the specified accuracy.
*/
NUMBER *
qroot(q1, q2, epsilon)
NUMBER *q1, *q2, *epsilon;
{
NUMBER *tmp1, *tmp2, *epsilon2;
int neg;

if (qisneg(q2) || qiszero(q2) || qisfrac(q2))
if (qiszero(q1) || qisone(q1) || qisone(q2))
if (qistwo(q2))
return qsqrt(q1, epsilon);
neg = qisneg(q1);
if (neg) {
if (iseven(q2->num))
error("Taking even root of negative number");
q1 = qabs(q1);
}
epsilon2 = qscale(epsilon, -4L);
tmp1 = qln(q1, epsilon2);
tmp2 = qdiv(tmp1, q2);
qfree(tmp1);
tmp1 = qexp(tmp2, epsilon);
qfree(tmp2);
qfree(epsilon2);
if (neg) {
tmp2 = qneg(tmp1);
qfree(tmp1);
tmp1 = tmp2;
}
return tmp1;
}

/*
* Calculate the hyperbolic cosine function with a relative accuracy less
* than epsilon.  This is defined by:
*	cosh(x) = (exp(x) + exp(-x)) / 2.
*/
NUMBER *
qcosh(q, epsilon)
NUMBER *q, *epsilon;
{
long scale;
FULL n;
FULL m;
long bits, bits2;
NUMBER *sum, *term, *qs, *epsilon2, *tmp;

if (qisneg(epsilon) || qiszero(epsilon))
error("Illegal epsilon value for exp");
if (qiszero(q))
epsilon = qscale(epsilon, -4L);
/*
* If the argument is larger than one, then divide it by a power of two
* so that it is one or less.  This will make the series converge quickly.
* We will extrapolate the result for the original argument afterwards.
*/
qs = qabs(q);
scale = zhighbit(q->num) - zhighbit(q->den) + 1;
if (scale < 0)
scale = 0;
if (scale > 0) {
if (scale > 100000)
error("Very large argument for exp");
tmp = qscale(qs, -scale);
qfree(qs);
qs = tmp;
tmp = qscale(epsilon, -scale);
qfree(epsilon);
epsilon = tmp;
}
epsilon2 = qscale(epsilon, -4L);
bits = qprecision(epsilon) + 1;
bits2 = bits + 10;
qfree(epsilon);
tmp = qsquare(qs);
qfree(qs);
qs = tmp;
/*
* Now use the Taylor series expansion to calculate the exponential.
* Keep using approximations so that the fractions don't get too large.
*/
n = 0;
while (qrel(term, epsilon2) > 0) {
m = ++n;
m *= ++n;
tmp = qmul(term, qs);
qfree(term);
term = qdivi(tmp, (long) m);
qfree(tmp);
tmp = qbround(term, bits2);
qfree(term);
term = tmp;
qfree(sum);
sum = qbround(tmp, bits2);
qfree(tmp);
}
qfree(term);
qfree(qs);
qfree(epsilon2);
/*
* Now bring the number back up into range to get the final result.
* This uses the formula:
*	cosh(2 * x) = 2 * cosh(x)^2 - 1.
*/
while (--scale >= 0) {
tmp = qsquare(sum);
qfree(sum);
sum = qscale(tmp, 1L);
qfree(tmp);
tmp = qdec(sum);
qfree(sum);
sum = qbround(tmp, bits2);
qfree(tmp);
}
tmp = qbround(sum, bits);
qfree(sum);
return tmp;
}

/*
* Calculate the hyperbolic sine with an accurary less than epsilon.
* This is calculated using the formula:
*	cosh(x)^2 - sinh(x)^2 = 1.
*/
NUMBER *
qsinh(q, epsilon)
NUMBER *q, *epsilon;
{
NUMBER *tmp1, *tmp2;

if (qisneg(epsilon) || qiszero(epsilon))
error("Illegal epsilon value for sinh");
if (qiszero(q))
epsilon = qscale(epsilon, -4L);
tmp1 = qcosh(q, epsilon);
tmp2 = qsquare(tmp1);
qfree(tmp1);
tmp1 = qdec(tmp2);
qfree(tmp2);
tmp2 = qsqrt(tmp1, epsilon);
qfree(tmp1);
if (qisneg(q)) {
tmp1 = qneg(tmp2);
qfree(tmp2);
tmp2 = tmp1;
}
qfree(epsilon);
return tmp2;
}

/*
* Calculate the hyperbolic tangent with an accurary less than epsilon.
* This is calculated using the formula:
*	tanh(x) = sinh(x) / cosh(x).
*/
NUMBER *
qtanh(q, epsilon)
NUMBER *q, *epsilon;
{
NUMBER *tmp1, *tmp2, *coshval;

if (qisneg(epsilon) || qiszero(epsilon))
error("Illegal epsilon value for tanh");
if (qiszero(q))
epsilon = qscale(epsilon, -4L);
coshval = qcosh(q, epsilon);
tmp2 = qsquare(coshval);
tmp1 = qdec(tmp2);
qfree(tmp2);
tmp2 = qsqrt(tmp1, epsilon);
qfree(tmp1);
if (qisneg(q)) {
tmp1 = qneg(tmp2);
qfree(tmp2);
tmp2 = tmp1;
}
qfree(epsilon);
tmp1 = qdiv(tmp2, coshval);
qfree(tmp2);
qfree(coshval);
return tmp1;
}

/*
* Compute the hyperbolic arccosine within the specified accuracy.
* This is calculated using the formula:
*	acosh(x) = ln(x + sqrt(x^2 - 1)).
*/
NUMBER *
qacosh(q, epsilon)
NUMBER *q, *epsilon;
{
NUMBER *tmp1, *tmp2, *epsilon2;

if (qisneg(epsilon) || qiszero(epsilon))
error("Illegal epsilon value for acosh");
if (qisone(q))
if (qreli(q, 1L) < 0)
error("Argument less than one for acosh");
epsilon2 = qscale(epsilon, -8L);
tmp1 = qsquare(q);
tmp2 = qdec(tmp1);
qfree(tmp1);
tmp1 = qsqrt(tmp2, epsilon2);
qfree(tmp2);
qfree(tmp1);
tmp1 = qln(tmp2, epsilon);
qfree(tmp2);
qfree(epsilon2);
return tmp1;
}

/*
* Compute the hyperbolic arcsine within the specified accuracy.
* This is calculated using the formula:
*	asinh(x) = ln(x + sqrt(x^2 + 1)).
*/
NUMBER *
qasinh(q, epsilon)
NUMBER *q, *epsilon;
{
NUMBER *tmp1, *tmp2, *epsilon2;

if (qisneg(epsilon) || qiszero(epsilon))
error("Illegal epsilon value for asinh");
if (qiszero(q))
epsilon2 = qscale(epsilon, -8L);
tmp1 = qsquare(q);
tmp2 = qinc(tmp1);
qfree(tmp1);
tmp1 = qsqrt(tmp2, epsilon2);
qfree(tmp2);
qfree(tmp1);
tmp1 = qln(tmp2, epsilon);
qfree(tmp2);
qfree(epsilon2);
return tmp1;
}

/*
* Compute the hyperbolic arctangent within the specified accuracy.
* This is calculated using the formula:
*	atanh(x) = ln((1 + u) / (1 - u)) / 2.
*/
NUMBER *
qatanh(q, epsilon)
NUMBER *q, *epsilon;
{
NUMBER *tmp1, *tmp2, *tmp3;

if (qisneg(epsilon) || qiszero(epsilon))
error("Illegal epsilon value for atanh");
if (qiszero(q))
if ((qreli(q, 1L) > 0) || (qreli(q, -1L) < 0))
error("Argument not between -1 and 1 for atanh");
tmp1 = qinc(q);
tmp2 = qsub(&_qone_, q);
tmp3 = qdiv(tmp1, tmp2);
qfree(tmp1);
qfree(tmp2);
tmp1 = qln(tmp3, epsilon);
qfree(tmp3);
tmp2 = qscale(tmp1, -1L);
qfree(tmp1);
return tmp2;
}

/* END CODE */
```