Commit | Line | Data |
---|---|---|
c4f36ac6 WJ |
1 | /* libmath.b for bc for minix. */ |
2 | ||
3 | /* This file is part of bc written for MINIX. | |
4 | Copyright (C) 1991 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, 675 Mass Ave, Cambridge, MA 02139, USA. | |
19 | ||
20 | You may contact the author by: | |
21 | e-mail: phil@cs.wwu.edu | |
22 | us-mail: Philip A. Nelson | |
23 | Computer Science Department, 9062 | |
24 | Western Washington University | |
25 | Bellingham, WA 98226-9062 | |
26 | ||
27 | *************************************************************************/ | |
28 | ||
29 | ||
30 | scale = 20 | |
31 | ||
32 | /* Uses the fact that e^x = (e^(x/2))^2 | |
33 | When x is small enough, we use the series: | |
34 | e^x = 1 + x + x^2/2! + x^3/3! + ... | |
35 | */ | |
36 | ||
37 | define e(x) { | |
38 | auto a, d, e, f, i, m, v, z | |
39 | ||
40 | /* Check the sign of x. */ | |
41 | if (x<0) { | |
42 | m = 1 | |
43 | x = -x | |
44 | } | |
45 | ||
46 | /* Precondition x. */ | |
47 | z = scale; | |
48 | scale = 4 + z + .44*x; | |
49 | while (x > 1) { | |
50 | f += 1; | |
51 | x /= 2; | |
52 | } | |
53 | ||
54 | /* Initialize the variables. */ | |
55 | v = 1+x | |
56 | a = x | |
57 | d = 1 | |
58 | ||
59 | for (i=2; 1; i++) { | |
60 | e = (a *= x) / (d *= i) | |
61 | if (e == 0) { | |
62 | if (f>0) while (f--) v = v*v; | |
63 | scale = z | |
64 | if (m) return (1/v); | |
65 | return (v/1); | |
66 | } | |
67 | v += e | |
68 | } | |
69 | } | |
70 | ||
71 | /* Natural log. Uses the fact that ln(x^2) = 2*ln(x) | |
72 | The series used is: | |
73 | ln(x) = 2(a+a^3/3+a^5/5+...) where a=(x-1)/(x+1) | |
74 | */ | |
75 | ||
76 | define l(x) { | |
77 | auto e, f, i, m, n, v, z | |
78 | ||
79 | /* return something for the special case. */ | |
80 | if (x <= 0) return (1 - 10^scale) | |
81 | ||
82 | /* Precondition x to make .5 < x < 2.0. */ | |
83 | z = scale; | |
84 | scale += 4; | |
85 | f = 2; | |
86 | i=0 | |
87 | while (x >= 2) { /* for large numbers */ | |
88 | f *= 2; | |
89 | x = sqrt(x); | |
90 | } | |
91 | while (x <= .5) { /* for small numbers */ | |
92 | f *= 2; | |
93 | x = sqrt(x); | |
94 | } | |
95 | ||
96 | /* Set up the loop. */ | |
97 | v = n = (x-1)/(x+1) | |
98 | m = n*n | |
99 | ||
100 | /* Sum the series. */ | |
101 | for (i=3; 1; i+=2) { | |
102 | e = (n *= m) / i | |
103 | if (e == 0) { | |
104 | v = f*v | |
105 | scale = z | |
106 | return (v/1) | |
107 | } | |
108 | v += e | |
109 | } | |
110 | } | |
111 | ||
112 | /* Sin(x) uses the standard series: | |
113 | sin(x) = x - x^3/3! + x^5/5! - x^7/7! ... */ | |
114 | ||
115 | define s(x) { | |
116 | auto e, i, m, n, s, v, z | |
117 | ||
118 | /* precondition x. */ | |
119 | z = scale | |
120 | scale = 1.1*z + 1; | |
121 | v = a(1) | |
122 | if (x < 0) { | |
123 | m = 1; | |
124 | x = -x; | |
125 | } | |
126 | scale = 0 | |
127 | n = (x / v + 2 )/4 | |
128 | x = x - 4*n*v | |
129 | if (n%2) x = -x | |
130 | ||
131 | /* Do the loop. */ | |
132 | scale = z + 2; | |
133 | v = e = x | |
134 | s = -x*x | |
135 | for (i=3; 1; i+=2) { | |
136 | e *= s/(i*(i-1)) | |
137 | if (e == 0) { | |
138 | scale = z | |
139 | if (m) return (-v/1); | |
140 | return (v/1); | |
141 | } | |
142 | v += e | |
143 | } | |
144 | } | |
145 | ||
146 | /* Cosine : cos(x) = sin(x+pi/2) */ | |
147 | define c(x) { | |
148 | auto v; | |
149 | scale += 1; | |
150 | v = s(x+a(1)*2); | |
151 | scale -= 1; | |
152 | return (v/1); | |
153 | } | |
154 | ||
155 | /* Arctan: Using the formula: | |
156 | atan(x) = atan(c) + atan((x-c)/(1+xc)) for a small c (.2 here) | |
157 | For under .2, use the series: | |
158 | atan(x) = x - x^3/3 + x^5/5 - x^7/7 + ... */ | |
159 | ||
160 | define a(x) { | |
161 | auto a, e, f, i, m, n, s, v, z | |
162 | ||
163 | /* Special case and for fast answers */ | |
164 | if (x==1) { | |
165 | if (scale <= 25) return (.7853981633974483096156608/1) | |
166 | if (scale <= 40) return (.7853981633974483096156608458198757210492/1) | |
167 | if (scale <= 60) \ | |
168 | return (.785398163397448309615660845819875721049292349843776455243736/1) | |
169 | } | |
170 | if (x==.2) { | |
171 | if (scale <= 25) return (.1973955598498807583700497/1) | |
172 | if (scale <= 40) return (.1973955598498807583700497651947902934475/1) | |
173 | if (scale <= 60) \ | |
174 | return (.197395559849880758370049765194790293447585103787852101517688/1) | |
175 | } | |
176 | ||
177 | /* Negative x? */ | |
178 | if (x<0) { | |
179 | m = 1; | |
180 | x = -x; | |
181 | } | |
182 | ||
183 | /* Save the scale. */ | |
184 | z = scale; | |
185 | ||
186 | /* Note: a and f are known to be zero due to being auto vars. */ | |
187 | /* Calculate atan of a known number. */ | |
188 | if (x > .2) { | |
189 | scale = z+4; | |
190 | a = a(.2); | |
191 | } | |
192 | ||
193 | /* Precondition x. */ | |
194 | scale = z+2; | |
195 | while (x > .2) { | |
196 | f += 1; | |
197 | x = (x-.2) / (1+x*.2); | |
198 | } | |
199 | ||
200 | /* Initialize the series. */ | |
201 | v = n = x; | |
202 | s = -x*x; | |
203 | ||
204 | /* Calculate the series. */ | |
205 | for (i=3; 1; i+=2) { | |
206 | e = (n *= s) / i; | |
207 | if (e == 0) { | |
208 | scale = z; | |
209 | if (m) return ((f*a+v)/-1); | |
210 | return ((f*a+v)/1); | |
211 | } | |
212 | v += e | |
213 | } | |
214 | } | |
215 | ||
216 | ||
217 | /* Bessel function of integer order. Uses the following: | |
218 | j(-n,x) = (-1)^n*j(n,x) | |
219 | 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)) | |
220 | - x^6/(2^6*3!*(n+1)*(n+2)*(n+3)) .... ) | |
221 | */ | |
222 | define j(n,x) { | |
223 | auto a, d, e, f, i, m, s, v, z | |
224 | ||
225 | /* Make n an integer and check for negative n. */ | |
226 | z = scale; | |
227 | scale = 0; | |
228 | n = n/1; | |
229 | if (n<0) { | |
230 | n = -n; | |
231 | if (n%2 == 1) m = 1; | |
232 | } | |
233 | ||
234 | /* Compute the factor of x^n/(2^n*n!) */ | |
235 | f = 1; | |
236 | for (i=2; i<=n; i++) f = f*i; | |
237 | scale = 1.5*z; | |
238 | f = x^n / 2^n / f; | |
239 | ||
240 | /* Initialize the loop .*/ | |
241 | v = e = 1; | |
242 | s = -x*x/4 | |
243 | scale = 1.5*z | |
244 | ||
245 | /* The Loop.... */ | |
246 | for (i=1; 1; i++) { | |
247 | e = e * s / i / (n+i); | |
248 | if (e == 0) { | |
249 | scale = z | |
250 | if (m) return (-f*v/1); | |
251 | return (f*v/1); | |
252 | } | |
253 | v += e; | |
254 | } | |
255 | } |