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

add new versin and vercos builtin functions

Browse Source 
Added new versin(x, [,eps]) for versed sine and vercos(x, [,eps])
for versed cosine.

Updated trig help files.
This commit is contained in:
Landon Curt Noll
2023-09-01 17:38:14 -07:00
parent 1c839dfede
commit b0a48a2b70
18 changed files with 493 additions and 52 deletions
Download Patch File Download Diff File Expand all files Collapse all files

5
CHANGES
View File

@@ -89,7 +89,7 @@ The following are the changes from calc version 2.14.3.5 to date:
Setting an invalid epsilon via the epsilon(value) or confiv("epsilon",
value) triggers an error. The epsilon value must be: 0 < epsilon < 1.
Add new logn(x, n [,eps]) builtin to compute logarithms to base n.
Added new logn(x, n [,eps]) builtin to compute logarithms to base n.
Verify that eps arguments (error tolerance arguments that override
the default epsilon value) to builtin functions have proper values.
@@ -99,6 +99,9 @@ The following are the changes from calc version 2.14.3.5 to date:
Document in help files for builtin functions that take eps arguments,
the LIMIT range for such eps values.
Added new versin(x, [,eps]) for versed sine and vercos(x, [,eps])
for versed cosine.
The following are the changes from calc version 2.14.3.4 to 2.14.3.5:

57
cal/regress.cal
View File

@@ -3498,15 +3498,68 @@ print '050: read -once test3400';
define test_trig()
{
local tnum; /* test number */
local pi; /* pi to 1e-20 precision */
local i;
print '3400: Beginning test_trig';
/* test 3401-3407 */
tnum = test3400(1, 3401);
vrfy(tnum == 3407, '3407: tnum == 3407');
vrfy(tnum++ == 3407, '3407: tnum == 3407');
print '3438: Ending test_trig';
/* test versed sine */
pi = pi(1e-20);
vrfy(round(versin(0.2, 1e-10), 10) == 0.0199334222,
strcat(str(tnum++),
': round(versin(0.2, 1e-10), 10) == 0.0199334222'));
vrfy(round(versin(3/7, 1e-10), 10) == 0.0904396483,
strcat(str(tnum++),
': round(versin(3/7, 1e-10), 10) == 0.0904396483'));
vrfy(round(versin(-31, 1e-10), 10) == 0.0852576422,
strcat(str(tnum++),
': round(versin(-31, 1e-10), 10) == 0.0852576422'));
vrfy(versin(pi/3, 1e-10) == 0.5,
strcat(str(tnum++), ': versin(pi/3, 1e-10) == 0.5'));
vrfy(versin(pi/2, 1e-10) == 1,
strcat(str(tnum++), ': versin(pi/2, 1e-10) == 1'));
vrfy(versin(pi, 1e-10) == 2,
strcat(str(tnum++), ': versin(pi, 1e-10) == 2'));
vrfy(versin(3*pi/2, 1e-10) == 1,
strcat(str(tnum++), ': versin(3*pi/2, 1e-10) == 1'));
vrfy(round(versin(1, 1e-10), 10) == 0.4596976941,
strcat(str(tnum++),
': round(versin(1, 1e-10), 10) == 0.4596976941'));
vrfy(round(versin(2 + 3i, 1e-10), 10) == 5.189625691+9.1092278938i,
strcat(str(tnum++),
': round(versin(2 + 3i, 1e-10), 10) == 5.189625691+9.1092278938i'));
/* test versed cosine */
pi = pi(1e-20);
vrfy(round(vercos(0.2, 1e-10), 10) == 0.8013306692,
strcat(str(tnum++),
': round(vercos(0.2, 1e-10), 10) == 0.8013306692'));
vrfy(round(vercos(3/7, 1e-10), 10) == 0.584428145,
strcat(str(tnum++),
': round(vercos(3/7, 1e-10), 10) == 0.584428145'));
vrfy(round(vercos(-31, 1e-10), 10) == 0.5959623547,
strcat(str(tnum++),
': round(vercos(-31, 1e-10), 10) == 0.5959623547'));
vrfy(vercos(pi/6, 1e-10) == 0.5,
strcat(str(tnum++), ': vercos(pi/6, 1e-10) == 0.5'));
vrfy(vercos(pi/2, 1e-10) == 0,
strcat(str(tnum++), ': vercos(pi/2, 1e-10) == 0'));
vrfy(vercos(pi, 1e-10) == 1,
strcat(str(tnum++), ': vercos(pi, 1e-10) == 1'));
vrfy(vercos(3*pi/2, 1e-10) == 2,
strcat(str(tnum++), ': vercos(3*pi/2, 1e-10) == 2'));
vrfy(round(vercos(1, 1e-10), 10) == 0.1585290152,
strcat(str(tnum++),
': round(vercos(1, 1e-10), 10) == 0.1585290152'));
vrfy(round(vercos(2 + 3i, 1e-10), 10) == -8.1544991469+4.16890696i,
strcat(str(tnum++),
': round(vercos(2 + 3i, 1e-10), 10) == -8.1544991469+4.16890696i'));
print strcat(str(tnum++), ': Ending test_trig');
}
print '051: parsed test_trig()';

6
calcerr.tbl
View File

@@ -550,3 +550,9 @@ E_LOGN_2 Non-numeric first argument for logn
E_LOGN_3 Cannot calculate logn for this value
E_LOGN_4 Cannot calculate logn for this log base
E_LOGN_5 Non-numeric second argument for logn
E_VERSIN1 Bad epsilon for versin
E_VERSIN2 Bad first argument for versin
E_VERSIN3 Too-large im(argument) for versin
E_VERCOS1 Bad epsilon for vercos
E_VERCOS2 Bad first argument for vercos
E_VERCOS3 Too-large im(argument) for vercos

2
cmath.h
View File

@@ -118,6 +118,8 @@ E_FUNC COMPLEX *c_asech(COMPLEX *c, NUMBER *epsilon);
E_FUNC COMPLEX *c_acsch(COMPLEX *c, NUMBER *epsilon);
E_FUNC COMPLEX *c_gd(COMPLEX *c, NUMBER *epsilon);
E_FUNC COMPLEX *c_agd(COMPLEX *c, NUMBER *epsilon);
E_FUNC COMPLEX *c_versin(COMPLEX *c, NUMBER *epsilon);
E_FUNC COMPLEX *c_vercos(COMPLEX *c, NUMBER *epsilon);

110
comfunc.c
View File

@@ -449,8 +449,8 @@ c_exp(COMPLEX *c, NUMBER *epsilon)
NUMBER *sin, *cos, *tmp1, *tmp2, *epsilon1;
long k, n;
if (qiszero(epsilon)) {
math_error("Zero epsilon for cexp");
if (check_epsilon(epsilon) == false) {
math_error("Invalid epsilon value for complex exp");
not_reached();
}
if (cisreal(c)) {
@@ -505,8 +505,8 @@ c_ln(COMPLEX *c, NUMBER *epsilon)
COMPLEX *r;
NUMBER *a2b2, *tmp1, *tmp2, *epsilon1;
if (ciszero(c)) {
math_error("logarithm of zero");
if (check_epsilon(epsilon) == false) {
math_error("Invalid epsilon value for complex ln");
not_reached();
}
if (cisone(c))
@@ -658,8 +658,8 @@ c_cos(COMPLEX *c, NUMBER *epsilon)
long n;
bool neg;
if (qiszero(epsilon)) {
math_error("Zero epsilon for ccos");
if (check_epsilon(epsilon) == false) {
math_error("Invalid epsilon value for complex cos");
not_reached();
}
n = qilog2(epsilon);
@@ -708,12 +708,13 @@ c_sin(COMPLEX *c, NUMBER *epsilon)
long n;
bool neg;
if (qiszero(epsilon)) {
math_error("Zero epsilon for csin");
if (check_epsilon(epsilon) == false) {
math_error("Invalid epsilon value for complex sin");
not_reached();
}
if (ciszero(c))
if (ciszero(c)) {
return clink(&_czero_);
}
n = qilog2(epsilon);
ctmp1 = comalloc();
neg = qisneg(c->imag);
@@ -725,8 +726,9 @@ c_sin(COMPLEX *c, NUMBER *epsilon)
ctmp2 = c_exp(ctmp1, epsilon1);
comfree(ctmp1);
qfree(epsilon1);
if (ctmp2 == NULL)
if (ctmp2 == NULL) {
return NULL;
}
if (ciszero(ctmp2)) {
comfree(ctmp2);
return clink(&_czero_);
@@ -1131,8 +1133,8 @@ c_polar(NUMBER *q1, NUMBER *q2, NUMBER *epsilon)
NUMBER *tmp, *cos, *sin;
long m, n;
if (qiszero(epsilon)) {
math_error("Zero epsilon for cpolar");
if (check_epsilon(epsilon) == false) {
math_error("Invalid epsilon value for complex polar");
not_reached();
}
if (qiszero(q1))
@@ -1173,13 +1175,13 @@ c_power(COMPLEX *c1, COMPLEX *c2, NUMBER *epsilon)
long k1, k2, k, m1, m2, m, n;
NUMBER *a2b2, *qtmp1, *qtmp2, *epsilon1;
if (qiszero(epsilon)) {
math_error("Zero epsilon for cpower");
if (check_epsilon(epsilon) == false) {
math_error("Invalid epsilon value for complex power");
not_reached();
}
if (ciszero(c1)) {
if (cisreal(c2) && qisneg(c2->real)) {
math_error ("Non-positive real exponent of zero");
math_error ("Non-positive real exponent of zero for complex power");
not_reached();
}
return clink(&_czero_);
@@ -1323,3 +1325,81 @@ c_ilog(COMPLEX *c, ZVALUE base)
qfree(qr);
return qi;
}
/*
* Calculate the complex versed sine within the specified accuracy.
*
* This uses the formula:
*
* versin(x) = 1 - cos(x)
*/
COMPLEX *
c_versin(COMPLEX *c, NUMBER *epsilon)
{
COMPLEX *r; /* return COMPLEX value */
COMPLEX *ctmp; /* complex cos(c) */
/*
* firewall
*/
if (c == NULL) {
math_error("%s: c is NULL", __func__);
not_reached();
}
if (check_epsilon(epsilon) == false) {
math_error("Invalid epsilon value for complex versin");
not_reached();
}
/*
* return r = 1 - cos(x)
*/
ctmp = c_cos(c, epsilon);
if (ctmp == NULL) {
math_error("Failed to compute complex cos for complex versin");
not_reached();
}
r = c_sub(&_cone_, ctmp);
comfree(ctmp);
return r;
}
/*
* Calculate the complex versed cosine within the specified accuracy.
*
* This uses the formula:
*
* vercos(x) = 1 - sin(x)
*/
COMPLEX *
c_vercos(COMPLEX *c, NUMBER *epsilon)
{
COMPLEX *r; /* return COMPLEX value */
COMPLEX *ctmp; /* complex sin(c) */
/*
* firewall
*/
if (c == NULL) {
math_error("%s: c is NULL", __func__);
not_reached();
}
if (check_epsilon(epsilon) == false) {
math_error("Invalid epsilon value for complex vercos");
not_reached();
}
/*
* return r = 1 - sin(x)
*/
ctmp = c_sin(c, epsilon);
if (ctmp == NULL) {
math_error("Failed to compute complex sin for complex vercos");
not_reached();
}
r = c_sub(&_cone_, ctmp);
comfree(ctmp);
return r;
}

107
func.c
View File

@@ -3219,6 +3219,7 @@ f_sinh(int count, VALUE **vals)
return result;
}
S_FUNC VALUE
f_cosh(int count, VALUE **vals)
{
@@ -3485,6 +3486,108 @@ f_csch(int count, VALUE **vals)
}
S_FUNC VALUE
f_versin(int count, VALUE **vals)
{
VALUE result;
COMPLEX *c;
NUMBER *eps;
/* initialize VALUE */
result.v_subtype = V_NOSUBTYPE;
/*
* set error tolerance for builtin function
*
* Use eps VALUE arg if given and value is in a valid range.
*/
eps = conf->epsilon;
if (count == 2) {
if (verify_eps(vals[1]) == false) {
return error_value(E_VERSIN1);
}
eps = vals[1]->v_num;
}
/*
* compute sine to a given error tolerance
*/
switch (vals[0]->v_type) {
case V_NUM:
result.v_num = qversin(vals[0]->v_num, eps);
result.v_type = V_NUM;
break;
case V_COM:
c = c_versin(vals[0]->v_com, eps);
if (c == NULL) {
return error_value(E_VERSIN3);
}
result.v_com = c;
result.v_type = V_COM;
if (cisreal(c)) {
result.v_num = qlink(c->real);
result.v_type = V_NUM;
comfree(c);
}
break;
default:
return error_value(E_VERSIN2);
}
return result;
}
S_FUNC VALUE
f_vercos(int count, VALUE **vals)
{
VALUE result;
COMPLEX *c;
NUMBER *eps;
/* initialize VALUE */
result.v_subtype = V_NOSUBTYPE;
/*
* set error tolerance for builtin function
*
* Use eps VALUE arg if given and value is in a valid range.
*/
eps = conf->epsilon;
if (count == 2) {
if (verify_eps(vals[1]) == false) {
return error_value(E_VERCOS1);
}
eps = vals[1]->v_num;
}
/*
* compute cosinr to a given error tolerance
*/
switch (vals[0]->v_type) {
case V_NUM:
result.v_num = qvercos(vals[0]->v_num, eps);
result.v_type = V_NUM;
break;
case V_COM:
c = c_vercos(vals[0]->v_com, eps);
if (c == NULL) {
return error_value(E_VERCOS3);
}
result.v_com = c;
result.v_type = V_COM;
if (cisreal(c)) {
result.v_num = qlink(c->real);
result.v_type = V_NUM;
comfree(c);
}
break;
default:
return error_value(E_VERCOS2);
}
return result;
}
S_FUNC VALUE
f_atan(int count, VALUE **vals)
{
@@ -11313,6 +11416,10 @@ STATIC CONST struct builtin builtins[] = {
"unget char read from file"},
{"usertime", 0, 0, 0, OP_NOP, {.numfunc_0 = f_usertime}, {.null = NULL},
"user mode CPU time in seconds"},
{"vercos", 1, 2, 0, OP_NOP, {.null = NULL}, {.valfunc_cnt = f_vercos},
"versed cosine of value a within accuracy b"},
{"versin", 1, 2, 0, OP_NOP, {.null = NULL}, {.valfunc_cnt = f_versin},
"versed sine of value a within accuracy b"},
{"version", 0, 0, 0, OP_NOP, {.null = NULL}, {.valfunc_0 = f_version},
"calc version string"},
{"xor", 1, IN, 0, OP_NOP, {.null = NULL}, {.valfunc_cnt = f_xor},

4
help/Makefile
View File

@@ -228,8 +228,8 @@ DETAIL_HELP= abs access acos acosh acot acoth acsc acsch address agd \
sin sinh size sizeof sleep sort sqrt srand srandom ssq stoponerror str \
strcasecmp strcat strcmp strcpy strerror strlen strncasecmp strncmp \
strncpy strpos strprintf strscan strscanf strtolower strtoupper substr \
sum swap system systime tail tan tanh test time trunc usertime version \
xor
sum swap system systime tail tan tanh test time trunc usertime vercos \
versin version xor
# This list is of files that are clones of DETAIL_HELP files. They are
# built from DETAIL_HELP files.

1
help/cos
View File

@@ -34,6 +34,7 @@ LINK LIBRARY
SEE ALSO
sin, tan, sec, csc, cot, epsilon
versin, vercos
## Copyright (C) 1999,2021,2023 Landon Curt Noll
##

2
help/cot
View File

@@ -23,9 +23,11 @@ LIMITS
LINK LIBRARY
NUMBER *qcot(NUMBER *x, NUMBER *eps)
COMPLEX *c_acot(COMPLEX *c, NUMBER *eps);
SEE ALSO
sin, cos, tan, sec, csc, epsilon
versin, vercos
## Copyright (C) 1999,2021,2023 Landon Curt Noll
##

1
help/csc
View File

@@ -26,6 +26,7 @@ LINK LIBRARY
SEE ALSO
sin, cos, tan, sec, cot, epsilon
versin, vercos
## Copyright (C) 1999,2023 Landon Curt Noll
##

1
help/sec
View File

@@ -27,6 +27,7 @@ LINK LIBRARY
SEE ALSO
sin, cos, tan, csc, cot, epsilon
versin, vercos
## Copyright (C) 1999,2023 Landon Curt Noll
##

1
help/sin
View File

@@ -34,6 +34,7 @@ LINK LIBRARY
SEE ALSO
cos, tan, sec, csc, cot, epsilon
versin, vercos
## Copyright (C) 1999,2021,2023 Landon Curt Noll
##

2
help/tan
View File

@@ -24,9 +24,11 @@ LIMITS
LINK LIBRARY
NUMBER *qtan(NUMBER *x, NUMBER *eps)
COMPLEX *c_atan(COMPLEX *c, NUMBER *eps);
SEE ALSO
sin, cos, sec, csc, cot, epsilon
versin, vercos
## Copyright (C) 1999,2023 Landon Curt Noll
##

67
help/vercos Normal file
View File

@@ -0,0 +1,67 @@
NAME
vercos - versed cosine
SYNOPSIS
vercos(x [,eps])
TYPES
x number (real or complex)
eps 0 < real < 1, defaults to epsilon()
return number
DESCRIPTION
Calculate the versed cosine of x to a multiple of eps with error less in
absolute value than .75 * eps.
The versed cosine function is sometimes called coversin, sometimes called cvs,
may be defined as:
vercos(x) = 1 - sin(x)
EXAMPLE
; print vercos(0.2), vercos(3/7), vercos(-31)
0.80133066920493878454 0.58442814500694799193 0.59596235467693499395
; print vercos(1, 1e-5), vercos(1, 1e-10), vercos(1, 1e-15), vercos(1, 1e-20)
0.15853 0.1585290152 0.158529015192104 0.15852901519210349335
; print vercos(2 + 3i, 1e-5), vercos(2 + 3i, 1e-10)
-8.1545+4.16891i -8.1544991469+4.16890696i
; pi = pi(1e-20)
; print vercos(pi/6, 1e-10), vercos(pi/2, 1e-10), vercos(pi, 1e-10), vercos(3*pi/2, 1e-10)
0.5 0 1 2
LIMITS
0 < eps < 1
LINK LIBRARY
NUMBER *qvercos(NUMBER *x, NUMBER *eps)
COMPLEX *c_vercos(COMPLEX *x, NUMBER *eps)
SEE ALSO
sin, cos, tan, sec, csc, cot, epsilon
versin
## Copyright (C) 2023 Landon Curt Noll
##
## 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: 2023/08/31 23:07:08
## File existed as early as: 2023
##
## chongo <was here> /\oo/\ http://www.isthe.com/chongo/
## Share and enjoy! :-) http://www.isthe.com/chongo/tech/comp/calc/

67
help/versin Normal file
View File

@@ -0,0 +1,67 @@
NAME
versin - versed sine
SYNOPSIS
versin(x [,eps])
TYPES
x number (real or complex)
eps 0 < real < 1, defaults to epsilon()
return number
DESCRIPTION
Calculate the versed sine of x to a multiple of eps with error less in
absolute value than .75 * eps.
The versed sine function is sometimes called vers, sometimes called ver,
may be defined as:
versin(x) = 1 - cos(x)
EXAMPLE
; print versin(0.2), versin(3/7), versin(-31)
0.01993342215875836888 0.09043964832583332597 0.08525764219546872104
; print versin(1, 1e-5), versin(1, 1e-10), versin(1, 1e-15), versin(1, 1e-20)
0.4597 0.4596976941 0.45969769413186 0.4596976941318602826
; print versin(2 + 3i, 1e-5), versin(2 + 3i, 1e-10)
5.18963+9.10923i 5.189625691+9.1092278938i
; pi = pi(1e-20)
; print versin(pi/3, 1e-10), versin(pi/2, 1e-10), versin(pi, 1e-10), versin(3*pi/2, 1e-10)
0.5 1 2 1
LIMITS
0 < eps < 1
LINK LIBRARY
NUMBER *qversin(NUMBER *x, NUMBER *eps)
COMPLEX *c_versin(COMPLEX *x, NUMBER *eps)
SEE ALSO
sin, cos, tan, sec, csc, cot, epsilon
vercos
## Copyright (C) 2023 Landon Curt Noll
##
## 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: 2023/08/31 23:05:28
## File existed as early as: 2023
##
## chongo <was here> /\oo/\ http://www.isthe.com/chongo/
## Share and enjoy! :-) http://www.isthe.com/chongo/tech/comp/calc/

50
qfunc.c
View File

@@ -45,7 +45,7 @@ STATIC long E_num;
/*
* verify_epsilon - verify that 0 < epsilon < 1
* check_epsilon - verify that 0 < epsilon < 1
*
* given:
* q epsilon or eps argument
@@ -65,27 +65,6 @@ check_epsilon(NUMBER *q)
}
/*
* verify_epsilon - verify that 0 < epsilon < 1
*
* This function is called via the OP_SETEPSILON op code, config("epsilon").
*
* If all is well, this function just returns. If the arg passed is
* out of range, then a math_error() is triggered causing this function
* to not return.
*/
void
verify_epsilon(NUMBER *q)
{
/* verify that 0 < epsilon < 1 */
if (check_epsilon(q) == false) {
math_error("Invalid value for epsilon: must be: 0 < epsilon < 1");
not_reached();
}
return;
}
/*
* Set the default epsilon for approximate calculations.
* This must be greater than zero.
@@ -98,7 +77,18 @@ setepsilon(NUMBER *q)
{
NUMBER *old;
verify_epsilon(q);
/*
* firewall
*/
if (q == NULL) {
math_error("q is NULL for %s", __func__);
not_reached();
}
if (check_epsilon(q) == false) {
math_error("Invalid value for epsilon: must be: 0 < epsilon < 1");
not_reached();
}
old = conf->epsilon;
conf->epsilonprec = qprecision(q);
conf->epsilon = qlink(q);
@@ -321,7 +311,7 @@ qlegtoleg(NUMBER *q, NUMBER *epsilon, bool wantneg)
ZVALUE num;
if (qiszero(epsilon)) {
math_error("Zero epsilon value for legtoleg");
math_error("Zero epsilon value for ltol");
not_reached();
}
if (qisunit(q))
@@ -334,7 +324,7 @@ qlegtoleg(NUMBER *q, NUMBER *epsilon, bool wantneg)
num = q->num;
num.sign = 0;
if (zrel(num, q->den) >= 0) {
math_error("Leg too large in legtoleg");
math_error("Leg too large for ltol");
not_reached();
}
qtmp1 = qsquare(q);
@@ -369,7 +359,7 @@ qsqrt(NUMBER *q1, NUMBER *epsilon, long rnd)
int sign;
if (qisneg(q1)) {
math_error("Square root of negative number");
math_error("Square root of negative number for qsqrt");
not_reached();
}
if (qiszero(q1))
@@ -475,7 +465,7 @@ qisqrt(NUMBER *q)
ZVALUE tmp;
if (qisneg(q)) {
math_error("Square root of negative number");
math_error("Square root of negative number for isqrt");
not_reached();
}
if (qiszero(q))
@@ -1910,15 +1900,15 @@ qispowerof2(NUMBER *q, NUMBER **qlog2)
/* firewall */
if (q == NULL) {
math_error("%s: q NULL", __func__);
math_error("%s: q is NULL", __func__);
not_reached();
}
if (qlog2 == NULL) {
math_error("%s: qlog2 NULL", __func__);
math_error("%s: qlog2 is NULL", __func__);
not_reached();
}
if (*qlog2 == NULL) {
math_error("%s: *qlog2 NULL", __func__);
math_error("%s: *qlog2 is NULL", __func__);
not_reached();
}

3
qmath.h
View File

@@ -168,7 +168,6 @@ E_FUNC long qplaces(NUMBER *q, ZVALUE base);
E_FUNC long qdecplaces(NUMBER *q);
E_FUNC long qdigits(NUMBER *q, ZVALUE base);
E_FUNC bool check_epsilon(NUMBER *q);
E_FUNC void verify_epsilon(NUMBER *q);
E_FUNC void setepsilon(NUMBER *q);
E_FUNC NUMBER *qbitvalue(long i);
E_FUNC NUMBER *qtenpow(long i);
@@ -223,6 +222,8 @@ E_FUNC NUMBER *qbern(ZVALUE z);
E_FUNC void qfreebern(void);
E_FUNC NUMBER *qeuler(ZVALUE z);
E_FUNC void qfreeeuler(void);
E_FUNC NUMBER *qversin(NUMBER *q, NUMBER *epsilon);
E_FUNC NUMBER *qvercos(NUMBER *q, NUMBER *epsilon);
/*

59
qtrans.c
View File

@@ -229,7 +229,6 @@ qcos(NUMBER *q, NUMBER *epsilon)
}
/*
* This calls qsincos() and discards the value of cos.
*/
@@ -1973,3 +1972,61 @@ qacoth(NUMBER *q, NUMBER *epsilon)
qfree(tmp);
return res;
}
/*
* versed sine - this calls qsincos() and discards the value of sin.
*
* This uses the formula: versin(x) = 1 - cos(x).
*/
NUMBER *
qversin(NUMBER *q, NUMBER *epsilon)
{
NUMBER *sin, *cos, *res;
NUMBER *versin;
long n;
if (qiszero(epsilon)) {
math_error("Zero epsilon value for %s", __func__);
not_reached();
}
n = -qilog2(epsilon);
if (qiszero(q) || n < 0)
return qlink(&_qzero_);
qsincos(q, n + 2, &sin, &cos);
qfree(sin);
versin = qsub(&_qone_, cos);
qfree(cos);
res = qmappr(versin, epsilon, 24);
qfree(versin);
return res;
}
/*
* versed cosine - this calls qsincos() and discards the value of cos.
*
* This uses the formula: vercos(x) = 1 - sin(x).
*/
NUMBER *
qvercos(NUMBER *q, NUMBER *epsilon)
{
NUMBER *sin, *cos, *res;
NUMBER *vercos;
long n;
if (qiszero(epsilon)) {
math_error("Zero epsilon value for %s", __func__);
not_reached();
}
n = -qilog2(epsilon);
if (qiszero(q) || n < 0)
return qlink(&_qzero_);
qsincos(q, n + 2, &sin, &cos);
qfree(cos);
vercos = qsub(&_qone_, sin);
qfree(sin);
res = qmappr(vercos, epsilon, 24);
qfree(vercos);
return res;
}