j1.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[] = "@(#)j1.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 j1(double x), y1(double x)
 * Bessel function of the first and second kinds of order zero.
 * Method -- j1(x):
 *	1. For tiny x, we use j1(x) = x/2 - x^3/16 + x^5/384 - ...
 *	2. Reduce x to |x| since j1(x)=-j1(-x),  and
 *	   for x in (0,2)
 *		j1(x) = x/2 + x*z*R0/S0,  where z = x*x;
 *	   (precision:  |j1/x - 1/2 - R0/S0 |<2**-61.51 )
 *	   for x in (2,inf)
 * 		j1(x) = sqrt(2/(pi*x))*(p1(x)*cos(x1)-q1(x)*sin(x1))
 * 		y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x1)+q1(x)*cos(x1))
 * 	   where x1 = x-3*pi/4. It is better to compute sin(x1),cos(x1)
 *	   as follows:
 *		cos(x1) =  cos(x)cos(3pi/4)+sin(x)sin(3pi/4)
 *			=  1/sqrt(2) * (sin(x) - cos(x))
 *		sin(x1) =  sin(x)cos(3pi/4)-cos(x)sin(3pi/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
 *		j1(nan)= nan
 *		j1(0) = 0
 *		j1(inf) = 0
 *
 * Method -- y1(x):
 *	1. screen out x<=0 cases: y1(0)=-inf, y1(x<0)=NaN
 *	2. For x<2.
 *	   Since
 *		y1(x) = 2/pi*(j1(x)*(ln(x/2)+Euler)-1/x-x/2+5/64*x^3-...)
 *	   therefore y1(x)-2/pi*j1(x)*ln(x)-1/x is an odd function.
 *	   We use the following function to approximate y1,
 *		y1(x) = x*U(z)/V(z) + (2/pi)*(j1(x)*ln(x)-1/x), z= x^2
 *	   where for x in [0,2] (abs err less than 2**-65.89)
 *		U(z) = u0 + u1*z + ... + u4*z^4
 *		V(z) = 1  + v1*z + ... + v5*z^5
 *	   Note: For tiny x, 1/x dominate y1 and hence
 *		y1(tiny) = -2/pi/tiny, (choose tiny<2**-54)
 *	3. For x>=2.
 * 		y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x1)+q1(x)*cos(x1))
 * 	   where x1 = x-3*pi/4. It is better to compute sin(x1),cos(x1)
 *	   by method mentioned above.
 */

#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 pone(double);
static double qone(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] */
static double const r00 =  -6.250000000000000020842322918309200910191e-0002;
static double const r01 =   1.407056669551897148204830386691427791200e-0003;
static double const r02 =  -1.599556310840356073980727783817809847071e-0005;
static double const r03 =   4.967279996095844750387702652791615403527e-0008;
static double const s01 =   1.915375995383634614394860200531091839635e-0002;
static double const s02 =   1.859467855886309024045655476348872850396e-0004;
static double const s03 =   1.177184640426236767593432585906758230822e-0006;
static double const s04 =   5.046362570762170559046714468225101016915e-0009;
static double const s05 =   1.235422744261379203512624973117299248281e-0011;

#define two_129	6.80564733841876926e+038	/* 2^129 */
#define two_m54	5.55111512312578270e-017	/* 2^-54 */
double
j1(double x)
{
	double z;
	double s;
	double c;
	double ss;
	double cc;
	double r;
	double u;
	double v;
	double y;

	y = fabs(x);
	if (!finite(x)) {		/* Inf or NaN */
		if (_IEEE && x != x)
			return(x);
		else
			return (copysign(x, zero));
	}
	y = fabs(x);
	if (y >= 2) {			/* |x| >= 2.0 */
		s = sin(y);
		c = cos(y);
		ss = -s-c;
		cc = s-c;
		if (y < .5*DBL_MAX) {  	/* make sure y+y not overflow */
		    z = cos(y+y);
		    if ((s*c)<zero) cc = z/ss;
		    else 	    ss = z/cc;
		}
	/*
	 * j1(x) = 1/sqrt(pi) * (P(1,x)*cc - Q(1,x)*ss) / sqrt(x)
	 * y1(x) = 1/sqrt(pi) * (P(1,x)*ss + Q(1,x)*cc) / sqrt(x)
	 */
#if !defined(vax) && !defined(tahoe)
		if (y > two_129)	 /* x > 2^129 */
			z = (invsqrtpi*cc)/sqrt(y);
		else
#endif /* defined(vax) || defined(tahoe) */
		{
		    u = pone(y); v = qone(y);
		    z = invsqrtpi*(u*cc-v*ss)/sqrt(y);
		}
		if (x < 0) return -z;
		else  	 return  z;
	}
	if (y < 7.450580596923828125e-009) {	/* |x|<2**-27 */
	    if(huge+x>one) return 0.5*x;/* inexact if x!=0 necessary */
	}
	z = x*x;
	r =  z*(r00+z*(r01+z*(r02+z*r03)));
	s =  one+z*(s01+z*(s02+z*(s03+z*(s04+z*s05))));
	r *= x;
	return (x*0.5+r/s);
}

static double const u0[5] = {
  -1.960570906462389484206891092512047539632e-0001,
   5.044387166398112572026169863174882070274e-0002,
  -1.912568958757635383926261729464141209569e-0003,
   2.352526005616105109577368905595045204577e-0005,
   -9.190991580398788465315411784276789663849e-0008,
};
static double const v0[5] = {
   1.991673182366499064031901734535479833387e-0002,
   2.025525810251351806268483867032781294682e-0004,
   1.356088010975162198085369545564475416398e-0006,
   6.227414523646214811803898435084697863445e-0009,
   1.665592462079920695971450872592458916421e-0011,
};

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

    /* if Y1(NaN) is NaN, Y1(-inf) is NaN, Y1(inf) is 0 */
	if (!finite(x)) {
		if (!_IEEE) return (infnan(EDOM));
		else if (x < 0)
			return(zero/zero);
		else if (x > 0)
			return (0);
		else
			return(x);
	}