"Fossies" - the Fresh Open Source Software Archive

Member "bc-1.06.95/bc/libmath.b" (5 Sep 2006, 7110 Bytes) of package /linux/misc/old/bc-1.06.95.tar.gz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) Limbo source code syntax highlighting (style: standard) with prefixed line numbers. Alternatively you can here view or download the uninterpreted source code file. See also the last Fossies "Diffs" side-by-side code changes report for "libmath.b": 1.06.95_vs_1.07.

    1 /*  This file is part of GNU bc.
    2 
    3     Copyright (C) 1991-1994, 1997, 2006 Free Software Foundation, Inc.
    4 
    5     This program is free software; you can redistribute it and/or modify
    6     it under the terms of the GNU General Public License as published by
    7     the Free Software Foundation; either version 2 of the License , or
    8     (at your option) any later version.
    9 
   10     This program is distributed in the hope that it will be useful,
   11     but WITHOUT ANY WARRANTY; without even the implied warranty of
   12     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   13     GNU General Public License for more details.
   14 
   15     You should have received a copy of the GNU General Public License
   16     along with this program; see the file COPYING.  If not, write to:
   17       The Free Software Foundation, Inc.
   18       Foundation, Inc.  51 Franklin Street, Fifth Floor,
   19       Boston, MA 02110-1301  USA
   20 
   21     You may contact the author by:
   22        e-mail:  philnelson@acm.org
   23       us-mail:  Philip A. Nelson
   24                 Computer Science Department, 9062
   25                 Western Washington University
   26                 Bellingham, WA 98226-9062
   27        
   28 *************************************************************************/
   29 
   30 /* libmath.b for bc.  */
   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, b, 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   /* b - holds the original ibase. */
   52 
   53   /* Non base 10 ibase? */
   54   if (ibase != A) {
   55      b = ibase;
   56      ibase = A;
   57      v = e(x);
   58      ibase = b;
   59      return (v);
   60   }
   61 
   62   /* Check the sign of x. */
   63   if (x<0) {
   64     m = 1
   65     x = -x
   66   } 
   67 
   68   /* Precondition x. */
   69   z = scale;
   70   n = 6 + z + .44*x;
   71   scale = scale(x)+1;
   72   while (x > 1) {
   73     f += 1;
   74     x /= 2;
   75     scale += 1;
   76   }
   77 
   78   /* Initialize the variables. */
   79   scale = n;
   80   v = 1+x
   81   a = x
   82   d = 1
   83 
   84   for (i=2; 1; i++) {
   85     e = (a *= x) / (d *= i)
   86     if (e == 0) {
   87       if (f>0) while (f--)  v = v*v;
   88       scale = z
   89       if (m) return (1/v);
   90       return (v/1);
   91     }
   92     v += e
   93   }
   94 }
   95 
   96 /* Natural log. Uses the fact that ln(x^2) = 2*ln(x)
   97     The series used is:
   98        ln(x) = 2(a+a^3/3+a^5/5+...) where a=(x-1)/(x+1)
   99 */
  100 
  101 define l(x) {
  102   auto b, e, f, i, m, n, v, z
  103 
  104   /* Non base 10 ibase? */
  105   if (ibase != A) {
  106      b = ibase;
  107      ibase = A;
  108      v = l(x);
  109      ibase = b;
  110      return (v);
  111   }
  112 
  113   /* return something for the special case. */
  114   if (x <= 0) return ((1 - 10^scale)/1)
  115 
  116   /* Precondition x to make .5 < x < 2.0. */
  117   z = scale;
  118   scale = 6 + scale;
  119   f = 2;
  120   i=0
  121   while (x >= 2) {  /* for large numbers */
  122     f *= 2;
  123     x = sqrt(x);
  124   }
  125   while (x <= .5) {  /* for small numbers */
  126     f *= 2;
  127     x = sqrt(x);
  128   }
  129 
  130   /* Set up the loop. */
  131   v = n = (x-1)/(x+1)
  132   m = n*n
  133 
  134   /* Sum the series. */
  135   for (i=3; 1; i+=2) {
  136     e = (n *= m) / i
  137     if (e == 0) {
  138       v = f*v
  139       scale = z
  140       return (v/1)
  141     }
  142     v += e
  143   }
  144 }
  145 
  146 /* Sin(x)  uses the standard series:
  147    sin(x) = x - x^3/3! + x^5/5! - x^7/7! ... */
  148 
  149 define s(x) {
  150   auto  b, e, i, m, n, s, v, z
  151 
  152   /* Non base 10 ibase? */
  153   if (ibase != A) {
  154      b = ibase;
  155      ibase = A;
  156      v = s(x);
  157      ibase = b;
  158      return (v);
  159   }
  160 
  161   /* precondition x. */
  162   z = scale 
  163   scale = 1.1*z + 2;
  164   v = a(1)
  165   if (x < 0) {
  166     m = 1;
  167     x = -x;
  168   }
  169   scale = 0
  170   n = (x / v + 2 )/4
  171   x = x - 4*n*v
  172   if (n%2) x = -x
  173 
  174   /* Do the loop. */
  175   scale = z + 2;
  176   v = e = x
  177   s = -x*x
  178   for (i=3; 1; i+=2) {
  179     e *= s/(i*(i-1))
  180     if (e == 0) {
  181       scale = z
  182       if (m) return (-v/1);
  183       return (v/1);
  184     }
  185     v += e
  186   }
  187 }
  188 
  189 /* Cosine : cos(x) = sin(x+pi/2) */
  190 define c(x) {
  191   auto b, v, z;
  192 
  193   /* Non base 10 ibase? */
  194   if (ibase != A) {
  195      b = ibase;
  196      ibase = A;
  197      v = c(x);
  198      ibase = b;
  199      return (v);
  200   }
  201 
  202   z = scale;
  203   scale = scale*1.2;
  204   v = s(x+a(1)*2);
  205   scale = z;
  206   return (v/1);
  207 }
  208 
  209 /* Arctan: Using the formula:
  210      atan(x) = atan(c) + atan((x-c)/(1+xc)) for a small c (.2 here)
  211    For under .2, use the series:
  212      atan(x) = x - x^3/3 + x^5/5 - x^7/7 + ...   */
  213 
  214 define a(x) {
  215   auto a, b, e, f, i, m, n, s, v, z
  216 
  217   /* a is the value of a(.2) if it is needed. */
  218   /* f is the value to multiply by a in the return. */
  219   /* e is the value of the current term in the series. */
  220   /* v is the accumulated value of the series. */
  221   /* m is 1 or -1 depending on x (-x -> -1).  results are divided by m. */
  222   /* i is the denominator value for series element. */
  223   /* n is the numerator value for the series element. */
  224   /* s is -x*x. */
  225   /* z is the saved user's scale. */
  226 
  227   /* Non base 10 ibase? */
  228   if (ibase != A) {
  229      b = ibase;
  230      ibase = A;
  231      v = a(x);
  232      ibase = b;
  233      return (v);
  234   }
  235 
  236   /* Negative x? */
  237   m = 1;
  238   if (x<0) {
  239     m = -1;
  240     x = -x;
  241   }
  242 
  243   /* Special case and for fast answers */
  244   if (x==1) {
  245     if (scale <= 25) return (.7853981633974483096156608/m)
  246     if (scale <= 40) return (.7853981633974483096156608458198757210492/m)
  247     if (scale <= 60) \
  248       return (.785398163397448309615660845819875721049292349843776455243736/m)
  249   }
  250   if (x==.2) {
  251     if (scale <= 25) return (.1973955598498807583700497/m)
  252     if (scale <= 40) return (.1973955598498807583700497651947902934475/m)
  253     if (scale <= 60) \
  254       return (.197395559849880758370049765194790293447585103787852101517688/m)
  255   }
  256 
  257 
  258   /* Save the scale. */
  259   z = scale;
  260 
  261   /* Note: a and f are known to be zero due to being auto vars. */
  262   /* Calculate atan of a known number. */ 
  263   if (x > .2)  {
  264     scale = z+5;
  265     a = a(.2);
  266   }
  267    
  268   /* Precondition x. */
  269   scale = z+3;
  270   while (x > .2) {
  271     f += 1;
  272     x = (x-.2) / (1+x*.2);
  273   }
  274 
  275   /* Initialize the series. */
  276   v = n = x;
  277   s = -x*x;
  278 
  279   /* Calculate the series. */
  280   for (i=3; 1; i+=2) {
  281     e = (n *= s) / i;
  282     if (e == 0) {
  283       scale = z;
  284       return ((f*a+v)/m);
  285     }
  286     v += e
  287   }
  288 }
  289 
  290 
  291 /* Bessel function of integer order.  Uses the following:
  292    j(-n,x) = (-1)^n*j(n,x) 
  293    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))
  294             - x^6/(2^6*3!*(n+1)*(n+2)*(n+3)) .... )
  295 */
  296 define j(n,x) {
  297   auto a, b, d, e, f, i, m, s, v, z
  298 
  299   /* Non base 10 ibase? */
  300   if (ibase != A) {
  301      b = ibase;
  302      ibase = A;
  303      v = j(n,x);
  304      ibase = b;
  305      return (v);
  306   }
  307 
  308   /* Make n an integer and check for negative n. */
  309   z = scale;
  310   scale = 0;
  311   n = n/1;
  312   if (n<0) {
  313     n = -n;
  314     if (n%2 == 1) m = 1;
  315   }
  316 
  317   /* Compute the factor of x^n/(2^n*n!) */
  318   f = 1;
  319   for (i=2; i<=n; i++) f = f*i;
  320   scale = 1.5*z;
  321   f = x^n / 2^n / f;
  322 
  323   /* Initialize the loop .*/
  324   v = e = 1;
  325   s = -x*x/4
  326   scale = 1.5*z + length(f) - scale(f);
  327 
  328   /* The Loop.... */
  329   for (i=1; 1; i++) {
  330     e =  e * s / i / (n+i);
  331     if (e == 0) {
  332        scale = z
  333        if (m) return (-f*v/1);
  334        return (f*v/1);
  335     }
  336     v += e;
  337   }
  338 }