1 /* libmath.b for GNU bc. */
3 /* This file is part of GNU bc.
4 Copyright (C) 1991, 1992, 1993, 1997 Free Software Foundation, Inc.
6 This program is free software; you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation; either version 2 of the License , or
9 (at your option) any later version.
11 This program is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
16 You should have received a copy of the GNU General Public License
17 along with this program; see the file COPYING. If not, write to
18 The Free Software Foundation, Inc.
19 59 Temple Place, Suite 330
22 You may contact the author by:
23 e-mail: philnelson@acm.org
24 us-mail: Philip A. Nelson
25 Computer Science Department, 9062
26 Western Washington University
27 Bellingham, WA 98226-9062
29 *************************************************************************/
34 /* Uses the fact that e^x = (e^(x/2))^2
35 When x is small enough, we use the series:
36 e^x = 1 + x + x^2/2! + x^3/3! + ...
40 auto a, d, e, f, i, m, n, v, z
42 /* a - holds x^y of x^y/y! */
44 /* e - is the value x^y/y! */
45 /* v - is the sum of the e's */
46 /* f - number of times x was divided by 2. */
47 /* m - is 1 if x was minus. */
48 /* i - iteration count. */
49 /* n - the scale to compute the sum. */
50 /* z - orignal scale. */
52 /* Check the sign of x. */
68 /* Initialize the variables. */
75 e = (a *= x) / (d *= i)
77 if (f>0) while (f--) v = v*v;
86 /* Natural log. Uses the fact that ln(x^2) = 2*ln(x)
88 ln(x) = 2(a+a^3/3+a^5/5+...) where a=(x-1)/(x+1)
92 auto e, f, i, m, n, v, z
94 /* return something for the special case. */
95 if (x <= 0) return ((1 - 10^scale)/1)
97 /* Precondition x to make .5 < x < 2.0. */
102 while (x >= 2) { /* for large numbers */
106 while (x <= .5) { /* for small numbers */
111 /* Set up the loop. */
115 /* Sum the series. */
127 /* Sin(x) uses the standard series:
128 sin(x) = x - x^3/3! + x^5/5! - x^7/7! ... */
131 auto e, i, m, n, s, v, z
133 /* precondition x. */
154 if (m) return (-v/1);
161 /* Cosine : cos(x) = sin(x+pi/2) */
171 /* Arctan: Using the formula:
172 atan(x) = atan(c) + atan((x-c)/(1+xc)) for a small c (.2 here)
173 For under .2, use the series:
174 atan(x) = x - x^3/3 + x^5/5 - x^7/7 + ... */
177 auto a, e, f, i, m, n, s, v, z
179 /* a is the value of a(.2) if it is needed. */
180 /* f is the value to multiply by a in the return. */
181 /* e is the value of the current term in the series. */
182 /* v is the accumulated value of the series. */
183 /* m is 1 or -1 depending on x (-x -> -1). results are divided by m. */
184 /* i is the denominator value for series element. */
185 /* n is the numerator value for the series element. */
187 /* z is the saved user's scale. */
196 /* Special case and for fast answers */
198 if (scale <= 25) return (.7853981633974483096156608/m)
199 if (scale <= 40) return (.7853981633974483096156608458198757210492/m)
201 return (.785398163397448309615660845819875721049292349843776455243736/m)
204 if (scale <= 25) return (.1973955598498807583700497/m)
205 if (scale <= 40) return (.1973955598498807583700497651947902934475/m)
207 return (.197395559849880758370049765194790293447585103787852101517688/m)
211 /* Save the scale. */
214 /* Note: a and f are known to be zero due to being auto vars. */
215 /* Calculate atan of a known number. */
221 /* Precondition x. */
225 x = (x-.2) / (1+x*.2);
228 /* Initialize the series. */
232 /* Calculate the series. */
244 /* Bessel function of integer order. Uses the following:
245 j(-n,x) = (-1)^n*j(n,x)
246 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))
247 - x^6/(2^6*3!*(n+1)*(n+2)*(n+3)) .... )
250 auto a, b, d, e, f, i, m, s, v, z
252 /* Make n an integer and check for negative n. */
265 /* Compute the factor of x^n/(2^n*n!) */
267 for (i=2; i<=n; i++) f = f*i;
271 /* Initialize the loop .*/
274 scale = 1.5*z + length(f) - scale(f);
278 e = e * s / i / (n+i);
282 if (m) return (-f*v/1);