home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / bc-1.03-base.tgz / bc-1.03-base.tar / fsf / bc / libmath.b < prev    next >
Text File  |  1994-08-07  |  6KB  |  280 lines
  1. /* libmath.b for bc for minix.  */
  2.  
  3. /*  This file is part of bc written for MINIX.
  4.     Copyright (C) 1991, 1992, 1993 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)
  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. }
  280.