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
20       Boston, MA 02111 USA
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 *************************************************************************/
32 scale = 20
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 */
39 define e(x) {
40   auto  a, d, e, f, i, m, n, v, z
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. */
52   /* Check the sign of x. */
53   if (x<0) {
54     m = 1
55     x = -x
56   }
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   }
68   /* Initialize the variables. */
69   scale = n;
70   v = 1+x
71   a = x
72   d = 1
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 }
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 */
91 define l(x) {
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. */
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   }
111   /* Set up the loop. */
112   v = n = (x-1)/(x+1)
113   m = n*n
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 }
127 /* Sin(x)  uses the standard series:
128    sin(x) = x - x^3/3! + x^5/5! - x^7/7! ... */
130 define s(x) {
131   auto  e, i, m, n, s, v, z
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
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 }
161 /* Cosine : cos(x) = sin(x+pi/2) */
162 define 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 }
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 + ...   */
176 define a(x) {
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. */
186   /* s is -x*x. */
187   /* z is the saved user's scale. */
189   /* Negative x? */
190   m = 1;
191   if (x<0) {
192     m = -1;
193     x = -x;
194   }
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   }
211   /* Save the scale. */
212   z = scale;
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   }
221   /* Precondition x. */
222   scale = z+3;
223   while (x > .2) {
224     f += 1;
225     x = (x-.2) / (1+x*.2);
226   }
228   /* Initialize the series. */
229   v = n = x;
230   s = -x*x;
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 }
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 */
249 define j(n,x) {
250   auto a, b, d, e, f, i, m, s, v, z
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   }
261   /* save ibase */
262   b = ibase;
263   ibase = A;
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;
271   /* Initialize the loop .*/
272   v = e = 1;
273   s = -x*x/4
274   scale = 1.5*z + length(f) - scale(f);
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 }