i86/elks/elkscmd/bc/libmath.b

   1 /* libmath.b for bc for minix.  */
   2
   3 /*  This file is part of bc written for MINIX.
   4     Copyright (C) 1991, 1992 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 }