1/* libmath.b for GNU bc.  */
2
3/*  This file is part of GNU bc.
4    Copyright (C) 1991, 1992, 1993, 1997 Free Software Foundation, Inc.
5
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.
10
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.
15
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
20      Boston, MA 02111 USA
21
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
28
29*************************************************************************/
30
31
32scale = 20
33
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! + ...
37*/
38
39define e(x) {
40  auto  a, d, e, f, i, m, n, v, z
41
42  /* a - holds x^y of x^y/y! */
43  /* d - holds 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. */
51
52  /* Check the sign of x. */
53  if (x<0) {
54    m = 1
55    x = -x
56  }
57
58  /* Precondition x. */
59  z = scale;
60  n = 6 + z + .44*x;
61  scale = scale(x)+1;
62  while (x > 1) {
63    f += 1;
64    x /= 2;
65    scale += 1;
66  }
67
68  /* Initialize the variables. */
69  scale = n;
70  v = 1+x
71  a = x
72  d = 1
73
74  for (i=2; 1; i++) {
75    e = (a *= x) / (d *= i)
76    if (e == 0) {
77      if (f>0) while (f--)  v = v*v;
78      scale = z
79      if (m) return (1/v);
80      return (v/1);
81    }
82    v += e
83  }
84}
85
86/* Natural log. Uses the fact that ln(x^2) = 2*ln(x)
87    The series used is:
88       ln(x) = 2(a+a^3/3+a^5/5+...) where a=(x-1)/(x+1)
89*/
90
91define l(x) {
92  auto e, f, i, m, n, v, z
93
94  /* return something for the special case. */
95  if (x <= 0) return ((1 - 10^scale)/1)
96
97  /* Precondition x to make .5 < x < 2.0. */
98  z = scale;
99  scale = 6 + scale;
100  f = 2;
101  i=0
102  while (x >= 2) {  /* for large numbers */
103    f *= 2;
104    x = sqrt(x);
105  }
106  while (x <= .5) {  /* for small numbers */
107    f *= 2;
108    x = sqrt(x);
109  }
110
111  /* Set up the loop. */
112  v = n = (x-1)/(x+1)
113  m = n*n
114
115  /* Sum the series. */
116  for (i=3; 1; i+=2) {
117    e = (n *= m) / i
118    if (e == 0) {
119      v = f*v
120      scale = z
121      return (v/1)
122    }
123    v += e
124  }
125}
126
127/* Sin(x)  uses the standard series:
128   sin(x) = x - x^3/3! + x^5/5! - x^7/7! ... */
129
130define s(x) {
131  auto  e, i, m, n, s, v, z
132
133  /* precondition x. */
134  z = scale
135  scale = 1.1*z + 2;
136  v = a(1)
137  if (x < 0) {
138    m = 1;
139    x = -x;
140  }
141  scale = 0
142  n = (x / v + 2 )/4
143  x = x - 4*n*v
144  if (n%2) x = -x
145
146  /* Do the loop. */
147  scale = z + 2;
148  v = e = x
149  s = -x*x
150  for (i=3; 1; i+=2) {
151    e *= s/(i*(i-1))
152    if (e == 0) {
153      scale = z
154      if (m) return (-v/1);
155      return (v/1);
156    }
157    v += e
158  }
159}
160
161/* Cosine : cos(x) = sin(x+pi/2) */
162define c(x) {
163  auto v, z;
164  z = scale;
165  scale = scale*1.2;
166  v = s(x+a(1)*2);
167  scale = z;
168  return (v/1);
169}
170
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 + ...   */
175
176define a(x) {
177  auto a, e, f, i, m, n, s, v, z
178
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. */
186  /* s is -x*x. */
187  /* z is the saved user's scale. */
188
189  /* Negative x? */
190  m = 1;
191  if (x<0) {
192    m = -1;
193    x = -x;
194  }
195
196  /* Special case and for fast answers */
197  if (x==1) {
198    if (scale <= 25) return (.7853981633974483096156608/m)
199    if (scale <= 40) return (.7853981633974483096156608458198757210492/m)
200    if (scale <= 60) \
201      return (.785398163397448309615660845819875721049292349843776455243736/m)
202  }
203  if (x==.2) {
204    if (scale <= 25) return (.1973955598498807583700497/m)
205    if (scale <= 40) return (.1973955598498807583700497651947902934475/m)
206    if (scale <= 60) \
207      return (.197395559849880758370049765194790293447585103787852101517688/m)
208  }
209
210
211  /* Save the scale. */
212  z = scale;
213
214  /* Note: a and f are known to be zero due to being auto vars. */
215  /* Calculate atan of a known number. */
216  if (x > .2)  {
217    scale = z+5;
218    a = a(.2);
219  }
220
221  /* Precondition x. */
222  scale = z+3;
223  while (x > .2) {
224    f += 1;
225    x = (x-.2) / (1+x*.2);
226  }
227
228  /* Initialize the series. */
229  v = n = x;
230  s = -x*x;
231
232  /* Calculate the series. */
233  for (i=3; 1; i+=2) {
234    e = (n *= s) / i;
235    if (e == 0) {
236      scale = z;
237      return ((f*a+v)/m);
238    }
239    v += e
240  }
241}
242
243
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)) .... )
248*/
249define j(n,x) {
250  auto a, b, d, e, f, i, m, s, v, z
251
252  /* Make n an integer and check for negative n. */
253  z = scale;
254  scale = 0;
255  n = n/1;
256  if (n<0) {
257    n = -n;
258    if (n%2 == 1) m = 1;
259  }
260
261  /* save ibase */
262  b = ibase;
263  ibase = A;
264
265  /* Compute the factor of x^n/(2^n*n!) */
266  f = 1;
267  for (i=2; i<=n; i++) f = f*i;
268  scale = 1.5*z;
269  f = x^n / 2^n / f;
270
271  /* Initialize the loop .*/
272  v = e = 1;
273  s = -x*x/4
274  scale = 1.5*z + length(f) - scale(f);
275
276  /* The Loop.... */
277  for (i=1; 1; i++) {
278    e =  e * s / i / (n+i);
279    if (e == 0) {
280       ibase = b;
281       scale = z
282       if (m) return (-f*v/1);
283       return (f*v/1);
284    }
285    v += e;
286  }
287}
288