1/* $NetBSD: libmath.b,v 1.1 2017/04/10 02:28:23 phil Exp $ */ 2 3/* 4 * Copyright (C) 1991-1994, 1997, 2006, 2008, 2012-2017 Free Software Foundation, Inc. 5 * Copyright (C) 2016-2017 Philip A. Nelson. 6 * All rights reserved. 7 * 8 * Redistribution and use in source and binary forms, with or without 9 * modification, are permitted provided that the following conditions 10 * are met: 11 * 12 * 1. Redistributions of source code must retain the above copyright 13 * notice, this list of conditions and the following disclaimer. 14 * 2. Redistributions in binary form must reproduce the above copyright 15 * notice, this list of conditions and the following disclaimer in the 16 * documentation and/or other materials provided with the distribution. 17 * 3. The names Philip A. Nelson and Free Software Foundation may not be 18 * used to endorse or promote products derived from this software 19 * without specific prior written permission. 20 * 21 * THIS SOFTWARE IS PROVIDED BY PHILIP A. NELSON ``AS IS'' AND ANY EXPRESS OR 22 * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES 23 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. 24 * IN NO EVENT SHALL PHILIP A. NELSON OR THE FREE SOFTWARE FOUNDATION BE 25 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR 26 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF 27 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS 28 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN 29 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) 30 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF 31 * THE POSSIBILITY OF SUCH DAMAGE. 32 */ 33 34/* libmath.b for bc. */ 35 36scale = 20 37 38/* Uses the fact that e^x = (e^(x/2))^2 39 When x is small enough, we use the series: 40 e^x = 1 + x + x^2/2! + x^3/3! + ... 41*/ 42 43define e(x) { 44 auto a, b, d, e, f, i, m, n, v, z 45 46 /* a - holds x^y of x^y/y! */ 47 /* d - holds y! */ 48 /* e - is the value x^y/y! */ 49 /* v - is the sum of the e's */ 50 /* f - number of times x was divided by 2. */ 51 /* m - is 1 if x was minus. */ 52 /* i - iteration count. */ 53 /* n - the scale to compute the sum. */ 54 /* z - orignal scale. */ 55 /* b - holds the original ibase. */ 56 57 /* Non base 10 ibase? */ 58 if (ibase != A) { 59 b = ibase; 60 ibase = A; 61 v = e(x); 62 ibase = b; 63 return (v); 64 } 65 66 /* Check the sign of x. */ 67 if (x<0) { 68 m = 1 69 x = -x 70 } 71 72 /* Precondition x. */ 73 z = scale; 74 n = 6 + z + .44*x; 75 scale = scale(x)+1; 76 while (x > 1) { 77 f += 1; 78 x /= 2; 79 scale += 1; 80 } 81 82 /* Initialize the variables. */ 83 scale = n; 84 v = 1+x 85 a = x 86 d = 1 87 88 for (i=2; 1; i++) { 89 e = (a *= x) / (d *= i) 90 if (e == 0) { 91 if (f>0) while (f--) v = v*v; 92 scale = z 93 if (m) return (1/v); 94 return (v/1); 95 } 96 v += e 97 } 98} 99 100/* Natural log. Uses the fact that ln(x^2) = 2*ln(x) 101 The series used is: 102 ln(x) = 2(a+a^3/3+a^5/5+...) where a=(x-1)/(x+1) 103*/ 104 105define l(x) { 106 auto b, e, f, i, m, n, v, z 107 108 /* Non base 10 ibase? */ 109 if (ibase != A) { 110 b = ibase; 111 ibase = A; 112 v = l(x); 113 ibase = b; 114 return (v); 115 } 116 117 /* return something for the special case. */ 118 if (x <= 0) return ((1 - 10^scale)/1) 119 120 /* Precondition x to make .5 < x < 2.0. */ 121 z = scale; 122 scale = 6 + scale; 123 f = 2; 124 i=0 125 while (x >= 2) { /* for large numbers */ 126 f *= 2; 127 x = sqrt(x); 128 } 129 while (x <= .5) { /* for small numbers */ 130 f *= 2; 131 x = sqrt(x); 132 } 133 134 /* Set up the loop. */ 135 v = n = (x-1)/(x+1) 136 m = n*n 137 138 /* Sum the series. */ 139 for (i=3; 1; i+=2) { 140 e = (n *= m) / i 141 if (e == 0) { 142 v = f*v 143 scale = z 144 return (v/1) 145 } 146 v += e 147 } 148} 149 150/* Sin(x) uses the standard series: 151 sin(x) = x - x^3/3! + x^5/5! - x^7/7! ... */ 152 153define s(x) { 154 auto b, e, i, m, n, s, v, z 155 156 /* Non base 10 ibase? */ 157 if (ibase != A) { 158 b = ibase; 159 ibase = A; 160 v = s(x); 161 ibase = b; 162 return (v); 163 } 164 165 /* precondition x. */ 166 z = scale 167 scale = 1.1*z + 2; 168 v = a(1) 169 if (x < 0) { 170 m = 1; 171 x = -x; 172 } 173 scale = 0 174 n = (x / v + 2 )/4 175 x = x - 4*n*v 176 if (n%2) x = -x 177 178 /* Do the loop. */ 179 scale = z + 2; 180 v = e = x 181 s = -x*x 182 for (i=3; 1; i+=2) { 183 e *= s/(i*(i-1)) 184 if (e == 0) { 185 scale = z 186 if (m) return (-v/1); 187 return (v/1); 188 } 189 v += e 190 } 191} 192 193/* Cosine : cos(x) = sin(x+pi/2) */ 194define c(x) { 195 auto b, v, z; 196 197 /* Non base 10 ibase? */ 198 if (ibase != A) { 199 b = ibase; 200 ibase = A; 201 v = c(x); 202 ibase = b; 203 return (v); 204 } 205 206 z = scale; 207 scale = scale*1.2; 208 v = s(x+a(1)*2); 209 scale = z; 210 return (v/1); 211} 212 213/* Arctan: Using the formula: 214 atan(x) = atan(c) + atan((x-c)/(1+xc)) for a small c (.2 here) 215 For under .2, use the series: 216 atan(x) = x - x^3/3 + x^5/5 - x^7/7 + ... */ 217 218define a(x) { 219 auto a, b, e, f, i, m, n, s, v, z 220 221 /* a is the value of a(.2) if it is needed. */ 222 /* f is the value to multiply by a in the return. */ 223 /* e is the value of the current term in the series. */ 224 /* v is the accumulated value of the series. */ 225 /* m is 1 or -1 depending on x (-x -> -1). results are divided by m. */ 226 /* i is the denominator value for series element. */ 227 /* n is the numerator value for the series element. */ 228 /* s is -x*x. */ 229 /* z is the saved user's scale. */ 230 231 /* Non base 10 ibase? */ 232 if (ibase != A) { 233 b = ibase; 234 ibase = A; 235 v = a(x); 236 ibase = b; 237 return (v); 238 } 239 240 /* Negative x? */ 241 m = 1; 242 if (x<0) { 243 m = -1; 244 x = -x; 245 } 246 247 /* Special case and for fast answers */ 248 if (x==1) { 249 if (scale <= 25) return (.7853981633974483096156608/m) 250 if (scale <= 40) return (.7853981633974483096156608458198757210492/m) 251 if (scale <= 60) \ 252 return (.785398163397448309615660845819875721049292349843776455243736/m) 253 } 254 if (x==.2) { 255 if (scale <= 25) return (.1973955598498807583700497/m) 256 if (scale <= 40) return (.1973955598498807583700497651947902934475/m) 257 if (scale <= 60) \ 258 return (.197395559849880758370049765194790293447585103787852101517688/m) 259 } 260 261 262 /* Save the scale. */ 263 z = scale; 264 265 /* Note: a and f are known to be zero due to being auto vars. */ 266 /* Calculate atan of a known number. */ 267 if (x > .2) { 268 scale = z+5; 269 a = a(.2); 270 } 271 272 /* Precondition x. */ 273 scale = z+3; 274 while (x > .2) { 275 f += 1; 276 x = (x-.2) / (1+x*.2); 277 } 278 279 /* Initialize the series. */ 280 v = n = x; 281 s = -x*x; 282 283 /* Calculate the series. */ 284 for (i=3; 1; i+=2) { 285 e = (n *= s) / i; 286 if (e == 0) { 287 scale = z; 288 return ((f*a+v)/m); 289 } 290 v += e 291 } 292} 293 294 295/* Bessel function of integer order. Uses the following: 296 j(-n,x) = (-1)^n*j(n,x) 297 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)) 298 - x^6/(2^6*3!*(n+1)*(n+2)*(n+3)) .... ) 299*/ 300define j(n,x) { 301 auto a, b, d, e, f, i, m, s, v, z 302 303 /* Non base 10 ibase? */ 304 if (ibase != A) { 305 b = ibase; 306 ibase = A; 307 v = j(n,x); 308 ibase = b; 309 return (v); 310 } 311 312 /* Make n an integer and check for negative n. */ 313 z = scale; 314 scale = 0; 315 n = n/1; 316 if (n<0) { 317 n = -n; 318 if (n%2 == 1) m = 1; 319 } 320 321 /* Compute the factor of x^n/(2^n*n!) */ 322 f = 1; 323 for (i=2; i<=n; i++) f = f*i; 324 scale = 1.5*z; 325 f = x^n / 2^n / f; 326 327 /* Initialize the loop .*/ 328 v = e = 1; 329 s = -x*x/4 330 scale = 1.5*z + length(f) - scale(f); 331 332 /* The Loop.... */ 333 for (i=1; 1; i++) { 334 e = e * s / i / (n+i); 335 if (e == 0) { 336 scale = z 337 if (m) return (-f*v/1); 338 return (f*v/1); 339 } 340 v += e; 341 } 342} 343