## "Fossies" - the Fresh Open Source Software Archive ### Member "bc-1.06.95/bc/libmath.b" (5 Sep 2006, 7110 Bytes) of package /linux/misc/old/bc-1.06.95.tar.gz:

As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) Limbo source code syntax highlighting (style: standard) with prefixed line numbers. Alternatively you can here view or download the uninterpreted source code file. See also the last Fossies "Diffs" side-by-side code changes report for "libmath.b": 1.06.95_vs_1.07.

```    1 /*  This file is part of GNU bc.
2
3     Copyright (C) 1991-1994, 1997, 2006 Free Software Foundation, Inc.
4
5     This program is free software; you can redistribute it and/or modify
7     the Free Software Foundation; either version 2 of the License , or
8     (at your option) any later version.
9
10     This program is distributed in the hope that it will be useful,
11     but WITHOUT ANY WARRANTY; without even the implied warranty of
12     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13     GNU General Public License for more details.
14
15     You should have received a copy of the GNU General Public License
16     along with this program; see the file COPYING.  If not, write to:
17       The Free Software Foundation, Inc.
18       Foundation, Inc.  51 Franklin Street, Fifth Floor,
19       Boston, MA 02110-1301  USA
20
21     You may contact the author by:
22        e-mail:  philnelson@acm.org
23       us-mail:  Philip A. Nelson
24                 Computer Science Department, 9062
25                 Western Washington University
26                 Bellingham, WA 98226-9062
27
28 *************************************************************************/
29
30 /* libmath.b for bc.  */
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, b, 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   /* b - holds the original ibase. */
52
53   /* Non base 10 ibase? */
54   if (ibase != A) {
55      b = ibase;
56      ibase = A;
57      v = e(x);
58      ibase = b;
59      return (v);
60   }
61
62   /* Check the sign of x. */
63   if (x<0) {
64     m = 1
65     x = -x
66   }
67
68   /* Precondition x. */
69   z = scale;
70   n = 6 + z + .44*x;
71   scale = scale(x)+1;
72   while (x > 1) {
73     f += 1;
74     x /= 2;
75     scale += 1;
76   }
77
78   /* Initialize the variables. */
79   scale = n;
80   v = 1+x
81   a = x
82   d = 1
83
84   for (i=2; 1; i++) {
85     e = (a *= x) / (d *= i)
86     if (e == 0) {
87       if (f>0) while (f--)  v = v*v;
88       scale = z
89       if (m) return (1/v);
90       return (v/1);
91     }
92     v += e
93   }
94 }
95
96 /* Natural log. Uses the fact that ln(x^2) = 2*ln(x)
97     The series used is:
98        ln(x) = 2(a+a^3/3+a^5/5+...) where a=(x-1)/(x+1)
99 */
100
101 define l(x) {
102   auto b, e, f, i, m, n, v, z
103
104   /* Non base 10 ibase? */
105   if (ibase != A) {
106      b = ibase;
107      ibase = A;
108      v = l(x);
109      ibase = b;
110      return (v);
111   }
112
113   /* return something for the special case. */
114   if (x <= 0) return ((1 - 10^scale)/1)
115
116   /* Precondition x to make .5 < x < 2.0. */
117   z = scale;
118   scale = 6 + scale;
119   f = 2;
120   i=0
121   while (x >= 2) {  /* for large numbers */
122     f *= 2;
123     x = sqrt(x);
124   }
125   while (x <= .5) {  /* for small numbers */
126     f *= 2;
127     x = sqrt(x);
128   }
129
130   /* Set up the loop. */
131   v = n = (x-1)/(x+1)
132   m = n*n
133
134   /* Sum the series. */
135   for (i=3; 1; i+=2) {
136     e = (n *= m) / i
137     if (e == 0) {
138       v = f*v
139       scale = z
140       return (v/1)
141     }
142     v += e
143   }
144 }
145
146 /* Sin(x)  uses the standard series:
147    sin(x) = x - x^3/3! + x^5/5! - x^7/7! ... */
148
149 define s(x) {
150   auto  b, e, i, m, n, s, v, z
151
152   /* Non base 10 ibase? */
153   if (ibase != A) {
154      b = ibase;
155      ibase = A;
156      v = s(x);
157      ibase = b;
158      return (v);
159   }
160
161   /* precondition x. */
162   z = scale
163   scale = 1.1*z + 2;
164   v = a(1)
165   if (x < 0) {
166     m = 1;
167     x = -x;
168   }
169   scale = 0
170   n = (x / v + 2 )/4
171   x = x - 4*n*v
172   if (n%2) x = -x
173
174   /* Do the loop. */
175   scale = z + 2;
176   v = e = x
177   s = -x*x
178   for (i=3; 1; i+=2) {
179     e *= s/(i*(i-1))
180     if (e == 0) {
181       scale = z
182       if (m) return (-v/1);
183       return (v/1);
184     }
185     v += e
186   }
187 }
188
189 /* Cosine : cos(x) = sin(x+pi/2) */
190 define c(x) {
191   auto b, v, z;
192
193   /* Non base 10 ibase? */
194   if (ibase != A) {
195      b = ibase;
196      ibase = A;
197      v = c(x);
198      ibase = b;
199      return (v);
200   }
201
202   z = scale;
203   scale = scale*1.2;
204   v = s(x+a(1)*2);
205   scale = z;
206   return (v/1);
207 }
208
209 /* Arctan: Using the formula:
210      atan(x) = atan(c) + atan((x-c)/(1+xc)) for a small c (.2 here)
211    For under .2, use the series:
212      atan(x) = x - x^3/3 + x^5/5 - x^7/7 + ...   */
213
214 define a(x) {
215   auto a, b, e, f, i, m, n, s, v, z
216
217   /* a is the value of a(.2) if it is needed. */
218   /* f is the value to multiply by a in the return. */
219   /* e is the value of the current term in the series. */
220   /* v is the accumulated value of the series. */
221   /* m is 1 or -1 depending on x (-x -> -1).  results are divided by m. */
222   /* i is the denominator value for series element. */
223   /* n is the numerator value for the series element. */
224   /* s is -x*x. */
225   /* z is the saved user's scale. */
226
227   /* Non base 10 ibase? */
228   if (ibase != A) {
229      b = ibase;
230      ibase = A;
231      v = a(x);
232      ibase = b;
233      return (v);
234   }
235
236   /* Negative x? */
237   m = 1;
238   if (x<0) {
239     m = -1;
240     x = -x;
241   }
242
243   /* Special case and for fast answers */
244   if (x==1) {
245     if (scale <= 25) return (.7853981633974483096156608/m)
246     if (scale <= 40) return (.7853981633974483096156608458198757210492/m)
247     if (scale <= 60) \
248       return (.785398163397448309615660845819875721049292349843776455243736/m)
249   }
250   if (x==.2) {
251     if (scale <= 25) return (.1973955598498807583700497/m)
252     if (scale <= 40) return (.1973955598498807583700497651947902934475/m)
253     if (scale <= 60) \
254       return (.197395559849880758370049765194790293447585103787852101517688/m)
255   }
256
257
258   /* Save the scale. */
259   z = scale;
260
261   /* Note: a and f are known to be zero due to being auto vars. */
262   /* Calculate atan of a known number. */
263   if (x > .2)  {
264     scale = z+5;
265     a = a(.2);
266   }
267
268   /* Precondition x. */
269   scale = z+3;
270   while (x > .2) {
271     f += 1;
272     x = (x-.2) / (1+x*.2);
273   }
274
275   /* Initialize the series. */
276   v = n = x;
277   s = -x*x;
278
279   /* Calculate the series. */
280   for (i=3; 1; i+=2) {
281     e = (n *= s) / i;
282     if (e == 0) {
283       scale = z;
284       return ((f*a+v)/m);
285     }
286     v += e
287   }
288 }
289
290
291 /* Bessel function of integer order.  Uses the following:
292    j(-n,x) = (-1)^n*j(n,x)
293    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))
294             - x^6/(2^6*3!*(n+1)*(n+2)*(n+3)) .... )
295 */
296 define j(n,x) {
297   auto a, b, d, e, f, i, m, s, v, z
298
299   /* Non base 10 ibase? */
300   if (ibase != A) {
301      b = ibase;
302      ibase = A;
303      v = j(n,x);
304      ibase = b;
305      return (v);
306   }
307
308   /* Make n an integer and check for negative n. */
309   z = scale;
310   scale = 0;
311   n = n/1;
312   if (n<0) {
313     n = -n;
314     if (n%2 == 1) m = 1;
315   }
316
317   /* Compute the factor of x^n/(2^n*n!) */
318   f = 1;
319   for (i=2; i<=n; i++) f = f*i;
320   scale = 1.5*z;
321   f = x^n / 2^n / f;
322
323   /* Initialize the loop .*/
324   v = e = 1;
325   s = -x*x/4
326   scale = 1.5*z + length(f) - scale(f);
327
328   /* The Loop.... */
329   for (i=1; 1; i++) {
330     e =  e * s / i / (n+i);
331     if (e == 0) {
332        scale = z
333        if (m) return (-f*v/1);
334        return (f*v/1);
335     }
336     v += e;
337   }
338 }
```