|
Packit |
70b277 |
/* This file is part of GNU bc.
|
|
Packit |
70b277 |
|
|
Packit |
70b277 |
Copyright (C) 1991-1994, 1997, 2006, 2008, 2012-2017 Free Software Foundation, Inc.
|
|
Packit |
70b277 |
|
|
Packit |
70b277 |
This program is free software; you can redistribute it and/or modify
|
|
Packit |
70b277 |
it under the terms of the GNU General Public License as published by
|
|
Packit |
70b277 |
the Free Software Foundation; either version 3 of the License , or
|
|
Packit |
70b277 |
(at your option) any later version.
|
|
Packit |
70b277 |
|
|
Packit |
70b277 |
This program is distributed in the hope that it will be useful,
|
|
Packit |
70b277 |
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
Packit |
70b277 |
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
Packit |
70b277 |
GNU General Public License for more details.
|
|
Packit |
70b277 |
|
|
Packit |
70b277 |
You should have received a copy of the GNU General Public License
|
|
Packit |
70b277 |
along with this program; see the file COPYING. If not, see
|
|
Packit |
70b277 |
<http://www.gnu.org/licenses>.
|
|
Packit |
70b277 |
|
|
Packit |
70b277 |
You may contact the author by:
|
|
Packit |
70b277 |
e-mail: philnelson@acm.org
|
|
Packit |
70b277 |
us-mail: Philip A. Nelson
|
|
Packit |
70b277 |
Computer Science Department, 9062
|
|
Packit |
70b277 |
Western Washington University
|
|
Packit |
70b277 |
Bellingham, WA 98226-9062
|
|
Packit |
70b277 |
|
|
Packit |
70b277 |
*************************************************************************/
|
|
Packit |
70b277 |
|
|
Packit |
70b277 |
/* libmath.b for bc. */
|
|
Packit |
70b277 |
|
|
Packit |
70b277 |
scale = 20
|
|
Packit |
70b277 |
|
|
Packit |
70b277 |
/* Uses the fact that e^x = (e^(x/2))^2
|
|
Packit |
70b277 |
When x is small enough, we use the series:
|
|
Packit |
70b277 |
e^x = 1 + x + x^2/2! + x^3/3! + ...
|
|
Packit |
70b277 |
*/
|
|
Packit |
70b277 |
|
|
Packit |
70b277 |
define e(x) {
|
|
Packit |
70b277 |
auto a, b, d, e, f, i, m, n, v, z
|
|
Packit |
70b277 |
|
|
Packit |
70b277 |
/* a - holds x^y of x^y/y! */
|
|
Packit |
70b277 |
/* d - holds y! */
|
|
Packit |
70b277 |
/* e - is the value x^y/y! */
|
|
Packit |
70b277 |
/* v - is the sum of the e's */
|
|
Packit |
70b277 |
/* f - number of times x was divided by 2. */
|
|
Packit |
70b277 |
/* m - is 1 if x was minus. */
|
|
Packit |
70b277 |
/* i - iteration count. */
|
|
Packit |
70b277 |
/* n - the scale to compute the sum. */
|
|
Packit |
70b277 |
/* z - orignal scale. */
|
|
Packit |
70b277 |
/* b - holds the original ibase. */
|
|
Packit |
70b277 |
|
|
Packit |
70b277 |
/* Non base 10 ibase? */
|
|
Packit |
70b277 |
if (ibase != A) {
|
|
Packit |
70b277 |
b = ibase;
|
|
Packit |
70b277 |
ibase = A;
|
|
Packit |
70b277 |
v = e(x);
|
|
Packit |
70b277 |
ibase = b;
|
|
Packit |
70b277 |
return (v);
|
|
Packit |
70b277 |
}
|
|
Packit |
70b277 |
|
|
Packit |
70b277 |
/* Check the sign of x. */
|
|
Packit |
70b277 |
if (x<0) {
|
|
Packit |
70b277 |
m = 1
|
|
Packit |
70b277 |
x = -x
|
|
Packit |
70b277 |
}
|
|
Packit |
70b277 |
|
|
Packit |
70b277 |
/* Precondition x. */
|
|
Packit |
70b277 |
z = scale;
|
|
Packit |
70b277 |
n = 6 + z + .44*x;
|
|
Packit |
70b277 |
scale = scale(x)+1;
|
|
Packit |
70b277 |
while (x > 1) {
|
|
Packit |
70b277 |
f += 1;
|
|
Packit |
70b277 |
x /= 2;
|
|
Packit |
70b277 |
scale += 1;
|
|
Packit |
70b277 |
}
|
|
Packit |
70b277 |
|
|
Packit |
70b277 |
/* Initialize the variables. */
|
|
Packit |
70b277 |
scale = n;
|
|
Packit |
70b277 |
v = 1+x
|
|
Packit |
70b277 |
a = x
|
|
Packit |
70b277 |
d = 1
|
|
Packit |
70b277 |
|
|
Packit |
70b277 |
for (i=2; 1; i++) {
|
|
Packit |
70b277 |
e = (a *= x) / (d *= i)
|
|
Packit |
70b277 |
if (e == 0) {
|
|
Packit |
70b277 |
if (f>0) while (f--) v = v*v;
|
|
Packit |
70b277 |
scale = z
|
|
Packit |
70b277 |
if (m) return (1/v);
|
|
Packit |
70b277 |
return (v/1);
|
|
Packit |
70b277 |
}
|
|
Packit |
70b277 |
v += e
|
|
Packit |
70b277 |
}
|
|
Packit |
70b277 |
}
|
|
Packit |
70b277 |
|
|
Packit |
70b277 |
/* Natural log. Uses the fact that ln(x^2) = 2*ln(x)
|
|
Packit |
70b277 |
The series used is:
|
|
Packit |
70b277 |
ln(x) = 2(a+a^3/3+a^5/5+...) where a=(x-1)/(x+1)
|
|
Packit |
70b277 |
*/
|
|
Packit |
70b277 |
|
|
Packit |
70b277 |
define l(x) {
|
|
Packit |
70b277 |
auto b, e, f, i, m, n, v, z
|
|
Packit |
70b277 |
|
|
Packit |
70b277 |
/* Non base 10 ibase? */
|
|
Packit |
70b277 |
if (ibase != A) {
|
|
Packit |
70b277 |
b = ibase;
|
|
Packit |
70b277 |
ibase = A;
|
|
Packit |
70b277 |
v = l(x);
|
|
Packit |
70b277 |
ibase = b;
|
|
Packit |
70b277 |
return (v);
|
|
Packit |
70b277 |
}
|
|
Packit |
70b277 |
|
|
Packit |
70b277 |
/* return something for the special case. */
|
|
Packit |
70b277 |
if (x <= 0) return ((1 - 10^scale)/1)
|
|
Packit |
70b277 |
|
|
Packit |
70b277 |
/* Precondition x to make .5 < x < 2.0. */
|
|
Packit |
70b277 |
z = scale;
|
|
Packit |
70b277 |
scale = 6 + scale;
|
|
Packit |
70b277 |
f = 2;
|
|
Packit |
70b277 |
i=0
|
|
Packit |
70b277 |
while (x >= 2) { /* for large numbers */
|
|
Packit |
70b277 |
f *= 2;
|
|
Packit |
70b277 |
x = sqrt(x);
|
|
Packit |
70b277 |
}
|
|
Packit |
70b277 |
while (x <= .5) { /* for small numbers */
|
|
Packit |
70b277 |
f *= 2;
|
|
Packit |
70b277 |
x = sqrt(x);
|
|
Packit |
70b277 |
}
|
|
Packit |
70b277 |
|
|
Packit |
70b277 |
/* Set up the loop. */
|
|
Packit |
70b277 |
v = n = (x-1)/(x+1)
|
|
Packit |
70b277 |
m = n*n
|
|
Packit |
70b277 |
|
|
Packit |
70b277 |
/* Sum the series. */
|
|
Packit |
70b277 |
for (i=3; 1; i+=2) {
|
|
Packit |
70b277 |
e = (n *= m) / i
|
|
Packit |
70b277 |
if (e == 0) {
|
|
Packit |
70b277 |
v = f*v
|
|
Packit |
70b277 |
scale = z
|
|
Packit |
70b277 |
return (v/1)
|
|
Packit |
70b277 |
}
|
|
Packit |
70b277 |
v += e
|
|
Packit |
70b277 |
}
|
|
Packit |
70b277 |
}
|
|
Packit |
70b277 |
|
|
Packit |
70b277 |
/* Sin(x) uses the standard series:
|
|
Packit |
70b277 |
sin(x) = x - x^3/3! + x^5/5! - x^7/7! ... */
|
|
Packit |
70b277 |
|
|
Packit |
70b277 |
define s(x) {
|
|
Packit |
70b277 |
auto b, e, i, m, n, s, v, z
|
|
Packit |
70b277 |
|
|
Packit |
70b277 |
/* Non base 10 ibase? */
|
|
Packit |
70b277 |
if (ibase != A) {
|
|
Packit |
70b277 |
b = ibase;
|
|
Packit |
70b277 |
ibase = A;
|
|
Packit |
70b277 |
v = s(x);
|
|
Packit |
70b277 |
ibase = b;
|
|
Packit |
70b277 |
return (v);
|
|
Packit |
70b277 |
}
|
|
Packit |
70b277 |
|
|
Packit |
70b277 |
/* precondition x. */
|
|
Packit |
70b277 |
z = scale
|
|
Packit |
70b277 |
scale = 1.1*z + 2;
|
|
Packit |
70b277 |
v = a(1)
|
|
Packit |
70b277 |
if (x < 0) {
|
|
Packit |
70b277 |
m = 1;
|
|
Packit |
70b277 |
x = -x;
|
|
Packit |
70b277 |
}
|
|
Packit |
70b277 |
scale = 0
|
|
Packit |
70b277 |
n = (x / v + 2 )/4
|
|
Packit |
70b277 |
x = x - 4*n*v
|
|
Packit |
70b277 |
if (n%2) x = -x
|
|
Packit |
70b277 |
|
|
Packit |
70b277 |
/* Do the loop. */
|
|
Packit |
70b277 |
scale = z + 2;
|
|
Packit |
70b277 |
v = e = x
|
|
Packit |
70b277 |
s = -x*x
|
|
Packit |
70b277 |
for (i=3; 1; i+=2) {
|
|
Packit |
70b277 |
e *= s/(i*(i-1))
|
|
Packit |
70b277 |
if (e == 0) {
|
|
Packit |
70b277 |
scale = z
|
|
Packit |
70b277 |
if (m) return (-v/1);
|
|
Packit |
70b277 |
return (v/1);
|
|
Packit |
70b277 |
}
|
|
Packit |
70b277 |
v += e
|
|
Packit |
70b277 |
}
|
|
Packit |
70b277 |
}
|
|
Packit |
70b277 |
|
|
Packit |
70b277 |
/* Cosine : cos(x) = sin(x+pi/2) */
|
|
Packit |
70b277 |
define c(x) {
|
|
Packit |
70b277 |
auto b, v, z;
|
|
Packit |
70b277 |
|
|
Packit |
70b277 |
/* Non base 10 ibase? */
|
|
Packit |
70b277 |
if (ibase != A) {
|
|
Packit |
70b277 |
b = ibase;
|
|
Packit |
70b277 |
ibase = A;
|
|
Packit |
70b277 |
v = c(x);
|
|
Packit |
70b277 |
ibase = b;
|
|
Packit |
70b277 |
return (v);
|
|
Packit |
70b277 |
}
|
|
Packit |
70b277 |
|
|
Packit |
70b277 |
z = scale;
|
|
Packit |
70b277 |
scale = scale*1.2;
|
|
Packit |
70b277 |
v = s(x+a(1)*2);
|
|
Packit |
70b277 |
scale = z;
|
|
Packit |
70b277 |
return (v/1);
|
|
Packit |
70b277 |
}
|
|
Packit |
70b277 |
|
|
Packit |
70b277 |
/* Arctan: Using the formula:
|
|
Packit |
70b277 |
atan(x) = atan(c) + atan((x-c)/(1+xc)) for a small c (.2 here)
|
|
Packit |
70b277 |
For under .2, use the series:
|
|
Packit |
70b277 |
atan(x) = x - x^3/3 + x^5/5 - x^7/7 + ... */
|
|
Packit |
70b277 |
|
|
Packit |
70b277 |
define a(x) {
|
|
Packit |
70b277 |
auto a, b, e, f, i, m, n, s, v, z
|
|
Packit |
70b277 |
|
|
Packit |
70b277 |
/* a is the value of a(.2) if it is needed. */
|
|
Packit |
70b277 |
/* f is the value to multiply by a in the return. */
|
|
Packit |
70b277 |
/* e is the value of the current term in the series. */
|
|
Packit |
70b277 |
/* v is the accumulated value of the series. */
|
|
Packit |
70b277 |
/* m is 1 or -1 depending on x (-x -> -1). results are divided by m. */
|
|
Packit |
70b277 |
/* i is the denominator value for series element. */
|
|
Packit |
70b277 |
/* n is the numerator value for the series element. */
|
|
Packit |
70b277 |
/* s is -x*x. */
|
|
Packit |
70b277 |
/* z is the saved user's scale. */
|
|
Packit |
70b277 |
|
|
Packit |
70b277 |
/* Non base 10 ibase? */
|
|
Packit |
70b277 |
if (ibase != A) {
|
|
Packit |
70b277 |
b = ibase;
|
|
Packit |
70b277 |
ibase = A;
|
|
Packit |
70b277 |
v = a(x);
|
|
Packit |
70b277 |
ibase = b;
|
|
Packit |
70b277 |
return (v);
|
|
Packit |
70b277 |
}
|
|
Packit |
70b277 |
|
|
Packit |
70b277 |
/* Negative x? */
|
|
Packit |
70b277 |
m = 1;
|
|
Packit |
70b277 |
if (x<0) {
|
|
Packit |
70b277 |
m = -1;
|
|
Packit |
70b277 |
x = -x;
|
|
Packit |
70b277 |
}
|
|
Packit |
70b277 |
|
|
Packit |
70b277 |
/* Special case and for fast answers */
|
|
Packit |
70b277 |
if (x==1) {
|
|
Packit |
70b277 |
if (scale <= 25) return (.7853981633974483096156608/m)
|
|
Packit |
70b277 |
if (scale <= 40) return (.7853981633974483096156608458198757210492/m)
|
|
Packit |
70b277 |
if (scale <= 60) \
|
|
Packit |
70b277 |
return (.785398163397448309615660845819875721049292349843776455243736/m)
|
|
Packit |
70b277 |
}
|
|
Packit |
70b277 |
if (x==.2) {
|
|
Packit |
70b277 |
if (scale <= 25) return (.1973955598498807583700497/m)
|
|
Packit |
70b277 |
if (scale <= 40) return (.1973955598498807583700497651947902934475/m)
|
|
Packit |
70b277 |
if (scale <= 60) \
|
|
Packit |
70b277 |
return (.197395559849880758370049765194790293447585103787852101517688/m)
|
|
Packit |
70b277 |
}
|
|
Packit |
70b277 |
|
|
Packit |
70b277 |
|
|
Packit |
70b277 |
/* Save the scale. */
|
|
Packit |
70b277 |
z = scale;
|
|
Packit |
70b277 |
|
|
Packit |
70b277 |
/* Note: a and f are known to be zero due to being auto vars. */
|
|
Packit |
70b277 |
/* Calculate atan of a known number. */
|
|
Packit |
70b277 |
if (x > .2) {
|
|
Packit |
70b277 |
scale = z+5;
|
|
Packit |
70b277 |
a = a(.2);
|
|
Packit |
70b277 |
}
|
|
Packit |
70b277 |
|
|
Packit |
70b277 |
/* Precondition x. */
|
|
Packit |
70b277 |
scale = z+3;
|
|
Packit |
70b277 |
while (x > .2) {
|
|
Packit |
70b277 |
f += 1;
|
|
Packit |
70b277 |
x = (x-.2) / (1+x*.2);
|
|
Packit |
70b277 |
}
|
|
Packit |
70b277 |
|
|
Packit |
70b277 |
/* Initialize the series. */
|
|
Packit |
70b277 |
v = n = x;
|
|
Packit |
70b277 |
s = -x*x;
|
|
Packit |
70b277 |
|
|
Packit |
70b277 |
/* Calculate the series. */
|
|
Packit |
70b277 |
for (i=3; 1; i+=2) {
|
|
Packit |
70b277 |
e = (n *= s) / i;
|
|
Packit |
70b277 |
if (e == 0) {
|
|
Packit |
70b277 |
scale = z;
|
|
Packit |
70b277 |
return ((f*a+v)/m);
|
|
Packit |
70b277 |
}
|
|
Packit |
70b277 |
v += e
|
|
Packit |
70b277 |
}
|
|
Packit |
70b277 |
}
|
|
Packit |
70b277 |
|
|
Packit |
70b277 |
|
|
Packit |
70b277 |
/* Bessel function of integer order. Uses the following:
|
|
Packit |
70b277 |
j(-n,x) = (-1)^n*j(n,x)
|
|
Packit |
70b277 |
j(n,x) = x^n/(2^n*n!) * (1 - x^2/(2^2*1!*(n+1)) + x^4/(2^4*2!*(n+1)*(n+2))
|
|
Packit |
70b277 |
- x^6/(2^6*3!*(n+1)*(n+2)*(n+3)) .... )
|
|
Packit |
70b277 |
*/
|
|
Packit |
70b277 |
define j(n,x) {
|
|
Packit |
70b277 |
auto a, b, d, e, f, i, m, s, v, z
|
|
Packit |
70b277 |
|
|
Packit |
70b277 |
/* Non base 10 ibase? */
|
|
Packit |
70b277 |
if (ibase != A) {
|
|
Packit |
70b277 |
b = ibase;
|
|
Packit |
70b277 |
ibase = A;
|
|
Packit |
70b277 |
v = j(n,x);
|
|
Packit |
70b277 |
ibase = b;
|
|
Packit |
70b277 |
return (v);
|
|
Packit |
70b277 |
}
|
|
Packit |
70b277 |
|
|
Packit |
70b277 |
/* Make n an integer and check for negative n. */
|
|
Packit |
70b277 |
z = scale;
|
|
Packit |
70b277 |
scale = 0;
|
|
Packit |
70b277 |
n = n/1;
|
|
Packit |
70b277 |
if (n<0) {
|
|
Packit |
70b277 |
n = -n;
|
|
Packit |
70b277 |
if (n%2 == 1) m = 1;
|
|
Packit |
70b277 |
}
|
|
Packit |
70b277 |
|
|
Packit |
70b277 |
/* Compute the factor of x^n/(2^n*n!) */
|
|
Packit |
70b277 |
f = 1;
|
|
Packit |
70b277 |
for (i=2; i<=n; i++) f = f*i;
|
|
Packit |
70b277 |
scale = 1.5*z;
|
|
Packit |
70b277 |
f = x^n / 2^n / f;
|
|
Packit |
70b277 |
|
|
Packit |
70b277 |
/* Initialize the loop .*/
|
|
Packit |
70b277 |
v = e = 1;
|
|
Packit |
70b277 |
s = -x*x/4
|
|
Packit |
70b277 |
scale = 1.5*z + length(f) - scale(f);
|
|
Packit |
70b277 |
|
|
Packit |
70b277 |
/* The Loop.... */
|
|
Packit |
70b277 |
for (i=1; 1; i++) {
|
|
Packit |
70b277 |
e = e * s / i / (n+i);
|
|
Packit |
70b277 |
if (e == 0) {
|
|
Packit |
70b277 |
scale = z
|
|
Packit |
70b277 |
if (m) return (-f*v/1);
|
|
Packit |
70b277 |
return (f*v/1);
|
|
Packit |
70b277 |
}
|
|
Packit |
70b277 |
v += e;
|
|
Packit |
70b277 |
}
|
|
Packit |
70b277 |
}
|