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