j0.c

newos is new operation system

/*-
 * Copyright (c) 1992, 1993
 *	The Regents of the University of California.  All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 * 3. All advertising materials mentioning features or use of this software
 *    must display the following acknowledgement:
 *	This product includes software developed by the University of
 *	California, Berkeley and its contributors.
 * 4. Neither the name of the University nor the names of its contributors
 *    may be used to endorse or promote products derived from this software
 *    without specific prior written permission.
 *
 * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
 * SUCH DAMAGE.
 */

#ifndef lint
static char sccsid[] = "@(#)j0.c	8.2 (Berkeley) 11/30/93";
#endif /* not lint */

/*
 * 16 December 1992
 * Minor modifications by Peter McIlroy to adapt non-IEEE architecture.
 */

/*
 * ====================================================
 * Copyright (C) 1992 by Sun Microsystems, Inc.
 *
 * Developed at SunPro, a Sun Microsystems, Inc. business.
 * Permission to use, copy, modify, and distribute this
 * software is freely granted, provided that this notice
 * is preserved.
 * ====================================================
 *
 * ******************* WARNING ********************
 * This is an alpha version of SunPro's FDLIBM (Freely
 * Distributable Math Library) for IEEE double precision
 * arithmetic. FDLIBM is a basic math library written
 * in C that runs on machines that conform to IEEE
 * Standard 754/854. This alpha version is distributed
 * for testing purpose. Those who use this software
 * should report any bugs to
 *
 *		[e-mail]
 *
 * -- K.C. Ng, Oct 12, 1992
 * ************************************************
 */

/* double j0(double x), y0(double x)
 * Bessel function of the first and second kinds of order zero.
 * Method -- j0(x):
 *	1. For tiny x, we use j0(x) = 1 - x^2/4 + x^4/64 - ...
 *	2. Reduce x to |x| since j0(x)=j0(-x),  and
 *	   for x in (0,2)
 *		j0(x) = 1-z/4+ z^2*R0/S0,  where z = x*x;
 *	   (precision:  |j0-1+z/4-z^2R0/S0 |<2**-63.67 )
 *	   for x in (2,inf)
 * 		j0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x0)-q0(x)*sin(x0))
 * 	   where x0 = x-pi/4. It is better to compute sin(x0),cos(x0)
 *	   as follow:
 *		cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4)
 *			= 1/sqrt(2) * (cos(x) + sin(x))
 *		sin(x0) = sin(x)cos(pi/4)-cos(x)sin(pi/4)
 *			= 1/sqrt(2) * (sin(x) - cos(x))
 * 	   (To avoid cancellation, use
 *		sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
 * 	    to compute the worse one.)
 *
 *	3 Special cases
 *		j0(nan)= nan
 *		j0(0) = 1
 *		j0(inf) = 0
 *
 * Method -- y0(x):
 *	1. For x<2.
 *	   Since
 *		y0(x) = 2/pi*(j0(x)*(ln(x/2)+Euler) + x^2/4 - ...)
 *	   therefore y0(x)-2/pi*j0(x)*ln(x) is an even function.
 *	   We use the following function to approximate y0,
 *		y0(x) = U(z)/V(z) + (2/pi)*(j0(x)*ln(x)), z= x^2
 *	   where
 *		U(z) = u0 + u1*z + ... + u6*z^6
 *		V(z) = 1  + v1*z + ... + v4*z^4
 *	   with absolute approximation error bounded by 2**-72.
 *	   Note: For tiny x, U/V = u0 and j0(x)~1, hence
 *		y0(tiny) = u0 + (2/pi)*ln(tiny), (choose tiny<2**-27)
 *	2. For x>=2.
 * 		y0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x0)+q0(x)*sin(x0))
 * 	   where x0 = x-pi/4. It is better to compute sin(x0),cos(x0)
 *	   by the method mentioned above.
 *	3. Special cases: y0(0)=-inf, y0(x<0)=NaN, y0(inf)=0.
 */

#include <math.h>
#include <float.h>
#if defined(vax) || defined(tahoe)
#define _IEEE	0
#else
#define _IEEE	1
#define infnan(x) (0.0)
#endif

static double pzero(double);
static double qzero(double);

static double const huge     = 1e300;
static double const zero     = 0.0;
static double const one      = 1.0;
static double const invsqrtpi= 5.641895835477562869480794515607725858441e-0001;
static double const tpi      = 0.636619772367581343075535053490057448;
 		/* R0/S0 on [0, 2.00] */
static double const r02 =   1.562499999999999408594634421055018003102e-0002;
static double const r03 =  -1.899792942388547334476601771991800712355e-0004;
static double const r04 =   1.829540495327006565964161150603950916854e-0006;
static double const r05 =  -4.618326885321032060803075217804816988758e-0009;
static double const s01 =   1.561910294648900170180789369288114642057e-0002;
static double const s02 =   1.169267846633374484918570613449245536323e-0004;
static double const s03 =   5.135465502073181376284426245689510134134e-0007;
static double const s04 =   1.166140033337900097836930825478674320464e-0009;

double
j0(double x)
{
	double z;
	double s;
	double c;
	double ss;
	double cc;
	double r;
	double u;
	double v;

	if (!finite(x)) {
		if (_IEEE) return one/(x*x);
		else return (0);
	}
	x = fabs(x);
	if (x >= 2.0) {	/* |x| >= 2.0 */
		s = sin(x);
		c = cos(x);
		ss = s-c;
		cc = s+c;
		if (x < .5 * DBL_MAX) {  /* make sure x+x not overflow */
		    z = -cos(x+x);
		    if ((s*c)<zero) cc = z/ss;
		    else 	    ss = z/cc;
		}
	/*
	 * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x)
	 * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x)
	 */
		if (_IEEE && x> 6.80564733841876927e+38) /* 2^129 */
			z = (invsqrtpi*cc)/sqrt(x);
		else {
		    u = pzero(x); v = qzero(x);
		    z = invsqrtpi*(u*cc-v*ss)/sqrt(x);
		}
		return z;
	}
	if (x < 1.220703125e-004) {		   /* |x| < 2**-13 */
	    if (huge+x > one) {			   /* raise inexact if x != 0 */
	        if (x < 7.450580596923828125e-009) /* |x|<2**-27 */
			return one;
	        else return (one - 0.25*x*x);
	    }
	}
	z = x*x;
	r =  z*(r02+z*(r03+z*(r04+z*r05)));
	s =  one+z*(s01+z*(s02+z*(s03+z*s04)));
	if (x < one) {			/* |x| < 1.00 */
	    return (one + z*(-0.25+(r/s)));
	} else {
	    u = 0.5*x;
	    return ((one+u)*(one-u)+z*(r/s));
	}
}

static double const u00 =  -7.380429510868722527422411862872999615628e-0002;
static double const u01 =   1.766664525091811069896442906220827182707e-0001;
static double const u02 =  -1.381856719455968955440002438182885835344e-0002;
static double const u03 =   3.474534320936836562092566861515617053954e-0004;
static double const u04 =  -3.814070537243641752631729276103284491172e-0006;
static double const u05 =   1.955901370350229170025509706510038090009e-0008;
static double const u06 =  -3.982051941321034108350630097330144576337e-0011;
static double const v01 =   1.273048348341237002944554656529224780561e-0002;
static double const v02 =   7.600686273503532807462101309675806839635e-0005;
static double const v03 =   2.591508518404578033173189144579208685163e-0007;
static double const v04 =   4.411103113326754838596529339004302243157e-0010;

double
y0(double x)
{
	double z;
	double s;
	double c;
	double ss;
	double cc;
	double u;
	double v;
    /* Y0(NaN) is NaN, y0(-inf) is Nan, y0(inf) is 0  */
	if (!finite(x)) {
		if (_IEEE)
			return (one/(x+x*x));
		else
			return (0);
	}
        if (x == 0) {
		if (_IEEE)	return (-one/zero);
		else		return(infnan(-ERANGE));
	}
        if (x<0) {
		if (_IEEE)	return (zero/zero);
		else		return (infnan(EDOM));
	}