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, 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, n, v, z 39: 40: /* a - holds x^y of x^y/y! */ 41: /* d - holds y! */ 42: /* e - is the value x^y/y! */ 43: /* v - is the sum of the e's */ 44: /* f - number of times x was divided by 2. */ 45: /* m - is 1 if x was minus. */ 46: /* i - iteration count. */ 47: /* n - the scale to compute the sum. */ 48: /* z - orignal scale. */ 49: 50: /* Check the sign of x. */ 51: if (x<0) { 52: m = 1 53: x = -x 54: } 55: 56: /* Precondition x. */ 57: z = scale; 58: n = 6 + z + .44*x; 59: scale = scale(x)+1; 60: while (x > 1) { 61: f += 1; 62: x /= 2; 63: scale += 1; 64: } 65: 66: /* Initialize the variables. */ 67: scale = n; 68: v = 1+x 69: a = x 70: d = 1 71: 72: for (i=2; 1; i++) { 73: e = (a *= x) / (d *= i) 74: if (e == 0) { 75: if (f>0) while (f--) v = v*v; 76: scale = z 77: if (m) return (1/v); 78: return (v/1); 79: } 80: v += e 81: } 82: } 83: 84: /* Natural log. Uses the fact that ln(x^2) = 2*ln(x) 85: The series used is: 86: ln(x) = 2(a+a^3/3+a^5/5+...) where a=(x-1)/(x+1) 87: */ 88: 89: define l(x) { 90: auto e, f, i, m, n, v, z 91: 92: /* return something for the special case. */ 93: if (x <= 0) return ((1 - 10^scale)/1) 94: 95: /* Precondition x to make .5 < x < 2.0. */ 96: z = scale; 97: scale = 6 + scale; 98: f = 2; 99: i=0 100: while (x >= 2) { /* for large numbers */ 101: f *= 2; 102: x = sqrt(x); 103: } 104: while (x <= .5) { /* for small numbers */ 105: f *= 2; 106: x = sqrt(x); 107: } 108: 109: /* Set up the loop. */ 110: v = n = (x-1)/(x+1) 111: m = n*n 112: 113: /* Sum the series. */ 114: for (i=3; 1; i+=2) { 115: e = (n *= m) / i 116: if (e == 0) { 117: v = f*v 118: scale = z 119: return (v/1) 120: } 121: v += e 122: } 123: } 124: 125: /* Sin(x) uses the standard series: 126: sin(x) = x - x^3/3! + x^5/5! - x^7/7! ... */ 127: 128: define s(x) { 129: auto e, i, m, n, s, v, z 130: 131: /* precondition x. */ 132: z = scale 133: scale = 1.1*z + 2; 134: v = a(1) 135: if (x < 0) { 136: m = 1; 137: x = -x; 138: } 139: scale = 0 140: n = (x / v + 2 )/4 141: x = x - 4*n*v 142: if (n%2) x = -x 143: 144: /* Do the loop. */ 145: scale = z + 2; 146: v = e = x 147: s = -x*x 148: for (i=3; 1; i+=2) { 149: e *= s/(i*(i-1)) 150: if (e == 0) { 151: scale = z 152: if (m) return (-v/1); 153: return (v/1); 154: } 155: v += e 156: } 157: } 158: 159: /* Cosine : cos(x) = sin(x+pi/2) */ 160: define c(x) { 161: auto v; 162: scale += 1; 163: v = s(x+a(1)*2); 164: scale -= 1; 165: return (v/1); 166: } 167: 168: /* Arctan: Using the formula: 169: atan(x) = atan(c) + atan((x-c)/(1+xc)) for a small c (.2 here) 170: For under .2, use the series: 171: atan(x) = x - x^3/3 + x^5/5 - x^7/7 + ... */ 172: 173: define a(x) { 174: auto a, e, f, i, m, n, s, v, z 175: 176: /* a is the value of a(.2) if it is needed. */ 177: /* f is the value to multiply by a in the return. */ 178: /* e is the value of the current term in the series. */ 179: /* v is the accumulated value of the series. */ 180: /* m is 1 or -1 depending on x (-x -> -1). results are divided by m. */ 181: /* i is the denominator value for series element. */ 182: /* n is the numerator value for the series element. */ 183: /* s is -x*x. */ 184: /* z is the saved user's scale. */ 185: 186: /* Negative x? */ 187: m = 1; 188: if (x<0) { 189: m = -1; 190: x = -x; 191: } 192: 193: /* Special case and for fast answers */ 194: if (x==1) { 195: if (scale <= 25) return (.7853981633974483096156608/m) 196: if (scale <= 40) return (.7853981633974483096156608458198757210492/m) 197: if (scale <= 60) \ 198: return (.785398163397448309615660845819875721049292349843776455243736/m) 199: } 200: if (x==.2) { 201: if (scale <= 25) return (.1973955598498807583700497/m) 202: if (scale <= 40) return (.1973955598498807583700497651947902934475/m) 203: if (scale <= 60) \ 204: return (.197395559849880758370049765194790293447585103787852101517688/m) 205: } 206: 207: 208: /* Save the scale. */ 209: z = scale; 210: 211: /* Note: a and f are known to be zero due to being auto vars. */ 212: /* Calculate atan of a known number. */ 213: if (x > .2) { 214: scale = z+5; 215: a = a(.2); 216: } 217: 218: /* Precondition x. */ 219: scale = z+3; 220: while (x > .2) { 221: f += 1; 222: x = (x-.2) / (1+x*.2); 223: } 224: 225: /* Initialize the series. */ 226: v = n = x; 227: s = -x*x; 228: 229: /* Calculate the series. */ 230: for (i=3; 1; i+=2) { 231: e = (n *= s) / i; 232: if (e == 0) { 233: scale = z; 234: return ((f*a+v)/m); 235: } 236: v += e 237: } 238: } 239: 240: 241: /* Bessel function of integer order. Uses the following: 242: j(-n,x) = (-1)^n*j(n,x) 243: 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)) 244: - x^6/(2^6*3!*(n+1)*(n+2)*(n+3)) .... ) 245: */ 246: define j(n,x) { 247: auto a, d, e, f, i, m, s, v, z 248: 249: /* Make n an integer and check for negative n. */ 250: z = scale; 251: scale = 0; 252: n = n/1; 253: if (n<0) { 254: n = -n; 255: if (n%2 == 1) m = 1; 256: } 257: 258: /* Compute the factor of x^n/(2^n*n!) */ 259: f = 1; 260: for (i=2; i<=n; i++) f = f*i; 261: scale = 1.5*z; 262: f = x^n / 2^n / f; 263: 264: /* Initialize the loop .*/ 265: v = e = 1; 266: s = -x*x/4 267: scale = 1.5*z 268: 269: /* The Loop.... */ 270: for (i=1; 1; i++) { 271: e = e * s / i / (n+i); 272: if (e == 0) { 273: scale = z 274: if (m) return (-f*v/1); 275: return (f*v/1); 276: } 277: v += e; 278: } 279: } |