/* libmath.b for bc for minix. */
/* This file is part of bc written for MINIX.
Copyright (C) 1991 Free Software Foundation, Inc.
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License , or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
You may contact the author by:
us-mail: Philip A. Nelson
Computer Science Department, 9062
Western Washington University
Bellingham, WA 98226-9062
*************************************************************************/
/* Uses the fact that e^x = (e^(x/2))^2
When x is small enough, we use the series:
e^x = 1 + x + x^2/2! + x^3/3! + ...
auto a, d, e, f, i, m, v, z
/* Check the sign of x. */
/* Initialize the variables. */
if (f>0) while (f--) v = v*v;
/* Natural log. Uses the fact that ln(x^2) = 2*ln(x)
ln(x) = 2(a+a^3/3+a^5/5+...) where a=(x-1)/(x+1)
/* return something for the special case. */
if (x <= 0) return (1 - 10^scale)
/* Precondition x to make .5 < x < 2.0. */
while (x >= 2) { /* for large numbers */
while (x <= .5) { /* for small numbers */
/* Sin(x) uses the standard series:
sin(x) = x - x^3/3! + x^5/5! - x^7/7! ... */
/* Cosine : cos(x) = sin(x+pi/2) */
/* Arctan: Using the formula:
atan(x) = atan(c) + atan((x-c)/(1+xc)) for a small c (.2 here)
For under .2, use the series:
atan(x) = x - x^3/3 + x^5/5 - x^7/7 + ... */
auto a, e, f, i, m, n, s, v, z
/* Special case and for fast answers */
if (scale <= 25) return (.7853981633974483096156608/1)
if (scale <= 40) return (.7853981633974483096156608458198757210492/1)
return (.785398163397448309615660845819875721049292349843776455243736/1)
if (scale <= 25) return (.1973955598498807583700497/1)
if (scale <= 40) return (.1973955598498807583700497651947902934475/1)
return (.197395559849880758370049765194790293447585103787852101517688/1)
/* Note: a and f are known to be zero due to being auto vars. */
/* Calculate atan of a known number. */
/* Initialize the series. */
/* Calculate the series. */
if (m) return ((f*a+v)/-1);
/* Bessel function of integer order. Uses the following:
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))
- x^6/(2^6*3!*(n+1)*(n+2)*(n+3)) .... )
auto a, d, e, f, i, m, s, v, z
/* Make n an integer and check for negative n. */
/* Compute the factor of x^n/(2^n*n!) */
for (i=2; i<=n; i++) f = f*i;
/* Initialize the loop .*/