Merge from vendor branch DIFFUTILS:
[dragonfly.git] / contrib / bc / bc / libmath.b
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
32 scale = 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
39 define 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
91 define 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
130 define 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) */
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 }
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
176 define 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 */
249 define 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 }