drygdryg/calc
1
0
Fork 0
mirror of https://github.com/lcn2/calc.git synced 2025-08-16 01:03:29 +03:00
Code Issues Packages Projects Releases Wiki Activity
Files
master
calc/qtrans.c
Landon Curt Noll db77e29a23 convert ASCII TABs to ASCII SPACEs 
Converted all ASCII tabs to ASCII spaces using a 8 character
tab stop, for all files, except for all Makefiles (plus rpm.mk).
The `git diff -w` reports no changes.
2024-07-11 22:03:52 -07:00

3694 lines
96 KiB
C
Raw Permalink Blame History

/*
* qtrans - transcendental functions for real numbers
*
* Copyright (C) 1999-2007,2021-2023 David I. Bell, Landon Curt Noll and Ernest Bowen
*
* Primary author: David I. Bell
*
* Calc is open software; you can redistribute it and/or modify it under
* the terms of the version 2.1 of the GNU Lesser General Public License
* as published by the Free Software Foundation.
*
* Calc is distributed in the hope that it will be useful, but WITHOUT
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
* or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General
* Public License for more details.
*
* A copy of version 2.1 of the GNU Lesser General Public License is
* distributed with calc under the filename COPYING-LGPL. You should have
* received a copy with calc; if not, write to Free Software Foundation, Inc.
* 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
*
* Under source code control: 1990/02/15 01:48:22
* File existed as early as: before 1990
*
* Share and enjoy! :-) http://www.isthe.com/chongo/tech/comp/calc/
*/
/*
* Transcendental functions for real numbers.
* These are sin, cos, exp, ln, power, cosh, sinh.
*/
#include "qmath.h"
#include "config.h"
#include "errtbl.h"
#include "banned.h" /* include after system header <> includes */
HALF _qlgenum_[] = { 36744 };
HALF _qlgeden_[] = { 25469 };
NUMBER _qlge_ = { { _qlgenum_, 1, 0 }, { _qlgeden_, 1, 0 }, 1, NULL };
/*
* cache the natural logarithm of 10 and 2
*/
STATIC NUMBER *ln_10 = NULL;
STATIC NUMBER *ln_10_epsilon = NULL;
STATIC NUMBER *ln_2 = NULL;
STATIC NUMBER *ln_2_epsilon = NULL;
STATIC NUMBER *ln_n = NULL;
STATIC NUMBER *ln_n_epsilon = NULL;
/*
* cache pi
*
* pivalue[LAST_PI_EPSILON] - last epsilon used to calculate pi
* pivalue[LAST_PI_VALUE] - last calculated pi
* given pivalue[LAST_PI_EPSILON] epsilon
* pivalue[LAST_PI_DIV_180_EPSILON] - last epsilon used to calculate pi/180
* pivalue[LAST_PI_DIV_180_VALUE] - last calculated pi/180 given
* pivalue[LAST_PI_DIV_180_EPSILON] epsilon
* pivalue[LAST_PI_DIV_200_EPSILON] - last epsilon used to calculate pi/200
* pivalue[LAST_PI_DIV_200_VALUE] - last calculated pi/200 given
* pivalue[LAST_PI_DIV_200_EPSILON] epsilon
*/
enum pi_cache {
LAST_PI_EPSILON = 0,
LAST_PI_VALUE,
LAST_PI_DIV_180_EPSILON,
LAST_PI_DIV_180_VALUE,
LAST_PI_DIV_200_EPSILON,
LAST_PI_DIV_200_VALUE,
PI_CACHE_LEN /* must be last */
};
STATIC NUMBER *pivalue[PI_CACHE_LEN] = {
NULL, /* LAST_PI_EPSILON */
NULL, /* LAST_PI_VALUE */
NULL, /* LAST_PI_DIV_180_EPSILON */
NULL, /* LAST_PI_DIV_180_VALUE */
NULL, /* LAST_PI_DIV_200_EPSILON */
NULL, /* LAST_PI_DIV_200_VALUE */
};
/*
* other static function declarations
*/
STATIC NUMBER *qexprel(NUMBER *q, long bitnum);
/*
* Evaluate and store in specified locations the sin and cos of a given
* number to accuracy corresponding to a specified number of binary digits.
*/
void
qsincos(NUMBER *q, long bitnum, NUMBER **vs, NUMBER **vc)
{
long n, m, k, h, s, t, d;
NUMBER *qtmp1, *qtmp2;
ZVALUE X, cossum, sinsum, mul, ztmp1, ztmp2, ztmp3;
qtmp1 = qqabs(q);
h = qilog2(qtmp1);
qfree(qtmp1);
k = bitnum + h + 1;
if (k < 0) {
*vs = qlink(&_qzero_);
*vc = qlink(&_qone_);
return;
}
s = k;
if (k) {
do {
t = s;
s = (s + k/s)/2;
}
while (t > s);
} /* s is int(sqrt(k)) */
s++;
if (s < -h)
s = -h;
n = h + s; /* n is number of squaring that will be required */
m = bitnum + n;
while (s > 0) { /* increasing m by ilog2(s) */
s >>= 1;
m++;
} /* m is working number of bits */
qtmp1 = qscale(q, m - n);
zquo(qtmp1->num, qtmp1->den, &X, conf->triground);
qfree(qtmp1);
if (ziszero(X)) {
zfree(X);
*vs = qlink(&_qzero_);
*vc = qlink(&_qone_);
return;
}
zbitvalue(m, &cossum);
zcopy(X, &sinsum);
zcopy(X, &mul);
d = 1;
for (;;) {
X.sign = !X.sign;
zmul(X, mul, &ztmp1);
zfree(X);
zshift(ztmp1, -m, &ztmp2);
zfree(ztmp1);
zdivi(ztmp2, ++d, &X);
zfree(ztmp2);
if (ziszero(X))
break;
zadd(cossum, X, &ztmp1);
zfree(cossum);
cossum = ztmp1;
zmul(X, mul, &ztmp1);
zfree(X);
zshift(ztmp1, -m, &ztmp2);
zfree(ztmp1);
zdivi(ztmp2, ++d, &X);
zfree(ztmp2);
if (ziszero(X))
break;
zadd(sinsum, X, &ztmp1);
zfree(sinsum);
sinsum = ztmp1;
}
zfree(X);
zfree(mul);
while (n-- > 0) {
zsquare(cossum, &ztmp1);
zsquare(sinsum, &ztmp2);
zsub(ztmp1, ztmp2, &ztmp3);
zfree(ztmp1);
zfree(ztmp2);
zmul(cossum, sinsum, &ztmp1);
zfree(cossum);
zfree(sinsum);
zshift(ztmp3, -m, &cossum);
zfree(ztmp3);
zshift(ztmp1, 1 - m, &sinsum);
zfree(ztmp1);
}
h = zlowbit(cossum);
qtmp1 = qalloc();
if (m > h) {
zshift(cossum, -h, &qtmp1->num);
zbitvalue(m - h, &qtmp1->den);
} else {
zshift(cossum, - m, &qtmp1->num);
}
zfree(cossum);
*vc = qtmp1;
h = zlowbit(sinsum);
qtmp2 = qalloc();
if (m > h) {
zshift(sinsum, -h, &qtmp2->num);
zbitvalue(m - h, &qtmp2->den);
} else {
zshift(sinsum, -m, &qtmp2->num);
}
zfree(sinsum);
*vs = qtmp2;
return;
}
/*
* Calculate the cosine of a number to a near multiple of epsilon.
* This calls qsincos() and discards the value of sin.
*/
NUMBER *
qcos(NUMBER *q, NUMBER *epsilon)
{
NUMBER *sin, *cos, *res;
long n;
if (qiszero(epsilon)) {
math_error("Zero epsilon value for cosine");
not_reached();
}
if (qiszero(q))
return qlink(&_qone_);
n = -qilog2(epsilon);
if (n < 0)
return qlink(&_qzero_);
qsincos(q, n + 2, &sin, &cos);
qfree(sin);
res = qmappr(cos, epsilon, conf->triground);
qfree(cos);
return res;
}
/*
* This calls qsincos() and discards the value of cos.
*/
NUMBER *
qsin(NUMBER *q, NUMBER *epsilon)
{
NUMBER *sin, *cos, *res;
long n;
if (qiszero(epsilon)) {
math_error("Zero epsilon value for sine");
not_reached();
}
n = -qilog2(epsilon);
if (qiszero(q) || n < 0)
return qlink(&_qzero_);
qsincos(q, n + 2, &sin, &cos);
qfree(cos);
res = qmappr(sin, epsilon, conf->triground);
qfree(sin);
return res;
}
/*
* Calculate the tangent function.
*/
NUMBER *
qtan(NUMBER *q, NUMBER *epsilon)
{
NUMBER *sin, *cos, *tan, *res;
long n, k, m;
if (qiszero(epsilon)) {
math_error("Zero epsilon value for tangent");
not_reached();
}
if (qiszero(q))
return qlink(q);
n = qilog2(epsilon);
k = (n > 0) ? 4 + n/2 : 4;
for (;;) {
qsincos(q, 2 * k - n, &sin, &cos);
if (qiszero(cos)) {
qfree(sin);
qfree(cos);
k = 2 * k - n + 4;
continue;
}
m = -qilog2(cos);
if (m < k)
break;
qfree(sin);
qfree(cos);
k = m + 1;
}
tan = qqdiv(sin, cos);
qfree(sin);
qfree(cos);
res = qmappr(tan, epsilon, conf->triground);
qfree(tan);
return res;
}
/*
* Calculate the cotangent function.
*/
NUMBER *
qcot(NUMBER *q, NUMBER *epsilon)
{
NUMBER *sin, *cos, *cot, *res;
long n, k, m;
if (qiszero(epsilon)) {
math_error("Zero epsilon value for cotangent");
not_reached();
}
if (qiszero(q)) {
math_error("Zero argument for cotangent");
not_reached();
}
k = -qilog2(q);
n = qilog2(epsilon);
if (k < 0)
k = (n > 0) ? n/2 : 0;
k += 4;
for (;;) {
qsincos(q, 2 * k - n, &sin, &cos);
if (qiszero(sin)) {
qfree(sin);
qfree(cos);
k = 2 * k - n + 4;
continue;
}
m = -qilog2(sin);
if (m < k)
break;
qfree(sin);
qfree(cos);
k = m + 1;
}
cot = qqdiv(cos, sin);
qfree(sin);
qfree(cos);
res = qmappr(cot, epsilon, conf->triground);
qfree(cot);
return res;
}
/*
* Calculate the secant function.
*/
NUMBER *
qsec(NUMBER *q, NUMBER *epsilon)
{
NUMBER *sin, *cos, *sec, *res;
long n, k, m;
if (qiszero(epsilon)) {
math_error("Zero epsilon value for secant");
not_reached();
}
if (qiszero(q))
return qlink(&_qone_);
n = qilog2(epsilon);
k = (n > 0) ? 4 + n/2 : 4;
for (;;) {
qsincos(q, 2 * k - n, &sin, &cos);
qfree(sin);
if (qiszero(cos)) {
qfree(cos);
k = 2 * k - n + 4;
continue;
}
m = -qilog2(cos);
if (m < k)
break;
qfree(cos);
k = m + 1;
}
sec = qinv(cos);
qfree(cos);
res = qmappr(sec, epsilon, conf->triground);
qfree(sec);
return res;
}
/*
* Calculate the cosecant function.
*/
NUMBER *
qcsc(NUMBER *q, NUMBER *epsilon)
{
NUMBER *sin, *cos, *csc, *res;
long n, k, m;
if (qiszero(epsilon)) {
math_error("Zero epsilon value for cosecant");
not_reached();
}
if (qiszero(q)) {
math_error("Zero argument for cosecant");
not_reached();
}
k = -qilog2(q);
n = qilog2(epsilon);
if (k < 0)
k = (n > 0) ? n/2 : 0;
k += 4;
for (;;) {
qsincos(q, 2 * k - n, &sin, &cos);
qfree(cos);
if (qiszero(sin)) {
qfree(sin);
k = 2 * k - n + 4;
continue;
}
m = -qilog2(sin);
if (m < k)
break;
qfree(sin);
k = m + 1;
}
csc = qinv(sin);
qfree(sin);
res = qmappr(csc, epsilon, conf->triground);
qfree(csc);
return res;
}
/*
* Calculate the arcsine function.
* The result is in the range -pi/2 to pi/2.
*/
NUMBER *
qasin(NUMBER *q, NUMBER *epsilon)
{
NUMBER *qtmp1, *qtmp2, *epsilon1;
ZVALUE ztmp;
bool neg;
FLAG r;
if (qiszero(epsilon)) {
math_error("Zero epsilon value for asin");
not_reached();
}
if (qiszero(q)) {
return qlink(&_qzero_);
}
ztmp = q->num;
neg = ztmp.sign;
ztmp.sign = 0;
r = zrel(ztmp, q->den);
if (r > 0) {
return NULL;
}
if (r == 0) {
epsilon1 = qscale(epsilon, 1L);
qtmp2 = qpi(epsilon1);
qtmp1 = qscale(qtmp2, -1L);
} else {
epsilon1 = qscale(epsilon, -2L);
qtmp1 = qalloc();
zsquare(q->num, &qtmp1->num);
zsquare(q->den, &ztmp);
zsub(ztmp, qtmp1->num, &qtmp1->den);
zfree(ztmp);
qtmp2 = qsqrt(qtmp1, epsilon1, conf->triground);
qfree(qtmp1);
qtmp1 = qatan(qtmp2, epsilon);
}
qfree(qtmp2);
qfree(epsilon1);
if (neg) {
qtmp2 = qneg(qtmp1);
qfree(qtmp1);
return(qtmp2);
}
return qtmp1;
}
/*
* Calculate the acos function.
* The result is in the range 0 to pi.
*/
NUMBER *
qacos(NUMBER *q, NUMBER *epsilon)
{
NUMBER *q1, *q2, *epsilon1;
ZVALUE z;
if (qiszero(epsilon)) {
math_error("Zero epsilon value for acos");
not_reached();
}
if (qisone(q))
return qlink(&_qzero_);
if (qisnegone(q))
return qpi(epsilon);
z = q->num;
z.sign = 0;
if (zrel(z, q->den) > 0)
return NULL;
epsilon1 = qscale(epsilon, -3L); /* ??? */
q1 = qalloc();
zsub(q->den, q->num, &q1->num);
zadd(q->den, q->num, &q1->den);
q2 = qsqrt(q1, epsilon1, conf->triground);
qfree(q1);
qfree(epsilon1);
epsilon1 = qscale(epsilon, -1L);
q1 = qatan(q2, epsilon1);
qfree(epsilon1);
qfree(q2);
q2 = qscale(q1, 1L);
qfree(q1)
return q2;
}
/*
* Calculate the arctangent function to the nearest or next to nearest
* multiple of epsilon. Algorithm uses
* atan(x) = 2 * atan(x/(1 + sqrt(1+x^2)))
* to reduce x to a small value and then
* atan(x) = x - x^3/3 + ...
*/
NUMBER *
qatan(NUMBER *q, NUMBER *epsilon)
{
long m, k, i;
FULL d;
ZVALUE X, D, DD, sum, mul, term, ztmp1, ztmp2;
NUMBER *qtmp, *res;
bool sign;
if (qiszero(epsilon)) {
math_error("Zero epsilon value for arctangent");
not_reached();
}
if (qiszero(q))
return qlink(&_qzero_);
m = 12 - qilog2(epsilon);
/* 4 bits for 4 doublings; 8 for rounding */
if (m < 8)
m = 8; /* m is number of working binary digits */
qtmp = qscale(q, m);
zquo(qtmp->num, qtmp->den, &X, conf->triground);
qfree(qtmp);
zbitvalue(m, &D); /* q has become X/D */
zsquare(D, &DD);
i = 4; /* maybe this should be larger */
while (i-- > 0 && !ziszero(X)) {
zsquare(X, &ztmp1);
zadd(ztmp1, DD, &ztmp2);
zfree(ztmp1);
zsqrt(ztmp2, &ztmp1, conf->triground);
zfree(ztmp2);
zadd(ztmp1, D, &ztmp2);
zfree(ztmp1);
zshift(X, m, &ztmp1);
zfree(X);
zquo(ztmp1, ztmp2, &X, conf->triground);
zfree(ztmp1);
zfree(ztmp2);
}
zfree(DD);
zfree(D);
if (ziszero(X)) {
zfree(X);
return qlink(&_qzero_);
}
zcopy(X, &sum);
zsquare(X, &ztmp1);
zshift(ztmp1, -m, &mul);
zfree(ztmp1);
d = 3;
sign = !X.sign;
for (;;) {
if (d > BASE) {
math_error("Too many terms required for atan");
not_reached();
}
zmul(X, mul, &ztmp1);
zfree(X);
zshift(ztmp1, -m, &X); /* X now (original X)^d */
zfree(ztmp1);
zdivi(X, d, &term);
if (ziszero(term)) {
zfree(term);
break;
}
term.sign = sign;
zadd(sum, term, &ztmp1);
zfree(sum);
zfree(term);
sum = ztmp1;
sign = !sign;
d += 2;
}
zfree(mul);
zfree(X);
qtmp = qalloc();
k = zlowbit(sum);
if (k) {
zshift(sum, -k, &qtmp->num);
zfree(sum);
} else {
qtmp->num = sum;
}
zbitvalue(m - 4 - k, &qtmp->den);
res = qmappr(qtmp, epsilon, conf->triground);
qfree(qtmp);
return res;
}
/*
* Inverse secant function
*/
NUMBER *
qasec(NUMBER *q, NUMBER *epsilon)
{
NUMBER *tmp, *res;
tmp = qinv(q);
res = qacos(tmp, epsilon);
qfree(tmp);
return res;
}
/*
* Inverse cosecant function
*/
NUMBER *
qacsc(NUMBER *q, NUMBER *epsilon)
{
NUMBER *tmp, *res;
tmp = qinv(q);
res = qasin(tmp, epsilon);
qfree(tmp);
return res;
}
/*
* Inverse cotangent function
*/
NUMBER *
qacot(NUMBER *q, NUMBER *epsilon)
{
NUMBER *tmp1, *tmp2, *tmp3, *epsilon1;
if (qiszero(epsilon)) {
math_error("Zero epsilon for acot");
not_reached();
}
if (qiszero(q)) {
epsilon1 = qscale(epsilon, 1L);
tmp1 = qpi(epsilon1);
qfree(epsilon1);
tmp2 = qscale(tmp1, -1L);
qfree(tmp1);
return tmp2;
}
tmp1 = qinv(q);
if (!qisneg(q)) {
tmp2 = qatan(tmp1, epsilon);
qfree(tmp1);
return tmp2;
}
epsilon1 = qscale(epsilon, -2L);
tmp2 = qatan(tmp1, epsilon1);
qfree(tmp1);
tmp1 = qpi(epsilon1);
qfree(epsilon1);
tmp3 = qqadd(tmp1, tmp2);
qfree(tmp1);
qfree(tmp2);
tmp1 = qmappr(tmp3, epsilon, conf->triground);
qfree(tmp3);
return tmp1;
}
/*
* Calculate the angle which is determined by the point (x,y).
* This is the same as atan(y/x) for positive x, but is continuous
* except for y = 0, x <= 0. By convention, y is the first argument.
* For all x, y, -pi < atan2 <= pi. For example, qatan2(1, -1) = 3/4 * pi.
*/
NUMBER *
qatan2(NUMBER *qy, NUMBER *qx, NUMBER *epsilon)
{
NUMBER *tmp1, *tmp2, *tmp3, *epsilon2;
if (qiszero(epsilon)) {
math_error("Zero epsilon value for atan2");
not_reached();
}
if (qiszero(qy) && qiszero(qx)) {
/* conform to 4.3BSD ANSI/IEEE 754-1985 math lib */
return qlink(&_qzero_);
}
/*
* 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))
return qlink(&_qzero_);
tmp1 = qqdiv(qy, qx);
tmp2 = qatan(tmp1, epsilon);
qfree(tmp1);
return tmp2;
}
/*
* The point is in the left half plane (x <= 0) with nonzero y.
* Calculate the angle by using the formula:
* atan2(y,x) = 2 * atan(sgn(y) * sqrt((x/y)^2 + 1) - x/y).
*/
epsilon2 = qscale(epsilon, -4L);
tmp1 = qqdiv(qx, qy);
tmp2 = qsquare(tmp1);
tmp3 = qqadd(tmp2, &_qone_);
qfree(tmp2);
tmp2 = qsqrt(tmp3, epsilon2, 24L | (qy->num.sign * 64));
qfree(tmp3);
tmp3 = qsub(tmp2, tmp1);
qfree(tmp2);
qfree(tmp1);
qfree(epsilon2);
epsilon2 = qscale(epsilon, -1L);
tmp1 = qatan(tmp3, epsilon2);
qfree(epsilon2);
qfree(tmp3);
tmp2 = qscale(tmp1, 1L);
qfree(tmp1);
return tmp2;
}
/*
* 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(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;
/* firewall */
if (qiszero(epsilon)) {
math_error("zero epsilon value for pi");
not_reached();
}
/* use pi cache if epsilon marches, else flush if needed */
if (pivalue[LAST_PI_EPSILON] != NULL &&
pivalue[LAST_PI_VALUE] != NULL &&
epsilon == pivalue[LAST_PI_EPSILON]) {
return qlink(pivalue[LAST_PI_VALUE]);
}
if (pivalue[LAST_PI_EPSILON] != NULL) {
qfree(pivalue[LAST_PI_EPSILON]);
}
if (pivalue[LAST_PI_VALUE] != NULL) {
qfree(pivalue[LAST_PI_VALUE]);
}
bits = -qilog2(epsilon) + 4;
if (bits < 4)
bits = 4;
comb = _one_;
itoz(5L, &sum);
N = 0;
shift = 4;
do {
t = 1 + (++N & 0x1);
(void) zdivi(comb, N / (3 - t), &tmp1);
zfree(comb);
zmuli(tmp1, t * (2 * N - 1), &comb);
zfree(tmp1);
zsquare(comb, &tmp1);
zmul(comb, tmp1, &tmp2);
zfree(tmp1);
zmuli(tmp2, 42 * N + 5, &tmp1);
zfree(tmp2);
zshift(sum, 12L, &tmp2);
zfree(sum);
zadd(tmp1, tmp2, &sum);
t = zhighbit(tmp1);
zfree(tmp1);
zfree(tmp2);
shift += 12;
} while ((shift - t) < bits);
zfree(comb);
qtmp.num = _one_;
qtmp.den = sum;
t1 = qscale(&qtmp, shift);
zfree(sum);
r = qmappr(t1, epsilon, conf->triground);
qfree(t1);
pivalue[LAST_PI_EPSILON] = qlink(epsilon);
pivalue[LAST_PI_VALUE] = qlink(r);
return r;
}
/*
* qpidiv180 - calculate pi / 180
*
* This function returns pi/180 as used to covert between radians and degrees.
*/
NUMBER *
qpidiv180(NUMBER *epsilon)
{
NUMBER *pi, *pidiv180;
/* firewall */
if (qiszero(epsilon)) {
math_error("zero epsilon value for qpidiv180");
not_reached();
}
/* use pi/180 cache if epsilon marches, else flush if needed */
if (pivalue[LAST_PI_DIV_180_EPSILON] != NULL &&
pivalue[LAST_PI_DIV_180_VALUE] != NULL &&
epsilon == pivalue[LAST_PI_DIV_180_EPSILON]) {
return qlink(pivalue[LAST_PI_DIV_180_VALUE]);
}
if (pivalue[LAST_PI_DIV_180_EPSILON] != NULL) {
qfree(pivalue[LAST_PI_DIV_180_EPSILON]);
}
if (pivalue[LAST_PI_DIV_180_VALUE] != NULL) {
qfree(pivalue[LAST_PI_DIV_180_VALUE]);
}
/* let qpi() returned cached pi or calculate new as needed */
pi = qpi(epsilon);
/* calculate pi/180 */
pidiv180 = qdivi(pi, 180);
/* cache epsilon and pi/180 */
pivalue[LAST_PI_DIV_180_EPSILON] = qlink(epsilon);
pivalue[LAST_PI_DIV_180_VALUE] = qlink(pidiv180);
/* return pi/180 */
return pidiv180;
}
/*
* qpidiv200 - calculate pi / 200
*
* This function returns pi/200 as used to covert between radians and gradians.
*/
NUMBER *
qpidiv200(NUMBER *epsilon)
{
NUMBER *pi, *pidiv200;
/* firewall */
if (qiszero(epsilon)) {
math_error("zero epsilon value for qpidiv200");
not_reached();
}
/* use pi/200 cache if epsilon marches, else flush if needed */
if (pivalue[LAST_PI_DIV_200_EPSILON] != NULL &&
pivalue[LAST_PI_DIV_200_VALUE] != NULL &&
epsilon == pivalue[LAST_PI_DIV_200_EPSILON]) {
return qlink(pivalue[LAST_PI_DIV_200_VALUE]);
}
if (pivalue[LAST_PI_DIV_200_EPSILON] != NULL) {
qfree(pivalue[LAST_PI_DIV_200_EPSILON]);
}
if (pivalue[LAST_PI_DIV_200_VALUE] != NULL) {
qfree(pivalue[LAST_PI_DIV_200_VALUE]);
}
/* let qpi() returned cached pi or calculate new as needed */
pi = qpi(epsilon);
/* calculate pi/200 */
pidiv200 = qdivi(pi, 200);
/* cache epsilon and pi/200 */
pivalue[LAST_PI_DIV_200_EPSILON] = qlink(epsilon);
pivalue[LAST_PI_DIV_200_VALUE] = qlink(pidiv200);
/* return pi/200 */
return pidiv200;
}
/*
* Calculate the exponential function to the nearest or next to nearest
* multiple of the positive number epsilon.
*/
NUMBER *
qexp(NUMBER *q, NUMBER *epsilon)
{
long m, n;
NUMBER *tmp1, *tmp2;
if (qiszero(epsilon)) {
math_error("Zero epsilon value for exp");
not_reached();
}
if (qiszero(q))
return qlink(&_qone_);
tmp1 = qmul(q, &_qlge_);
m = qtoi(tmp1); /* exp(q) < 2^(m+1) or m == MAXLONG */
qfree(tmp1);
if (m > (1 << 30))
return NULL;
n = qilog2(epsilon); /* 2^n <= epsilon < 2^(n+1) */
if (m < n)
return qlink(&_qzero_);
tmp1 = qqabs(q);
tmp2 = qexprel(tmp1, m - n + 1);
qfree(tmp1);
if (tmp2 == NULL)
return NULL;
if (qisneg(q)) {
tmp1 = qinv(tmp2);
qfree(tmp2);
tmp2 = tmp1;
}
tmp1 = qmappr(tmp2, epsilon, conf->triground);
qfree(tmp2);
return tmp1;
}
/*
* Calculate the exponential function with relative error corresponding
* to a specified number of significant bits
* Requires *q >= 0, bitnum >= 0.
* This returns NULL if more than 2^30 working bits would be required.
*/
S_FUNC NUMBER *
qexprel(NUMBER *q, long bitnum)
{
long n, m, k, h, s, t, d;
NUMBER *qtmp1;
ZVALUE X, B, sum, term, ztmp1, ztmp2;
h = qilog2(q);
k = bitnum + h + 1;
if (k < 0)
return qlink(&_qone_);
s = k;
if (k) {
do {
t = s;
s = (s + k/s)/2;
}
while (t > s);
} /* s is int(sqrt(k)) */
s++;
if (s < -h)
s = -h;
n = h + s; /* n is number of squarings that will be required */
m = bitnum + n;
if (m > (1 << 30))
return NULL;
while (s > 0) { /* increasing m by ilog2(s) */
s >>= 1;
m++;
} /* m is working number of bits */
qtmp1 = qscale(q, m - n);
zquo(qtmp1->num, qtmp1->den, &X, conf->triground);
qfree(qtmp1);
if (ziszero(X)) {
zfree(X);
return qlink(&_qone_);
}
zbitvalue(m, &sum);
zcopy(X, &term);
d = 1;
do {
zadd(sum, term, &ztmp1);
zfree(sum);
sum = ztmp1;
zmul(term, X, &ztmp1);
zfree(term);
zshift(ztmp1, -m, &ztmp2);
zfree(ztmp1);
zdivi(ztmp2, ++d, &term);
zfree(ztmp2);
}
while (!ziszero(term));
zfree(term);
zfree(X);
k = 0;
zbitvalue(2 * m + 1, &B);
while (n-- > 0) {
k *= 2;
zsquare(sum, &ztmp1);
zfree(sum);
if (zrel(ztmp1, B) >= 0) {
zshift(ztmp1, -m - 1, &sum);
k++;
} else {
zshift(ztmp1, -m, &sum);
}
zfree(ztmp1);
}
zfree(B);
h = zlowbit(sum);
qtmp1 = qalloc();
if (m > h + k) {
zshift(sum, -h, &qtmp1->num);
zbitvalue(m - h - k, &qtmp1->den);
}
else
zshift(sum, k - m, &qtmp1->num);
zfree(sum);
return qtmp1;
}
/*
* Calculate the natural logarithm of a number accurate to the specified
* positive epsilon.
*/
NUMBER *
qln(NUMBER *q, NUMBER *epsilon)
{
long m, n, k, h, d;
ZVALUE term, sum, mul, pow, X, D, B, ztmp;
NUMBER *qtmp, *res;
bool neg;
if (qiszero(q)) {
math_error("ln of 0");
not_reached();
}
if (qiszero(epsilon)) {
math_error("Zero epsilon value for ln");
not_reached();
}
if (qisunit(q))
return qlink(&_qzero_);
q = qqabs(q); /* Ignore sign of q */
neg = (zrel(q->num, q->den) < 0);
if (neg) {
qtmp = qinv(q);
qfree(q);
q = qtmp;
}
k = qilog2(q);
m = -qilog2(epsilon); /* m will be number of working bits */
if (m < 0)
m = 0;
h = k;
while (h > 0) {
h /= 2;
m++; /* Add 1 for each sqrt until X < 2 */
}
m += 18; /* 8 more sqrts, 8 for rounding, 2 for epsilon/4 */
qtmp = qscale(q, m - k);
zquo(qtmp->num, qtmp->den, &X, conf->triground);
qfree(q);
qfree(qtmp);
zbitvalue(m, &D); /* Now "q" = X/D */
zbitvalue(m - 8, &ztmp);
zadd(D, ztmp, &B); /* Will take sqrts until X <= B */
zfree(ztmp);
n = 1; /* n is to count 1 + number of sqrts */
while (k > 0 || zrel(X, B) > 0) {
n++;
zshift(X, m + (k & 1), &ztmp);
zfree(X);
zsqrt(ztmp, &X, conf->triground);
zfree(ztmp)
k /= 2;
}
zfree(B);
zsub(X, D, &pow); /* pow, mul used as tmps */
zadd(X, D, &mul);
zfree(X);
zfree(D);
zshift(pow, m, &ztmp);
zfree(pow);
zquo(ztmp, mul, &pow, conf->triground); /* pow now (X - D)/(X + D) */
zfree(ztmp);
zfree(mul);
zcopy(pow, &sum); /* pow is first term of sum */
zsquare(pow, &ztmp);
zshift(ztmp, -m, &mul); /* mul is now multiplier for powers */
zfree(ztmp);
d = 1;
for (;;) {
zmul(pow, mul, &ztmp);
zfree(pow);
zshift(ztmp, -m, &pow);
zfree(ztmp);
d += 2;
zdivi(pow, d, &term); /* Round down div should be round off */
if (ziszero(term)) {
zfree(term);
break;
}
zadd(sum, term, &ztmp);
zfree(term);
zfree(sum);
sum = ztmp;
}
zfree(pow);
zfree(mul);
qtmp = qalloc(); /* qtmp is to be 2^n * sum / 2^m */
k = zlowbit(sum);
sum.sign = neg;
if (k + n >= m) {
zshift(sum, n - m, &qtmp->num);
zfree(sum);
} else {
if (k) {
zshift(sum, -k, &qtmp->num);
zfree(sum);
} else {
qtmp->num = sum;
}
zbitvalue(m - k - n, &qtmp->den);
}
res = qmappr(qtmp, epsilon, conf->triground);
qfree(qtmp);
return res;
}
/*
* Calculate the base 10 logarithm
*
* log(q) = ln(q) / ln(10)
*/
NUMBER *
qlog(NUMBER *q, NUMBER *epsilon)
{
int need_new_ln_10 = true; /* false => use cached ln_10 value */
NUMBER *ln_q; /* ln(x) */
NUMBER *ret; /* base 10 logarithm of x */
/* firewall */
if (qiszero(q)) {
math_error("log of 0");
not_reached();
}
if (qiszero(epsilon)) {
math_error("Zero epsilon value for log");
not_reached();
}
/*
* shortcut for small integer powers of 10
*/
if (qisint(q) && qispos(q) && !zge8192b(q->num) && ziseven(q->num)) {
bool is_10_power; /* true ==> q is a power of 10 */
long ilog_10; /* integer log base 10 */
/* try for a quick small power of 10 log */
ilog_10 = zlog10(q->num, &is_10_power );
if (is_10_power == true) {
/* is small power of 10, return log */
return itoq(ilog_10);
}
/* q is an even integer that is not a power of 10 */
}
/*
* compute ln(c) first
*/
ln_q = qln(q, epsilon);
/* quick return for log(1) == 0 */
if (qiszero(ln_q)) {
return ln_q;
}
/*
* save epsilon for ln(10) if needed
*/
if (ln_10_epsilon == NULL) {
/* first time call */
ln_10_epsilon = qcopy(epsilon);
} else if (qcmp(ln_10_epsilon, epsilon) == true) {
/* replaced cached value with epsilon arg */
qfree(ln_10_epsilon);
ln_10_epsilon = qcopy(epsilon);
} else if (ln_10 != NULL) {
/* the previously computed ln(2) is OK to use */
need_new_ln_10 = false;
}
/*
* compute ln(10) if needed
*/
if (need_new_ln_10 == true) {
if (ln_10 != NULL) {
qfree(ln_10);
}
ln_10 = qln(&_qten_, ln_10_epsilon);
}
/*
* return ln(q) / ln(10)
*/
ret = qqdiv(ln_q, ln_10);
qfree(ln_q);
return ret;
}
/*
* Calculate the base 2 logarithm
*
* log(q) = ln(q) / ln(2)
*/
NUMBER *
qlog2(NUMBER *q, NUMBER *epsilon)
{
int need_new_ln_2 = true; /* false => use cached ln_2 value */
NUMBER *ln_q; /* ln(x) */
NUMBER *ret; /* base 2 logarithm of x */
/* firewall */
if (qiszero(q)) {
math_error("log2 of 0");
not_reached();
}
if (qiszero(epsilon)) {
math_error("Zero epsilon value for log2");
not_reached();
}
/*
* special case: q is integer power of 2
*/
ret = qalloc();
if (qispowerof2(q, &ret)) {
return ret;
}
qfree(ret);
/*
* compute ln(c) first
*/
ln_q = qln(q, epsilon);
/* quick return for ln(1) == 0 */
if (qiszero(ln_q)) {
return ln_q;
}
/*
* save epsilon for ln(2) if needed
*/
if (ln_2_epsilon == NULL) {
/* first time call */
ln_2_epsilon = qcopy(epsilon);
} else if (qcmp(ln_2_epsilon, epsilon) == true) {
/* replaced cached value with epsilon arg */
qfree(ln_2_epsilon);
ln_2_epsilon = qcopy(epsilon);
} else if (ln_2 != NULL) {
/* the previously computed ln(2) is OK to use */
need_new_ln_2 = false;
}
/*
* compute ln(2) if needed
*/
if (need_new_ln_2 == true) {
if (ln_2 != NULL) {
qfree(ln_2);
}
ln_2 = qln(&_qtwo_, ln_2_epsilon);
}
/*
* return ln(q) / ln(2)
*/
ret = qqdiv(ln_q, ln_2);
qfree(ln_q);
return ret;
}
/*
* Calculate the base n logarithm
*
* logn(q, n) = ln(q) / ln(n)
*/
NUMBER *
qlogn(NUMBER *q, NUMBER *n, NUMBER *epsilon)
{
int need_new_ln_n = true; /* false => use cached ln_n value */
NUMBER *ln_q; /* ln(x) */
NUMBER *ret; /* base 2 logarithm of x */
/* firewall */
if (qiszero(q)) {
math_error("log base n of 0");
not_reached();
}
if (qiszero(epsilon)) {
math_error("Zero epsilon value for logn");
not_reached();
}
if (qiszero(n)) {
math_error("invalid logarithm base of 0 for logn");
not_reached();
}
if (qisone(n)) {
math_error("invalid logarithm base of 1 for logn");
not_reached();
}
/*
* special case: q is integer power of 2
*/
ret = qalloc();
if (qispowerof2(q, &ret)) {
return ret;
}
/* XXX - deal with n is integer power of 2 - XXX */
qfree(ret);
/*
* compute ln(q) first
*/
ln_q = qln(q, epsilon);
/* quick return for ln(1) == 0 */
if (qiszero(ln_q)) {
return ln_q;
}
/*
* save epsilon for ln(n) if needed
*/
if (ln_n_epsilon == NULL) {
/* first time call */
ln_n_epsilon = qcopy(epsilon);
} else if (qcmp(ln_n_epsilon, epsilon) == true) {
/* replaced cached value with epsilon arg */
qfree(ln_n_epsilon);
ln_n_epsilon = qcopy(epsilon);
} else if (ln_n != NULL) {
/* the previously computed ln(2) is OK to use */
need_new_ln_n = false;
}
/*
* compute ln(n) if needed
*/
if (need_new_ln_n == true) {
if (ln_n != NULL) {
qfree(ln_n);
}
ln_n = qln(&_qtwo_, ln_n_epsilon);
}
/*
* return ln(q) / ln(2)
*/
ret = qqdiv(ln_q, ln_n);
qfree(ln_q);
return ret;
}
/*
* Calculate the result of raising one number to the power of another.
* The result is calculated to the nearest or next to nearest multiple of
* epsilon.
*/
NUMBER *
qpower(NUMBER *q1, NUMBER *q2, NUMBER *epsilon)
{
NUMBER *tmp1, *tmp2, *epsilon2;
NUMBER *q1tmp, *q2tmp;
long m, n;
if (qiszero(epsilon)) {
math_error("Zero epsilon for power");
not_reached();
}
if (qiszero(q1) && qisneg(q2)) {
math_error("Negative power of zero");
not_reached();
}
if (qiszero(q2) || qisone(q1))
return qlink(&_qone_);
if (qiszero(q1))
return qlink(&_qzero_);
if (qisneg(q1)) {
math_error("Negative base for qpower");
not_reached();
}
if (qisone(q2))
return qmappr(q1, epsilon, conf->triground);
if (zrel(q1->num, q1->den) < 0) {
q1tmp = qinv(q1);
q2tmp = qneg(q2);
} else {
q1tmp = qlink(q1);
q2tmp = qlink(q2);
}
if (qisone(q2tmp)) {
qfree(q2tmp);
q2tmp = qmappr(q1tmp, epsilon, conf->triground);
qfree(q1tmp);
return q2tmp;
}
m = qilog2(q1tmp);
n = qilog2(epsilon);
if (qisneg(q2tmp)) {
if (m > 0) {
tmp1 = itoq(m);
tmp2 = qmul(tmp1, q2tmp);
m = qtoi(tmp2);
} else {
tmp1 = qdec(q1tmp);
tmp2 = qqdiv(tmp1, q1tmp);
qfree(tmp1);
tmp1 = qmul(tmp2, q2tmp);
qfree(tmp2);
tmp2 = qmul(tmp1, &_qlge_);
m = qtoi(tmp2);
}
} else {
if (m > 0) {
tmp1 = itoq(m + 1);
tmp2 = qmul(tmp1, q2tmp);
m = qtoi(tmp2);
} else {
tmp1 = qdec(q1tmp);
tmp2 = qmul(tmp1, q2tmp);
qfree(tmp1);
tmp1 = qmul(tmp2, &_qlge_);
m = qtoi(tmp1);
}
}
qfree(tmp1);
qfree(tmp2);
if (m > (1 << 30)) {
qfree(q1tmp);
qfree(q2tmp);
return NULL;
}
m += 1;
if (m < n) {
qfree(q1tmp);
qfree(q2tmp);
return qlink(&_qzero_);
}
tmp1 = qqdiv(epsilon, q2tmp);
tmp2 = qscale(tmp1, -m - 4);
epsilon2 = qqabs(tmp2);
qfree(tmp1);
qfree(tmp2);
tmp1 = qln(q1tmp, epsilon2);
qfree(epsilon2);
tmp2 = qmul(tmp1, q2tmp);
qfree(tmp1);
qfree(q1tmp);
qfree(q2tmp);
if (qisneg(tmp2)) {
tmp1 = qneg(tmp2);
qfree(tmp2);
tmp2 = qexprel(tmp1, m - n + 3);
if (tmp2 == NULL) {
qfree(tmp1);
return NULL;
}
qfree(tmp1);
tmp1 = qinv(tmp2);
} else {
tmp1 = qexprel(tmp2, m - n + 3) ;
if (tmp1 == NULL) {
qfree(tmp2);
return NULL;
}
}
qfree(tmp2);
tmp2 = qmappr(tmp1, epsilon, conf->triground);
qfree(tmp1);
return tmp2;
}
/*
* Calculate the K-th root of a number to within the specified accuracy.
*/
NUMBER *
qroot(NUMBER *q1, NUMBER *q2, NUMBER *epsilon)
{
NUMBER *tmp1, *tmp2;
int neg;
if (qiszero(epsilon)) {
math_error("Zero epsilon for root");
not_reached();
}
if (qisneg(q2) || qiszero(q2) || qisfrac(q2)) {
math_error("Taking bad root of number");
not_reached();
}
if (qiszero(q1) || qisone(q1) || qisone(q2))
return qlink(q1);
if (qistwo(q2))
return qsqrt(q1, epsilon, conf->triground);
neg = qisneg(q1);
if (neg) {
if (ziseven(q2->num)) {
math_error("Taking even root of negative number");
not_reached();
}
q1 = qqabs(q1);
}
tmp2 = qinv(q2);
tmp1 = qpower(q1, tmp2, epsilon);
qfree(tmp2);
if (tmp1 == NULL)
return NULL;
if (neg) {
tmp2 = qneg(tmp1);
qfree(tmp1);
tmp1 = tmp2;
}
return tmp1;
}
/* Calculate the hyperbolic cosine function to the nearest or next to
* nearest multiple of epsilon.
* This is calculated using cosh(x) = (exp(x) + 1/exp(x))/2;
*/
NUMBER *
qcosh(NUMBER *q, NUMBER *epsilon)
{
NUMBER *tmp1, *tmp2, *tmp3, *epsilon1;
epsilon1 = qscale(epsilon, -2);
tmp1 = qqabs(q);
tmp2 = qexp(tmp1, epsilon1);
qfree(tmp1);
qfree(epsilon1);
if (tmp2 == NULL)
return NULL;
tmp1 = qinv(tmp2);
tmp3 = qqadd(tmp1, tmp2);
qfree(tmp1)
qfree(tmp2)
tmp1 = qscale(tmp3, -1);
qfree(tmp3);
tmp2 = qmappr(tmp1, epsilon, conf->triground);
qfree(tmp1);
return tmp2;
}
/*
* Calculate the hyperbolic sine to the nearest or next to nearest
* multiple of epsilon.
* This is calculated using sinh(x) = (exp(x) - 1/exp(x))/2.
*/
NUMBER *
qsinh(NUMBER *q, NUMBER *epsilon)
{
NUMBER *tmp1, *tmp2, *tmp3, *epsilon1;
if (qiszero(q))
return qlink(&_qzero_);
epsilon1 = qscale(epsilon, -3);
tmp1 = qqabs(q);
tmp2 = qexp(tmp1, epsilon1);
qfree(tmp1);
qfree(epsilon1);
if (tmp2 == NULL)
return NULL;
tmp1 = qinv(tmp2);
tmp3 = qispos(q) ? qsub(tmp2, tmp1) : qsub(tmp1, tmp2);
qfree(tmp1)
qfree(tmp2)
tmp1 = qscale(tmp3, -1);
qfree(tmp3);
tmp2 = qmappr(tmp1, epsilon, conf->triground);
qfree(tmp1);
return tmp2;
}
/*
* Calculate the hyperbolic tangent to the nearest or next to nearest
* multiple of epsilon.
* This is calculated using the formulae:
* tanh(x) = 1 or -1 for very large abs(x)
* tanh(x) = (+ or -)(1 - 2 * exp(2 * x)) for large abx(x)
* tanh(x) = (exp(2*x) - 1)/(exp(2*x) + 1) otherwise
*/
NUMBER *
qtanh(NUMBER *q, NUMBER *epsilon)
{
NUMBER *tmp1, *tmp2, *tmp3;
long n, m;
n = qilog2(epsilon);
if (n > 0 || qiszero(q))
return qlink(&_qzero_);
n = -n;
tmp1 = qqabs(q);
tmp2 = qmul(tmp1, &_qlge_);
m = qtoi(tmp2); /* exp(|q|) < 2^(m+1) or m == MAXLONG */
qfree(tmp2);
if (m > 1 + n/2) {
qfree(tmp1);
return q->num.sign ? qlink(&_qnegone_) : qlink(&_qone_);
}
tmp2 = qscale(tmp1, 1);
qfree(tmp1);
tmp1 = qexprel(tmp2, 2 + n);
qfree(tmp2);
if (m > 1 + n/4) {
tmp2 = qqdiv(&_qtwo_, tmp1);
qfree(tmp1);
tmp1 = qsub(&_qone_, tmp2);
qfree(tmp2);
} else {
tmp2 = qdec(tmp1);
tmp3 = qinc(tmp1);
qfree(tmp1);
tmp1 = qqdiv(tmp2, tmp3);
qfree(tmp2);
qfree(tmp3);
}
tmp2 = qmappr(tmp1, epsilon, conf->triground);
qfree(tmp1);
if (qisneg(q)) {
tmp1 = qneg(tmp2);
qfree(tmp2);
return tmp1;
}
return tmp2;
}
/*
* Hyperbolic cotangent.
* Calculated using coth(x) = 1 + 2/(exp(2*x) - 1)
*/
NUMBER *
qcoth(NUMBER *q, NUMBER *epsilon)
{
NUMBER *tmp1, *tmp2, *res;
long n, k;
if (qiszero(epsilon)) {
math_error("Zero epsilon value for coth");
not_reached();
}
if (qiszero(q)) {
math_error("Zero argument for coth");
not_reached();
}
tmp1 = qscale(q, 1);
tmp2 = qqabs(tmp1);
qfree(tmp1);
k = qilog2(tmp2);
n = qilog2(epsilon);
if (k > 0) {
tmp1 = qmul(&_qlge_, tmp2);
k = qtoi(tmp1);
qfree(tmp1);
} else {
k = 2 * k;
}
k = 4 - k - n;
if (k < 4)
k = 4;
tmp1 = qexprel(tmp2, k);
qfree(tmp2);
tmp2 = qdec(tmp1);
qfree(tmp1);
if (qiszero(tmp2)) {
math_error("This should not happen ????");
not_reached();
}
tmp1 = qinv(tmp2);
qfree(tmp2);
tmp2 = qscale(tmp1, 1);
qfree(tmp1);
tmp1 = qinc(tmp2);
qfree(tmp2);
if (qisneg(q)) {
tmp2 = qneg(tmp1);
qfree(tmp1);
tmp1 = tmp2;
}
res = qmappr(tmp1, epsilon, conf->triground);
qfree(tmp1);
return res;
}
NUMBER *
qsech(NUMBER *q, NUMBER *epsilon)
{
NUMBER *tmp1, *tmp2, *tmp3, *res;
long n, k;
if (qiszero(epsilon)) {
math_error("Zero epsilon value for sech");
not_reached();
}
if (qiszero(q))
return qlink(&_qone_);
tmp1 = qqabs(q);
k = 0;
if (zrel(tmp1->num, tmp1->den) >= 0) {
tmp2 = qmul(&_qlge_, tmp1);
k = qtoi(tmp2);
qfree(tmp2);
}
n = qilog2(epsilon);
if (k + n > 1) {
qfree(tmp1);
return qlink(&_qzero_);
}
tmp2 = qexprel(tmp1, 4 - k - n);
qfree(tmp1);
tmp1 = qinv(tmp2);
tmp3 = qqadd(tmp1, tmp2);
qfree(tmp1);
qfree(tmp2);
tmp1 = qinv(tmp3);
qfree(tmp3);
tmp2 = qscale(tmp1, 1);
qfree(tmp1);
res = qmappr(tmp2, epsilon, conf->triground);
qfree(tmp2);
return res;
}
NUMBER *
qcsch(NUMBER *q, NUMBER *epsilon)
{
NUMBER *tmp1, *tmp2, *tmp3, *res;
long n, k;
if (qiszero(epsilon)) {
math_error("Zero epsilon value for csch");
not_reached();
}
if (qiszero(q)) {
math_error("Zero argument for csch");
not_reached();
}
n = qilog2(epsilon);
tmp1 = qqabs(q);
if (zrel(tmp1->num, tmp1->den) >= 0) {
tmp2 = qmul(&_qlge_, tmp1);
k = qtoi(tmp2);
qfree(tmp2);
} else {
k = 2 * qilog2(tmp1);
}
if (k + n >= 1) {
qfree(tmp1);
return qlink(&_qzero_);
}
tmp2 = qexprel(tmp1, 4 - k - n);
qfree(tmp1);
tmp1 = qinv(tmp2);
if (qisneg(q))
tmp3 = qsub(tmp1, tmp2);
else
tmp3 = qsub(tmp2, tmp1);
qfree(tmp1);
qfree(tmp2);
tmp1 = qinv(tmp3);
qfree(tmp3)
tmp2 = qscale(tmp1, 1);
qfree(tmp1);
res = qmappr(tmp2, epsilon, conf->triground);
qfree(tmp2);
return res;
}
/*
* Compute the hyperbolic arccosine within the specified accuracy.
* This is calculated using the formula:
* acosh(x) = ln(x + sqrt(x^2 - 1)).
*/
NUMBER *
qacosh(NUMBER *q, NUMBER *epsilon)
{
NUMBER *tmp1, *tmp2, *epsilon1;
long n;
if (qiszero(epsilon)) {
math_error("Zero epsilon value for acosh");
not_reached();
}
if (qisone(q))
return qlink(&_qzero_);
if (zrel(q->num, q->den) < 0)
return NULL;
n = qilog2(epsilon);
epsilon1 = qbitvalue(n - 3);
tmp1 = qsquare(q);
tmp2 = qdec(tmp1);
qfree(tmp1);
tmp1 = qsqrt(tmp2, epsilon1, conf->triground);
qfree(tmp2);
tmp2 = qqadd(tmp1, q);
qfree(tmp1);
tmp1 = qln(tmp2, epsilon1);
qfree(tmp2);
qfree(epsilon1);
tmp2 = qmappr(tmp1, epsilon, conf->triground);
qfree(tmp1);
return tmp2;
}
/*
* Compute the hyperbolic arcsine within the specified accuracy.
* This is calculated using the formula:
* asinh(x) = ln(x + sqrt(x^2 + 1)).
*/
NUMBER *
qasinh(NUMBER *q, NUMBER *epsilon)
{
NUMBER *tmp1, *tmp2, *epsilon1;
long n;
bool neg;
if (qiszero(epsilon)) {
math_error("Zero epsilon value for asinh");
not_reached();
}
if (qiszero(q))
return qlink(&_qzero_);
neg = qisneg(q);
q = qqabs(q);
n = qilog2(epsilon);
epsilon1 = qbitvalue(n - 3);
tmp1 = qsquare(q);
tmp2 = qinc(tmp1);
qfree(tmp1);
tmp1 = qsqrt(tmp2, epsilon1, conf->triground);
qfree(tmp2);
tmp2 = qqadd(tmp1, q);
qfree(tmp1);
tmp1 = qln(tmp2, epsilon1);
qfree(tmp2);
qfree(q);
qfree(epsilon1);
tmp2 = qmappr(tmp1, epsilon, conf->triground);
if (neg) {
tmp1 = qneg(tmp2);
qfree(tmp2);
return tmp1;
}
return tmp2;
}
/*
* Compute the hyperbolic arctangent within the specified accuracy.
* This is calculated using the formula:
* atanh(x) = ln((1 + x) / (1 - x)) / 2.
*/
NUMBER *
qatanh(NUMBER *q, NUMBER *epsilon)
{
NUMBER *tmp1, *tmp2, *tmp3, *epsilon1;
ZVALUE z;
if (qiszero(epsilon)) {
math_error("Zero epsilon value for atanh");
not_reached();
}
if (qiszero(q))
return qlink(&_qzero_);
z = q->num;
z.sign = 0;
if (zrel(z, q->den) >= 0)
return NULL;
tmp1 = qinc(q);
tmp2 = qsub(&_qone_, q);
tmp3 = qqdiv(tmp1, tmp2);
qfree(tmp1);
qfree(tmp2);
epsilon1 = qscale(epsilon, 1L);
tmp1 = qln(tmp3, epsilon1);
qfree(tmp3);
tmp2 = qscale(tmp1, -1L);
qfree(tmp1);
qfree(epsilon1);
return tmp2;
}
/*
* Inverse hyperbolic secant function
*/
NUMBER *
qasech(NUMBER *q, NUMBER *epsilon)
{
NUMBER *tmp, *res;
tmp = qinv(q);
res = qacosh(tmp, epsilon);
qfree(tmp);
return res;
}
/*
* Inverse hyperbolic cosecant function
*/
NUMBER *
qacsch(NUMBER *q, NUMBER *epsilon)
{
NUMBER *tmp, *res;
tmp = qinv(q);
res = qasinh(tmp, epsilon);
qfree(tmp);
return res;
}
/*
* Inverse hyperbolic cotangent function
*/
NUMBER *
qacoth(NUMBER *q, NUMBER *epsilon)
{
NUMBER *tmp, *res;
tmp = qinv(q);
res = qatanh(tmp, epsilon);
qfree(tmp);
return res;
}
/*
* qversin - versed trigonometric sine
*
* This uses the formula:
*
* versin(x) = 1 - cos(x)
*
* given:
* q real value to pass to the trig function
* epsilon error tolerance / precision for trig calculation
*
* returns:
* real value result of trig function on q with error epsilon
*/
NUMBER *
qversin(NUMBER *q, NUMBER *epsilon)
{
NUMBER *res; /* inverse trig value result */
NUMBER *qtmp; /* argument to inverse trig function */
/*
* firewall
*/
if (q == NULL) {
math_error("q is NULL for %s", __func__);
not_reached();
}
if (check_epsilon(epsilon) == false) {
math_error("Invalid epsilon arg for %s", __func__);
not_reached();
}
/*
* calculate trig function value
*/
qtmp = qcos(q, epsilon);
res = qsub(&_qone_, qtmp);
qfree(qtmp);
/*
* return trigonometric result
*/
return res;
}
/*
* qaversin_or_NULL - return NULL or non-complex inverse versed trigonometric sine
*
* This uses the formula:
*
* aversin(x) = acos(1 - x)
*
* given:
* q real value to pass to the trig function
* epsilon error tolerance / precision for trig calculation
*
* returns:
* != NULL ==> real value result of trig function on q with error epsilon,
* NULL ==> trig function value cannot be expressed as a NUMBER
*
* NOTE: If this function returns NULL, consider calling the equivalent
* COMPLEX function from comfunc.c. See the help file for the
* related builtin for details.
*/
NUMBER *
qaversin_or_NULL(NUMBER *q, NUMBER *epsilon)
{
NUMBER *res; /* inverse trig value result */
NUMBER *qtmp; /* argument to inverse trig function */
/*
* firewall
*/
if (q == NULL) {
math_error("q is NULL for %s", __func__);
not_reached();
}
if (check_epsilon(epsilon) == false) {
math_error("Invalid epsilon arg for %s", __func__);
not_reached();
}
/*
* calculate inverse trig function value
*/
qtmp = qsub(&_qone_, q);
res = qacos(qtmp, epsilon);
qfree(qtmp);
if (res == NULL) {
return NULL;
}
/*
* return inverse trigonometric result
*/
return res;
}
/*
* qaversin - non-complex inverse versed trigonometric sine
*
* This uses the formula:
*
* aversin(x) = acos(1 - x)
*
* given:
* q real value to pass to the trig function
* epsilon error tolerance / precision for trig calculation
*
* returns:
* real value result of trig function on q with error epsilon
*/
NUMBER *
qaversin(NUMBER *q, NUMBER *epsilon)
{
NUMBER *res; /* inverse trig value result */
/*
* firewall
*/
if (q == NULL) {
math_error("q is NULL for %s", __func__);
not_reached();
}
if (check_epsilon(epsilon) == false) {
math_error("Invalid epsilon arg for %s", __func__);
not_reached();
}
/*
* calculate inverse trig function value
*/
res = qaversin_or_NULL(q, epsilon);
if (res == NULL) {
math_error("cannot compute inverse cosine for aversin");
not_reached();
}
/*
* return inverse trigonometric result
*/
return res;
}
/*
* qcoversin - coversed trigonometric sine
*
* This uses the formula:
*
* coversin((x) = 1 - sin(x)
*
* given:
* q real value to pass to the trig function
* epsilon error tolerance / precision for trig calculation
*
* returns:
* real value result of trig function on q with error epsilon
*/
NUMBER *
qcoversin(NUMBER *q, NUMBER *epsilon)
{
NUMBER *res; /* inverse trig value result */
NUMBER *qtmp; /* argument to inverse trig function */
/*
* firewall
*/
if (q == NULL) {
math_error("q is NULL for %s", __func__);
not_reached();
}
if (check_epsilon(epsilon) == false) {
math_error("Invalid epsilon arg for %s", __func__);
not_reached();
}
/*
* calculate trig function value
*/
qtmp = qsin(q, epsilon);
res = qsub(&_qone_, qtmp);
qfree(qtmp);
/*
* return trigonometric result
*/
return res;
}
/*
* qacoversin_or_NULL - return NULL or non-complex inverse coversed trigonometric sine
*
* This uses the formula:
*
* acoversin(x) = asin(1 - x)
*
* given:
* q real value to pass to the trig function
* epsilon error tolerance / precision for trig calculation
*
* returns:
* real value result of trig function on q with error epsilon
*
* returns:
* != NULL ==> real value result of trig function on q with error epsilon,
* NULL ==> trig function value cannot be expressed as a NUMBER
*
* NOTE: If this function returns NULL, consider calling the equivalent
* COMPLEX function from comfunc.c. See the help file for the
* related builtin for details.
*/
NUMBER *
qacoversin_or_NULL(NUMBER *q, NUMBER *epsilon)
{
NUMBER *res; /* inverse trig value result */
NUMBER *qtmp; /* argument to inverse trig function */
/*
* firewall
*/
if (q == NULL) {
math_error("q is NULL for %s", __func__);
not_reached();
}
if (check_epsilon(epsilon) == false) {
math_error("Invalid epsilon arg for %s", __func__);
not_reached();
}
/*
* calculate inverse trig function value
*/
qtmp = qsub(&_qone_, q);
res = qasin(qtmp, epsilon);
qfree(qtmp);
if (res == NULL) {
return NULL;
}
/*
* return inverse trigonometric result
*/
return res;
}
/*
* qacoversin - non-complex inverse coversed trigonometric sine
*
* This uses the formula:
*
* acoversin(x) = asin(1 - x)
*
* given:
* q real value to pass to the trig function
* epsilon error tolerance / precision for trig calculation
*
* returns:
* real value result of trig function on q with error epsilon
*/
NUMBER *
qacoversin(NUMBER *q, NUMBER *epsilon)
{
NUMBER *res; /* inverse trig value result */
/*
* firewall
*/
if (q == NULL) {
math_error("q is NULL for %s", __func__);
not_reached();
}
if (check_epsilon(epsilon) == false) {
math_error("Invalid epsilon arg for %s", __func__);
not_reached();
}
/*
* calculate inverse trig function value
*/
res = qacoversin_or_NULL(q, epsilon);
if (res == NULL) {
math_error("cannot compute inverse sine for acoversin");
not_reached();
}
/*
* return inverse trigonometric result
*/
return res;
}
/*
* qvercos - versed trigonometric cosine
*
* This uses the formula:
*
* vercos(x) = 1 + cos(x)
*
* given:
* q real value to pass to the trig function
* epsilon error tolerance / precision for trig calculation
*
* returns:
* real value result of trig function on q with error epsilon
*/
NUMBER *
qvercos(NUMBER *q, NUMBER *epsilon)
{
NUMBER *res; /* inverse trig value result */
NUMBER *qtmp; /* argument to inverse trig function */
/*
* firewall
*/
if (q == NULL) {
math_error("q is NULL for %s", __func__);
not_reached();
}
if (check_epsilon(epsilon) == false) {
math_error("Invalid epsilon arg for %s", __func__);
not_reached();
}
/*
* calculate trig function value
*/
qtmp = qcos(q, epsilon);
res = qqadd(&_qone_, qtmp);
qfree(qtmp);
/*
* return trigonometric result
*/
return res;
}
/*
* qavercos_or_NULL - return NULL or non-complex inverse versed trigonometric cosine
*
* This uses the formula:
*
* avercos(x) = acos(x - 1)
*
* given:
* q real value to pass to the trig function
* epsilon error tolerance / precision for trig calculation
*
* returns:
* != NULL ==> real value result of trig function on q with error epsilon,
* NULL ==> trig function value cannot be expressed as a NUMBER
*
* NOTE: If this function returns NULL, consider calling the equivalent
* COMPLEX function from comfunc.c. See the help file for the
* related builtin for details.
*/
NUMBER *
qavercos_or_NULL(NUMBER *q, NUMBER *epsilon)
{
NUMBER *res; /* inverse trig value result */
NUMBER *qtmp; /* argument to inverse trig function */
/*
* firewall
*/
if (q == NULL) {
math_error("q is NULL for %s", __func__);
not_reached();
}
if (check_epsilon(epsilon) == false) {
math_error("Invalid epsilon arg for %s", __func__);
not_reached();
}
/*
* calculate inverse trig function value
*/
qtmp = qsub(q, &_qone_);
res = qacos(qtmp, epsilon);
qfree(qtmp);
if (res == NULL) {
return NULL;
}
/*
* return inverse trigonometric result
*/
return res;
}
/*
* qavercos - non-complex inverse versed trigonometric cosine
*
* This uses the formula:
*
* avercos(x) = acos(x - 1)
*
* given:
* q real value to pass to the trig function
* epsilon error tolerance / precision for trig calculation
*
* returns:
* real value result of trig function on q with error epsilon
*/
NUMBER *
qavercos(NUMBER *q, NUMBER *epsilon)
{
NUMBER *res; /* inverse trig value result */
/*
* firewall
*/
if (q == NULL) {
math_error("q is NULL for %s", __func__);
not_reached();
}
if (check_epsilon(epsilon) == false) {
math_error("Invalid epsilon arg for %s", __func__);
not_reached();
}
/*
* calculate inverse trig function value
*/
res = qavercos_or_NULL(q, epsilon);
if (res == NULL) {
math_error("cannot compute inverse cosine for avercos");
not_reached();
}
/*
* return inverse trigonometric result
*/
return res;
}
/*
* qcovercos - coversed trigonometric cosine
*
* This uses the formula:
*
* covercos((x) = 1 + sin(x)
*
* given:
* q real value to pass to the trig function
* epsilon error tolerance / precision for trig calculation
*
* returns:
* real value result of trig function on q with error epsilon
*/
NUMBER *
qcovercos(NUMBER *q, NUMBER *epsilon)
{
NUMBER *res; /* inverse trig value result */
NUMBER *qtmp; /* argument to inverse trig function */
/*
* firewall
*/
if (q == NULL) {
math_error("q is NULL for %s", __func__);
not_reached();
}
if (check_epsilon(epsilon) == false) {
math_error("Invalid epsilon arg for %s", __func__);
not_reached();
}
/*
* calculate trig function value
*/
qtmp = qsin(q, epsilon);
res = qqadd(&_qone_, qtmp);
qfree(qtmp);
/*
* return trigonometric result
*/
return res;
}
/*
* qacovercos_or_NULL - return NULL or non-complex inverse coversed trigonometric cosine
*
* This uses the formula:
*
* acovercos(x) = asin(x - 1)
*
* given:
* q real value to pass to the trig function
* epsilon error tolerance / precision for trig calculation
*
* returns:
* real value result of trig function on q with error epsilon
*
* returns:
* != NULL ==> real value result of trig function on q with error epsilon,
* NULL ==> trig function value cannot be expressed as a NUMBER
*
* NOTE: If this function returns NULL, consider calling the equivalent
* COMPLEX function from comfunc.c. See the help file for the
* related builtin for details.
*/
NUMBER *
qacovercos_or_NULL(NUMBER *q, NUMBER *epsilon)
{
NUMBER *res; /* inverse trig value result */
NUMBER *qtmp; /* argument to inverse trig function */
/*
* firewall
*/
if (q == NULL) {
math_error("q is NULL for %s", __func__);
not_reached();
}
if (check_epsilon(epsilon) == false) {
math_error("Invalid epsilon arg for %s", __func__);
not_reached();
}
/*
* calculate inverse trig function value
*/
qtmp = qsub(&_qone_, q);
res = qasin(qtmp, epsilon);
qfree(qtmp);
if (res == NULL) {
return NULL;
}
/*
* return inverse trigonometric result
*/
return res;
}
/*
* qacovercos - non-complex inverse coversed trigonometric cosine
*
* This uses the formula:
*
* acovercos(x) = asin(x - 1)
*
* given:
* q real value to pass to the trig function
* epsilon error tolerance / precision for trig calculation
*
* returns:
* real value result of trig function on q with error epsilon
*/
NUMBER *
qacovercos(NUMBER *q, NUMBER *epsilon)
{
NUMBER *res; /* inverse trig value result */
/*
* firewall
*/
if (q == NULL) {
math_error("q is NULL for %s", __func__);
not_reached();
}
if (check_epsilon(epsilon) == false) {
math_error("Invalid epsilon arg for %s", __func__);
not_reached();
}
/*
* calculate inverse trig function value
*/
res = qacovercos_or_NULL(q, epsilon);
if (res == NULL) {
math_error("cannot compute inverse sine for acovercos");
not_reached();
}
/*
* return inverse trigonometric result
*/
return res;
}
/*
* qhaversin - half versed trigonometric sine
*
* This uses the formula:
*
* haversin(x) = versin(x) / 2
*
* given:
* q real value to pass to the trig function
* epsilon error tolerance / precision for trig calculation
*
* returns:
* real value result of trig function on q with error epsilon
*/
NUMBER *
qhaversin(NUMBER *q, NUMBER *epsilon)
{
NUMBER *res; /* inverse trig value result */
NUMBER *qtmp; /* argument to inverse trig function */
/*
* firewall
*/
if (q == NULL) {
math_error("q is NULL for %s", __func__);
not_reached();
}
if (check_epsilon(epsilon) == false) {
math_error("Invalid epsilon arg for %s", __func__);
not_reached();
}
/*
* calculate trig function value
*/
qtmp = qversin(q, epsilon);
res = qdivi(qtmp, 2);
qfree(qtmp);
/*
* return trigonometric result
*/
return res;
}
/*
* qahaversin_or_NULL - return NULL or non-complex inverse half versed trigonometric sine
*
* This uses the formula:
*
* ahaversin(x) = acos(1 - 2*x)
*
* given:
* q real value to pass to the trig function
* epsilon error tolerance / precision for trig calculation
*
* returns:
* != NULL ==> real value result of trig function on q with error epsilon,
* NULL ==> trig function value cannot be expressed as a NUMBER
*
* NOTE: If this function returns NULL, consider calling the equivalent
* COMPLEX function from comfunc.c. See the help file for the
* related builtin for details.
*/
NUMBER *
qahaversin_or_NULL(NUMBER *q, NUMBER *epsilon)
{
NUMBER *res; /* inverse trig value result */
NUMBER *qtmp; /* argument to inverse trig function */
NUMBER *x2; /* twice x */
/*
* firewall
*/
if (q == NULL) {
math_error("q is NULL for %s", __func__);
not_reached();
}
if (check_epsilon(epsilon) == false) {
math_error("Invalid epsilon arg for %s", __func__);
not_reached();
}
/*
* calculate inverse trig function value
*/
x2 = qmuli(q, 2);
qtmp = qsub(&_qone_, x2);
qfree(x2);
res = qacos(qtmp, epsilon);
qfree(qtmp);
if (res == NULL) {
return NULL;
}
/*
* return inverse trigonometric result
*/
return res;
}
/*
* qahaversin - non-complex inverse half versed trigonometric sine
*
* This uses the formula:
*
* ahaversin(x) = acos(1 - 2*x)
*
* given:
* q real value to pass to the trig function
* epsilon error tolerance / precision for trig calculation
*
* returns:
* real value result of trig function on q with error epsilon
*/
NUMBER *
qahaversin(NUMBER *q, NUMBER *epsilon)
{
NUMBER *res; /* inverse trig value result */
/*
* firewall
*/
if (q == NULL) {
math_error("q is NULL for %s", __func__);
not_reached();
}
if (check_epsilon(epsilon) == false) {
math_error("Invalid epsilon arg for %s", __func__);
not_reached();
}
/*
* calculate inverse trig function value
*/
res = qahaversin_or_NULL(q, epsilon);
if (res == NULL) {
math_error("cannot compute inverse cosine for ahaversin");
not_reached();
}
/*
* return inverse trigonometric result
*/
return res;
}
/*
* qhacoversin - coversed trigonometric sine
*
* This uses the formula:
*
* hacoversin((x) = coversin(x) / 2
*
* given:
* q real value to pass to the trig function
* epsilon error tolerance / precision for trig calculation
*
* returns:
* real value result of trig function on q with error epsilon
*/
NUMBER *
qhacoversin(NUMBER *q, NUMBER *epsilon)
{
NUMBER *res; /* inverse trig value result */
NUMBER *qtmp; /* argument to inverse trig function */
/*
* firewall
*/
if (q == NULL) {
math_error("q is NULL for %s", __func__);
not_reached();
}
if (check_epsilon(epsilon) == false) {
math_error("Invalid epsilon arg for %s", __func__);
not_reached();
}
/*
* calculate trig function value
*/
qtmp = qcoversin(q, epsilon);
res = qdivi(qtmp, 2);
qfree(qtmp);
/*
* return trigonometric result
*/
return res;
}
/*
* qahacoversin_or_NULL - return NULL or non-complex inverse coversed trigonometric sine
*
* This uses the formula:
*
* ahacoversin(x) = asin(1 - 2*x)
*
* given:
* q real value to pass to the trig function
* epsilon error tolerance / precision for trig calculation
*
* returns:
* real value result of trig function on q with error epsilon
*
* returns:
* != NULL ==> real value result of trig function on q with error epsilon,
* NULL ==> trig function value cannot be expressed as a NUMBER
*
* NOTE: If this function returns NULL, consider calling the equivalent
* COMPLEX function from comfunc.c. See the help file for the
* related builtin for details.
*/
NUMBER *
qahacoversin_or_NULL(NUMBER *q, NUMBER *epsilon)
{
NUMBER *res; /* inverse trig value result */
NUMBER *qtmp; /* argument to inverse trig function */
/*
* firewall
*/
if (q == NULL) {
math_error("q is NULL for %s", __func__);
not_reached();
}
if (check_epsilon(epsilon) == false) {
math_error("Invalid epsilon arg for %s", __func__);
not_reached();
}
/*
* calculate inverse trig function value
*/
qtmp = qsub(&_qone_, q);
res = qasin(qtmp, epsilon);
qfree(qtmp);
if (res == NULL) {
return NULL;
}
/*
* return inverse trigonometric result
*/
return res;
}
/*
* qahacoversin - non-complex inverse coversed trigonometric sine
*
* This uses the formula:
*
* ahacoversin(x) = asin(1 - 2*x)
*
* given:
* q real value to pass to the trig function
* epsilon error tolerance / precision for trig calculation
*
* returns:
* real value result of trig function on q with error epsilon
*/
NUMBER *
qahacoversin(NUMBER *q, NUMBER *epsilon)
{
NUMBER *res; /* inverse trig value result */
/*
* firewall
*/
if (q == NULL) {
math_error("q is NULL for %s", __func__);
not_reached();
}
if (check_epsilon(epsilon) == false) {
math_error("Invalid epsilon arg for %s", __func__);
not_reached();
}
/*
* calculate inverse trig function value
*/
res = qahacoversin_or_NULL(q, epsilon);
if (res == NULL) {
math_error("cannot compute inverse sine for ahacoversin");
not_reached();
}
/*
* return inverse trigonometric result
*/
return res;
}
/*
* qhavercos - half versed trigonometric cosine
*
* This uses the formula:
*
* havercos(x) = vercos(x) / 2
*
* given:
* q real value to pass to the trig function
* epsilon error tolerance / precision for trig calculation
*
* returns:
* real value result of trig function on q with error epsilon
*/
NUMBER *
qhavercos(NUMBER *q, NUMBER *epsilon)
{
NUMBER *res; /* inverse trig value result */
NUMBER *qtmp; /* argument to inverse trig function */
/*
* firewall
*/
if (q == NULL) {
math_error("q is NULL for %s", __func__);
not_reached();
}
if (check_epsilon(epsilon) == false) {
math_error("Invalid epsilon arg for %s", __func__);
not_reached();
}
/*
* calculate trig function value
*/
qtmp = qvercos(q, epsilon);
res = qdivi(qtmp, 2);
qfree(qtmp);
/*
* return trigonometric result
*/
return res;
}
/*
* qahavercos_or_NULL - return NULL or non-complex inverse half versed trigonometric cosine
*
* This uses the formula:
*
* ahavercos(x) = acos(1 - 2*x)
*
* given:
* q real value to pass to the trig function
* epsilon error tolerance / precision for trig calculation
*
* returns:
* != NULL ==> real value result of trig function on q with error epsilon,
* NULL ==> trig function value cannot be expressed as a NUMBER
*
* NOTE: If this function returns NULL, consider calling the equivalent
* COMPLEX function from comfunc.c. See the help file for the
* related builtin for details.
*/
NUMBER *
qahavercos_or_NULL(NUMBER *q, NUMBER *epsilon)
{
NUMBER *res; /* inverse trig value result */
NUMBER *qtmp; /* argument to inverse trig function */
NUMBER *x2; /* twice x */
/*
* firewall
*/
if (q == NULL) {
math_error("q is NULL for %s", __func__);
not_reached();
}
if (check_epsilon(epsilon) == false) {
math_error("Invalid epsilon arg for %s", __func__);
not_reached();
}
/*
* calculate inverse trig function value
*/
x2 = qmuli(q, 2);
qtmp = qsub(&_qone_, x2);
qfree(x2);
res = qacos(qtmp, epsilon);
qfree(qtmp);
if (res == NULL) {
return NULL;
}
/*
* return inverse trigonometric result
*/
return res;
}
/*
* qahavercos - non-complex inverse half versed trigonometric cosine
*
* This uses the formula:
*
* ahavercos(x) = acos(1 - 2*x)
*
* given:
* q real value to pass to the trig function
* epsilon error tolerance / precision for trig calculation
*
* returns:
* real value result of trig function on q with error epsilon
*/
NUMBER *
qahavercos(NUMBER *q, NUMBER *epsilon)
{
NUMBER *res; /* inverse trig value result */
/*
* firewall
*/
if (q == NULL) {
math_error("q is NULL for %s", __func__);
not_reached();
}
if (check_epsilon(epsilon) == false) {
math_error("Invalid epsilon arg for %s", __func__);
not_reached();
}
/*
* calculate inverse trig function value
*/
res = qahaversin_or_NULL(q, epsilon);
if (res == NULL) {
math_error("cannot compute inverse cocosine for ahavercos");
not_reached();
}
/*
* return inverse trigonometric result
*/
return res;
}
/*
* qhacovercos - half coversed trigonometric cosine
*
* This uses the formula:
*
* hacovercos((x) = covercos(x) / 2
*
* given:
* q real value to pass to the trig function
* epsilon error tolerance / precision for trig calculation
*
* returns:
* real value result of trig function on q with error epsilon
*/
NUMBER *
qhacovercos(NUMBER *q, NUMBER *epsilon)
{
NUMBER *res; /* inverse trig value result */
NUMBER *qtmp; /* argument to inverse trig function */
/*
* firewall
*/
if (q == NULL) {
math_error("q is NULL for %s", __func__);
not_reached();
}
if (check_epsilon(epsilon) == false) {
math_error("Invalid epsilon arg for %s", __func__);
not_reached();
}
/*
* calculate trig function value
*/
qtmp = qcovercos(q, epsilon);
res = qdivi(qtmp, 2);
qfree(qtmp);
/*
* return trigonometric result
*/
return res;
}
/*
* qahacovercos_or_NULL - return NULL or non-complex inverse half coversed trigonometric cosine
*
* This uses the formula:
*
* ahacovercos(x) = acos(2*x - 1)
*
* given:
* q real value to pass to the trig function
* epsilon error tolerance / precision for trig calculation
*
* returns:
* real value result of trig function on q with error epsilon
*
* returns:
* != NULL ==> real value result of trig function on q with error epsilon,
* NULL ==> trig function value cannot be expressed as a NUMBER
*
* NOTE: If this function returns NULL, consider calling the equivalent
* COMPLEX function from comfunc.c. See the help file for the
* related builtin for details.
*/
NUMBER *
qahacovercos_or_NULL(NUMBER *q, NUMBER *epsilon)
{
NUMBER *res; /* inverse trig value result */
NUMBER *qtmp; /* argument to inverse trig function */
/*
* firewall
*/
if (q == NULL) {
math_error("q is NULL for %s", __func__);
not_reached();
}
if (check_epsilon(epsilon) == false) {
math_error("Invalid epsilon arg for %s", __func__);
not_reached();
}
/*
* calculate inverse trig function value
*/
qtmp = qsub(&_qone_, q);
res = qacos(qtmp, epsilon);
qfree(qtmp);
if (res == NULL) {
return NULL;
}
/*
* return inverse trigonometric result
*/
return res;
}
/*
* qahacovercos - non-complex inverse half coversed trigonometric cosine
*
* This uses the formula:
*
* ahacovercos(x) = acos(2*x - 1)
*
* given:
* q real value to pass to the trig function
* epsilon error tolerance / precision for trig calculation
*
* returns:
* real value result of trig function on q with error epsilon
*/
NUMBER *
qahacovercos(NUMBER *q, NUMBER *epsilon)
{
NUMBER *res; /* inverse trig value result */
/*
* firewall
*/
if (q == NULL) {
math_error("q is NULL for %s", __func__);
not_reached();
}
if (check_epsilon(epsilon) == false) {
math_error("Invalid epsilon arg for %s", __func__);
not_reached();
}
/*
* calculate inverse trig function value
*/
res = qahacoversin_or_NULL(q, epsilon);
if (res == NULL) {
math_error("cannot compute inverse cosine for ahacovercos");
not_reached();
}
/*
* return inverse trigonometric result
*/
return res;
}
/*
* qexsec - exterior trigonometric secant
*
* This uses the formula:
*
* exsec(x) = sec(x) - 1
*
* given:
* q real value to pass to the trig function
* epsilon error tolerance / precision for trig calculation
*
* returns:
* real value result of trig function on q with error epsilon
*/
NUMBER *
qexsec(NUMBER *q, NUMBER *epsilon)
{
NUMBER *res; /* inverse trig value result */
NUMBER *qtmp; /* argument to inverse trig function */
/*
* firewall
*/
if (q == NULL) {
math_error("q is NULL for %s", __func__);
not_reached();
}
if (check_epsilon(epsilon) == false) {
math_error("Invalid epsilon arg for %s", __func__);
not_reached();
}
/*
* calculate trig function value
*/
qtmp = qsec(q, epsilon);
res = qsub(qtmp, &_qone_);
qfree(qtmp);
/*
* return trigonometric result
*/
return res;
}
/*
* qaexsec_or_NULL - return NULL or non-complex inverse versed trigonometric sine
*
* This uses the formula:
*
* aexsec(x) = asec(x + 1)
*
* given:
* q real value to pass to the trig function
* epsilon error tolerance / precision for trig calculation
*
* returns:
* != NULL ==> real value result of trig function on q with error epsilon,
* NULL ==> trig function value cannot be expressed as a NUMBER
*
* NOTE: If this function returns NULL, consider calling the equivalent
* COMPLEX function from comfunc.c. See the help file for the
* related builtin for details.
*/
NUMBER *
qaexsec_or_NULL(NUMBER *q, NUMBER *epsilon)
{
NUMBER *res; /* inverse trig value result */
NUMBER *qtmp; /* argument to inverse trig function */
/*
* firewall
*/
if (q == NULL) {
math_error("q is NULL for %s", __func__);
not_reached();
}
if (check_epsilon(epsilon) == false) {
math_error("Invalid epsilon arg for %s", __func__);
not_reached();
}
/*
* calculate inverse trig function value
*/
qtmp = qaddi(q, 1);
res = qasec(qtmp, epsilon);
qfree(qtmp);
if (res == NULL) {
return NULL;
}
/*
* return inverse trigonometric result
*/
return res;
}
/*
* qaexsec - non-complex inverse versed trigonometric sine
*
* This uses the formula:
*
* aexsec(x) = asec(x + 1)
*
* given:
* q real value to pass to the trig function
* epsilon error tolerance / precision for trig calculation
*
* returns:
* real value result of trig function on q with error epsilon
*/
NUMBER *
qaexsec(NUMBER *q, NUMBER *epsilon)
{
NUMBER *res; /* inverse trig value result */
/*
* firewall
*/
if (q == NULL) {
math_error("q is NULL for %s", __func__);
not_reached();
}
if (check_epsilon(epsilon) == false) {
math_error("Invalid epsilon arg for %s", __func__);
not_reached();
}
/*
* calculate inverse trig function value
*/
res = qaexsec_or_NULL(q, epsilon);
if (res == NULL) {
math_error("cannot compute inverse cosine for aexsec");
not_reached();
}
/*
* return inverse trigonometric result
*/
return res;
}
/*
* qexcsc - coversed trigonometric sine
*
* This uses the formula:
*
* excsc(x) = csc(x) - 1
*
* given:
* q real value to pass to the trig function
* epsilon error tolerance / precision for trig calculation
*
* returns:
* real value result of trig function on q with error epsilon
*/
NUMBER *
qexcsc(NUMBER *q, NUMBER *epsilon)
{
NUMBER *res; /* inverse trig value result */
NUMBER *qtmp; /* argument to inverse trig function */
/*
* firewall
*/
if (q == NULL) {
math_error("q is NULL for %s", __func__);
not_reached();
}
if (check_epsilon(epsilon) == false) {
math_error("Invalid epsilon arg for %s", __func__);
not_reached();
}
/*
* calculate trig function value
*/
qtmp = qcsc(q, epsilon);
res = qsub(qtmp, &_qone_);
qfree(qtmp);
/*
* return trigonometric result
*/
return res;
}
/*
* qaexcsc_or_NULL - return NULL or non-complex inverse coversed trigonometric sine
*
* This uses the formula:
*
* aexcsc(x) = acsc(x + 1)
*
* given:
* q real value to pass to the trig function
* epsilon error tolerance / precision for trig calculation
*
* returns:
* real value result of trig function on q with error epsilon
*
* returns:
* != NULL ==> real value result of trig function on q with error epsilon,
* NULL ==> trig function value cannot be expressed as a NUMBER
*
* NOTE: If this function returns NULL, consider calling the equivalent
* COMPLEX function from comfunc.c. See the help file for the
* related builtin for details.
*/
NUMBER *
qaexcsc_or_NULL(NUMBER *q, NUMBER *epsilon)
{
NUMBER *res; /* inverse trig value result */
NUMBER *qtmp; /* argument to inverse trig function */
/*
* firewall
*/
if (q == NULL) {
math_error("q is NULL for %s", __func__);
not_reached();
}
if (check_epsilon(epsilon) == false) {
math_error("Invalid epsilon arg for %s", __func__);
not_reached();
}
/*
* calculate inverse trig function value
*/
qtmp = qaddi(q, 1);
res = qacsc(qtmp, epsilon);
qfree(qtmp);
if (res == NULL) {
return NULL;
}
/*
* return inverse trigonometric result
*/
return res;
}
/*
* qaexcsc - non-complex inverse coversed trigonometric sine
*
* This uses the formula:
*
* aexcsc(x) = acsc(x + 1)
*
* given:
* q real value to pass to the trig function
* epsilon error tolerance / precision for trig calculation
*
* returns:
* real value result of trig function on q with error epsilon
*/
NUMBER *
qaexcsc(NUMBER *q, NUMBER *epsilon)
{
NUMBER *res; /* inverse trig value result */
/*
* firewall
*/
if (q == NULL) {
math_error("q is NULL for %s", __func__);
not_reached();
}
if (check_epsilon(epsilon) == false) {
math_error("Invalid epsilon arg for %s", __func__);
not_reached();
}
/*
* calculate inverse trig function value
*/
res = qaexcsc_or_NULL(q, epsilon);
if (res == NULL) {
math_error("cannot compute inverse sine for aexcsc");
not_reached();
}
/*
* return inverse trigonometric result
*/
return res;
}
/*
* qcrd - trigonometric chord of a unit circle
*
* This uses the formula:
*
* crd(x) = 2 * sin(x / 2)
*
* given:
* q real value to pass to the trig function
* epsilon error tolerance / precision for trig calculation
*
* returns:
* real value result of trig function on q with error epsilon
*/
NUMBER *
qcrd(NUMBER *q, NUMBER *epsilon)
{
NUMBER *res; /* inverse trig value result */
NUMBER *qdiv2; /* q/2 */
NUMBER *qtmp; /* argument to inverse trig function */
/*
* firewall
*/
if (q == NULL) {
math_error("q is NULL for %s", __func__);
not_reached();
}
if (check_epsilon(epsilon) == false) {
math_error("Invalid epsilon arg for %s", __func__);
not_reached();
}
/*
* calculate trig function value
*/
qdiv2 = qdivi(q, 2);
qtmp = qsin(qdiv2, epsilon);
qfree(qdiv2);
res = qmuli(qtmp, 2);
qfree(qtmp);
/*
* return trigonometric result
*/
return res;
}
/*
* qacrd_or_NULL - return NULL or non-complex inverse trigonometric chord of a unit circle
*
* This uses the formula:
*
* acrd(x) = 2 * asin(x / 2)
*
* given:
* q real value to pass to the trig function
* epsilon error tolerance / precision for trig calculation
*
* returns:
* real value result of trig function on q with error epsilon
*
* returns:
* != NULL ==> real value result of trig function on q with error epsilon,
* NULL ==> trig function value cannot be expressed as a NUMBER
*
* NOTE: If this function returns NULL, consider calling the equivalent
* COMPLEX function from comfunc.c. See the help file for the
* related builtin for details.
*/
NUMBER *
qacrd_or_NULL(NUMBER *q, NUMBER *epsilon)
{
NUMBER *res; /* inverse trig value result */
NUMBER *qdiv2; /* q/2 */
NUMBER *qtmp; /* argument to inverse trig function */
/*
* firewall
*/
if (q == NULL) {
math_error("q is NULL for %s", __func__);
not_reached();
}
if (check_epsilon(epsilon) == false) {
math_error("Invalid epsilon arg for %s", __func__);
not_reached();
}
/*
* calculate inverse trig function value
*/
qdiv2 = qdivi(q, 2);
qtmp = qasin(qdiv2, epsilon);
qfree(qdiv2);
if (qtmp == NULL) {
return NULL;
}
res = qmuli(qtmp, 2);
qfree(qtmp);
/*
* return inverse trigonometric result
*/
return res;
}
/*
* qacrd - non-complex inverse trigonometric chord of a unit circle
*
* This uses the formula:
*
* acrd(x) = 2 * asin(x / 2)
*
* given:
* q real value to pass to the trig function
* epsilon error tolerance / precision for trig calculation
*
* returns:
* real value result of trig function on q with error epsilon
*/
NUMBER *
qacrd(NUMBER *q, NUMBER *epsilon)
{
NUMBER *res; /* inverse trig value result */
/*
* firewall
*/
if (q == NULL) {
math_error("q is NULL for %s", __func__);
not_reached();
}
if (check_epsilon(epsilon) == false) {
math_error("Invalid epsilon arg for %s", __func__);
not_reached();
}
/*
* calculate inverse trig function value
*/
res = qacrd_or_NULL(q, epsilon);
if (res == NULL) {
math_error("cannot compute inverse sine for acrd");
not_reached();
}
/*
* return inverse trigonometric result
*/
return res;
}
/*
* qcas - trigonometric chord of a unit circle
*
* This uses the formula:
*
* cas(x) = cos(x) + sin(x)
*
* given:
* q real value to pass to the trig function
* epsilon error tolerance / precision for trig calculation
*
* returns:
* real value result of trig function on q with error epsilon
*/
NUMBER *
qcas(NUMBER *q, NUMBER *epsilon)
{
NUMBER *sin; /* sin(x) */
NUMBER *tsin; /* sin(x) rounded to nearest epsilon multiple */
NUMBER *cos; /* cos(x) */
NUMBER *tcos; /* cos(x) rounded to nearest epsilon multiple */
NUMBER *res;
long n;
/*
* firewall
*/
if (qiszero(epsilon)) {
math_error("Zero epsilon value for cosine");
not_reached();
}
/*
* case 0: quick return 1
*/
if (qiszero(q)) {
return qlink(&_qone_);
}
/*
* case epsilon > 1: quick return 0
*/
n = -qilog2(epsilon);
if (n < 0) {
return qlink(&_qzero_);
}
/*
* compute cosine and sine
*/
qsincos(q, n + 2, &sin, &cos);
tcos = qmappr(cos, epsilon, conf->triground);
qfree(cos);
tsin = qmappr(sin, epsilon, conf->triground);
qfree(sin);
res = qqadd(tcos, tsin);
qfree(tcos);
qfree(tsin);
/*
* return trigonometric result
*/
return res;
}
Reference in New Issue View Git Blame Copy Permalink