contrib/bc/bc/libmath.b

   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
  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, 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
  91 define 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
 130 define 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) */
 162 define 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
 176 define 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 */
 249 define 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 }