[unix-history] / usr / src / usr.lib / libom / j1.c

/*	@(#)j1.c	4.1	12/25/82	*/

/*
	floating point Bessel's function
	of the first and second kinds
	of order one

	j1(x) returns the value of J1(x)
	for all real values of x.

	There are no error returns.
	Calls sin, cos, sqrt.

	There is a niggling bug in J1 which
	causes errors up to 2e-16 for x in the
	interval [-8,8].
	The bug is caused by an inappropriate order
	of summation of the series.  rhm will fix it
	someday.

	Coefficients are from Hart & Cheney.
	#6050 (20.98D)
	#6750 (19.19D)
	#7150 (19.35D)

	y1(x) returns the value of Y1(x)
	for positive real values of x.
	For x<=0, error number EDOM is set and a
	large negative value is returned.

	Calls sin, cos, sqrt, log, j1.

	The values of Y1 have not been checked
	to more than ten places.

	Coefficients are from Hart & Cheney.
	#6447 (22.18D)
	#6750 (19.19D)
	#7150 (19.35D)
*/

#include <math.h>
#include <errno.h>

int	errno;
static double pzero, qzero;
static double tpi	= .6366197723675813430755350535e0;
static double pio4	= .7853981633974483096156608458e0;
static double p1[] = {
	0.581199354001606143928050809e21,
	-.6672106568924916298020941484e20,
	0.2316433580634002297931815435e19,
	-.3588817569910106050743641413e17,
	0.2908795263834775409737601689e15,
	-.1322983480332126453125473247e13,
	0.3413234182301700539091292655e10,
	-.4695753530642995859767162166e7,
	0.2701122710892323414856790990e4,
};
static double q1[] = {
	0.1162398708003212287858529400e22,
	0.1185770712190320999837113348e20,
	0.6092061398917521746105196863e17,
	0.2081661221307607351240184229e15,
	0.5243710262167649715406728642e12,
	0.1013863514358673989967045588e10,
	0.1501793594998585505921097578e7,
	0.1606931573481487801970916749e4,
	1.0,
};
static double p2[] = {
	-.4435757816794127857114720794e7,
	-.9942246505077641195658377899e7,
	-.6603373248364939109255245434e7,
	-.1523529351181137383255105722e7,
	-.1098240554345934672737413139e6,
	-.1611616644324610116477412898e4,
	0.0,
};
static double q2[] = {
	-.4435757816794127856828016962e7,
	-.9934124389934585658967556309e7,
	-.6585339479723087072826915069e7,
	-.1511809506634160881644546358e7,
	-.1072638599110382011903063867e6,
	-.1455009440190496182453565068e4,
	1.0,
};
static double p3[] = {
	0.3322091340985722351859704442e5,
	0.8514516067533570196555001171e5,
	0.6617883658127083517939992166e5,
	0.1849426287322386679652009819e5,
	0.1706375429020768002061283546e4,
	0.3526513384663603218592175580e2,
	0.0,
};
static double q3[] = {
	0.7087128194102874357377502472e6,
	0.1819458042243997298924553839e7,
	0.1419460669603720892855755253e7,
	0.4002944358226697511708610813e6,
	0.3789022974577220264142952256e5,
	0.8638367769604990967475517183e3,
	1.0,
};
static double p4[] = {
	-.9963753424306922225996744354e23,
	0.2655473831434854326894248968e23,
	-.1212297555414509577913561535e22,
	0.2193107339917797592111427556e20,
	-.1965887462722140658820322248e18,
	0.9569930239921683481121552788e15,
	-.2580681702194450950541426399e13,
	0.3639488548124002058278999428e10,
	-.2108847540133123652824139923e7,
	0.0,
};
static double q4[] = {
	0.5082067366941243245314424152e24,
	0.5435310377188854170800653097e22,
	0.2954987935897148674290758119e20,
	0.1082258259408819552553850180e18,
	0.2976632125647276729292742282e15,
	0.6465340881265275571961681500e12,
	0.1128686837169442121732366891e10,
	0.1563282754899580604737366452e7,
	0.1612361029677000859332072312e4,
	1.0,
};

double
j1(arg) double arg;{
	double xsq, n, d, x;
	double sin(), cos(), sqrt();
	int i;

	x = arg;
	if(x < 0.) x = -x;
	if(x > 8.){
		asympt(x);
		n = x - 3.*pio4;
		n = sqrt(tpi/x)*(pzero*cos(n) - qzero*sin(n));
		if(arg <0.) n = -n;
		return(n);
	}
	xsq = x*x;
	for(n=0,d=0,i=8;i>=0;i--){
		n = n*xsq + p1[i];
		d = d*xsq + q1[i];
	}
	return(arg*n/d);
}

double
y1(arg) double arg;{
	double xsq, n, d, x;
	double sin(), cos(), sqrt(), log(), j1();
	int i;

	errno = 0;
	x = arg;
	if(x <= 0.){
		errno = EDOM;
		return(-HUGE);
	}
	if(x > 8.){
		asympt(x);
		n = x - 3*pio4;
		return(sqrt(tpi/x)*(pzero*sin(n) + qzero*cos(n)));
	}
	xsq = x*x;
	for(n=0,d=0,i=9;i>=0;i--){
		n = n*xsq + p4[i];
		d = d*xsq + q4[i];
	}
	return(x*n/d + tpi*(j1(x)*log(x)-1./x));
}

static
asympt(arg) double arg;{
	double zsq, n, d;
	int i;
	zsq = 64./(arg*arg);
	for(n=0,d=0,i=6;i>=0;i--){
		n = n*zsq + p2[i];
		d = d*zsq + q2[i];
	}
	pzero = n/d;
	for(n=0,d=0,i=6;i>=0;i--){
		n = n*zsq + p3[i];
		d = d*zsq + q3[i];
	}
	qzero = (8./arg)*(n/d);
}
Commit	Line	Data
e804469b	1	/* @(#)j1.c 4.1 12/25/82 */
3c2bcac9 SL	2
	3	/*
	4	floating point Bessel's function
	5	of the first and second kinds
	6	of order one
	7
	8	j1(x) returns the value of J1(x)
	9	for all real values of x.
	10
	11	There are no error returns.
	12	Calls sin, cos, sqrt.
	13
	14	There is a niggling bug in J1 which
	15	causes errors up to 2e-16 for x in the
	16	interval [-8,8].
	17	The bug is caused by an inappropriate order
	18	of summation of the series. rhm will fix it
	19	someday.
	20
	21	Coefficients are from Hart & Cheney.
	22	#6050 (20.98D)
	23	#6750 (19.19D)
	24	#7150 (19.35D)
	25
	26	y1(x) returns the value of Y1(x)
	27	for positive real values of x.
	28	For x<=0, error number EDOM is set and a
	29	large negative value is returned.
	30
	31	Calls sin, cos, sqrt, log, j1.
	32
	33	The values of Y1 have not been checked
	34	to more than ten places.
	35
	36	Coefficients are from Hart & Cheney.
	37	#6447 (22.18D)
	38	#6750 (19.19D)
	39	#7150 (19.35D)
	40	*/
	41
	42	#include <math.h>
	43	#include <errno.h>
	44
	45	int errno;
	46	static double pzero, qzero;
	47	static double tpi = .6366197723675813430755350535e0;
	48	static double pio4 = .7853981633974483096156608458e0;
	49	static double p1[] = {
	50	0.581199354001606143928050809e21,
	51	-.6672106568924916298020941484e20,
	52	0.2316433580634002297931815435e19,
	53	-.3588817569910106050743641413e17,
	54	0.2908795263834775409737601689e15,
	55	-.1322983480332126453125473247e13,
	56	0.3413234182301700539091292655e10,
	57	-.4695753530642995859767162166e7,
	58	0.2701122710892323414856790990e4,
	59	};
	60	static double q1[] = {
	61	0.1162398708003212287858529400e22,
	62	0.1185770712190320999837113348e20,
	63	0.6092061398917521746105196863e17,
	64	0.2081661221307607351240184229e15,
	65	0.5243710262167649715406728642e12,
66	0.1013863514358673989967045588e10,
67	0.1501793594998585505921097578e7,
68	0.1606931573481487801970916749e4,
69	1.0,
70	};
71	static double p2[] = {
72	-.4435757816794127857114720794e7,
73	-.9942246505077641195658377899e7,
74	-.6603373248364939109255245434e7,
75	-.1523529351181137383255105722e7,
76	-.1098240554345934672737413139e6,
77	-.1611616644324610116477412898e4,
78	0.0,
79	};
80	static double q2[] = {
81	-.4435757816794127856828016962e7,
82	-.9934124389934585658967556309e7,
83	-.6585339479723087072826915069e7,
84	-.1511809506634160881644546358e7,
85	-.1072638599110382011903063867e6,
86	-.1455009440190496182453565068e4,
87	1.0,
88	};
89	static double p3[] = {
90	0.3322091340985722351859704442e5,
91	0.8514516067533570196555001171e5,
92	0.6617883658127083517939992166e5,
93	0.1849426287322386679652009819e5,
94	0.1706375429020768002061283546e4,
95	0.3526513384663603218592175580e2,
96	0.0,
97	};
98	static double q3[] = {
99	0.7087128194102874357377502472e6,
100	0.1819458042243997298924553839e7,
101	0.1419460669603720892855755253e7,
102	0.4002944358226697511708610813e6,
103	0.3789022974577220264142952256e5,
104	0.8638367769604990967475517183e3,
105	1.0,
106	};
107	static double p4[] = {
108	-.9963753424306922225996744354e23,
109	0.2655473831434854326894248968e23,
110	-.1212297555414509577913561535e22,
111	0.2193107339917797592111427556e20,
112	-.1965887462722140658820322248e18,
113	0.9569930239921683481121552788e15,
114	-.2580681702194450950541426399e13,
115	0.3639488548124002058278999428e10,
116	-.2108847540133123652824139923e7,
117	0.0,
118	};
119	static double q4[] = {
120	0.5082067366941243245314424152e24,
121	0.5435310377188854170800653097e22,
122	0.2954987935897148674290758119e20,
123	0.1082258259408819552553850180e18,
124	0.2976632125647276729292742282e15,
125	0.6465340881265275571961681500e12,
126	0.1128686837169442121732366891e10,
127	0.1563282754899580604737366452e7,
128	0.1612361029677000859332072312e4,
129	1.0,
130	};
131
132	double
133	j1(arg) double arg;{
134	double xsq, n, d, x;
135	double sin(), cos(), sqrt();
136	int i;
137
138	x = arg;
139	if(x < 0.) x = -x;
140	if(x > 8.){
141	asympt(x);
142	n = x - 3.*pio4;
143	n = sqrt(tpi/x)(pzerocos(n) - qzero*sin(n));
144	if(arg <0.) n = -n;
145	return(n);
146	}
147	xsq = x*x;
148	for(n=0,d=0,i=8;i>=0;i--){
149	n = n*xsq + p1[i];
150	d = d*xsq + q1[i];
151	}
152	return(arg*n/d);
153	}
154
155	double
156	y1(arg) double arg;{
157	double xsq, n, d, x;
158	double sin(), cos(), sqrt(), log(), j1();
159	int i;
160
161	errno = 0;
162	x = arg;
163	if(x <= 0.){
164	errno = EDOM;
165	return(-HUGE);
166	}
167	if(x > 8.){
168	asympt(x);
169	n = x - 3*pio4;
170	return(sqrt(tpi/x)(pzerosin(n) + qzero*cos(n)));
171	}
172	xsq = x*x;
173	for(n=0,d=0,i=9;i>=0;i--){
174	n = n*xsq + p4[i];
175	d = d*xsq + q4[i];
176	}
177	return(xn/d + tpi(j1(x)*log(x)-1./x));
178	}
179
180	static
181	asympt(arg) double arg;{
182	double zsq, n, d;
183	int i;
184	zsq = 64./(arg*arg);
185	for(n=0,d=0,i=6;i>=0;i--){
186	n = n*zsq + p2[i];
187	d = d*zsq + q2[i];
188	}
189	pzero = n/d;
190	for(n=0,d=0,i=6;i>=0;i--){
191	n = n*zsq + p3[i];
192	d = d*zsq + q3[i];
193	}
194	qzero = (8./arg)*(n/d);
195	}