From b0a48a2b70af0aff123b8f51dee7a86e0f99b1ce Mon Sep 17 00:00:00 2001 From: Landon Curt Noll Date: Fri, 1 Sep 2023 17:38:14 -0700 Subject: [PATCH] add new versin and vercos builtin functions Added new versin(x, [,eps]) for versed sine and vercos(x, [,eps]) for versed cosine. Updated trig help files. --- CHANGES | 5 ++- cal/regress.cal | 57 ++++++++++++++++++++++++- calcerr.tbl | 6 +++ cmath.h | 2 + comfunc.c | 110 +++++++++++++++++++++++++++++++++++++++++------- func.c | 107 ++++++++++++++++++++++++++++++++++++++++++++++ help/Makefile | 4 +- help/cos | 1 + help/cot | 2 + help/csc | 1 + help/sec | 1 + help/sin | 1 + help/tan | 2 + help/vercos | 67 +++++++++++++++++++++++++++++ help/versin | 67 +++++++++++++++++++++++++++++ qfunc.c | 50 +++++++++------------- qmath.h | 3 +- qtrans.c | 59 +++++++++++++++++++++++++- 18 files changed, 493 insertions(+), 52 deletions(-) create mode 100644 help/vercos create mode 100644 help/versin diff --git a/CHANGES b/CHANGES index 3b84686..0cce4d2 100644 --- a/CHANGES +++ b/CHANGES @@ -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: diff --git a/cal/regress.cal b/cal/regress.cal index 54e5a10..8044520 100644 --- a/cal/regress.cal +++ b/cal/regress.cal @@ -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()'; diff --git a/calcerr.tbl b/calcerr.tbl index 5d98e8c..d5c1560 100644 --- a/calcerr.tbl +++ b/calcerr.tbl @@ -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 diff --git a/cmath.h b/cmath.h index 1996fec..2abeefc 100644 --- a/cmath.h +++ b/cmath.h @@ -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); diff --git a/comfunc.c b/comfunc.c index d93d551..50dec9a 100644 --- a/comfunc.c +++ b/comfunc.c @@ -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; +} diff --git a/func.c b/func.c index 6dcd44a..2f8a990 100644 --- a/func.c +++ b/func.c @@ -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}, diff --git a/help/Makefile b/help/Makefile index d65fc89..d9d6e93 100644 --- a/help/Makefile +++ b/help/Makefile @@ -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. diff --git a/help/cos b/help/cos index 2376cb3..6a7c921 100644 --- a/help/cos +++ b/help/cos @@ -34,6 +34,7 @@ LINK LIBRARY SEE ALSO sin, tan, sec, csc, cot, epsilon + versin, vercos ## Copyright (C) 1999,2021,2023 Landon Curt Noll ## diff --git a/help/cot b/help/cot index d6a47c2..414d758 100644 --- a/help/cot +++ b/help/cot @@ -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 ## diff --git a/help/csc b/help/csc index 9c93928..4b95241 100644 --- a/help/csc +++ b/help/csc @@ -26,6 +26,7 @@ LINK LIBRARY SEE ALSO sin, cos, tan, sec, cot, epsilon + versin, vercos ## Copyright (C) 1999,2023 Landon Curt Noll ## diff --git a/help/sec b/help/sec index eb7696e..b563824 100644 --- a/help/sec +++ b/help/sec @@ -27,6 +27,7 @@ LINK LIBRARY SEE ALSO sin, cos, tan, csc, cot, epsilon + versin, vercos ## Copyright (C) 1999,2023 Landon Curt Noll ## diff --git a/help/sin b/help/sin index 97b102e..e8bfade 100644 --- a/help/sin +++ b/help/sin @@ -34,6 +34,7 @@ LINK LIBRARY SEE ALSO cos, tan, sec, csc, cot, epsilon + versin, vercos ## Copyright (C) 1999,2021,2023 Landon Curt Noll ## diff --git a/help/tan b/help/tan index d66166d..5b17f82 100644 --- a/help/tan +++ b/help/tan @@ -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 ## diff --git a/help/vercos b/help/vercos new file mode 100644 index 0000000..4cac2e2 --- /dev/null +++ b/help/vercos @@ -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 /\oo/\ http://www.isthe.com/chongo/ +## Share and enjoy! :-) http://www.isthe.com/chongo/tech/comp/calc/ diff --git a/help/versin b/help/versin new file mode 100644 index 0000000..9ae00d6 --- /dev/null +++ b/help/versin @@ -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 /\oo/\ http://www.isthe.com/chongo/ +## Share and enjoy! :-) http://www.isthe.com/chongo/tech/comp/calc/ diff --git a/qfunc.c b/qfunc.c index d5e43d5..6a8f1c2 100644 --- a/qfunc.c +++ b/qfunc.c @@ -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(); } diff --git a/qmath.h b/qmath.h index 99c02aa..21f9532 100644 --- a/qmath.h +++ b/qmath.h @@ -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); /* diff --git a/qtrans.c b/qtrans.c index c977739..179a164 100644 --- a/qtrans.c +++ b/qtrans.c @@ -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; +}