![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to j0.c CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [CSRG BSD Unix] / 43BSDReno / lib / libm / common_source|
Return to j0.c CVS log | Up to [CSRG BSD Unix] / 43BSDReno / lib / libm / common_source |
1.1 root 1: /*
2: * Copyright (c) 1985 Regents of the University of California.
3: * All rights reserved. The Berkeley software License Agreement
4: * specifies the terms and conditions for redistribution.
5: */
6:
7: #ifndef lint
8: static char sccsid[] = "@(#)j0.c 5.3 (Berkeley) 9/22/88";
9: #endif /* not lint */
10:
11: /*
12: floating point Bessel's function
13: of the first and second kinds
14: of order zero
15:
16: j0(x) returns the value of J0(x)
17: for all real values of x.
18:
19: There are no error returns.
20: Calls sin, cos, sqrt.
21:
22: There is a niggling bug in J0 which
23: causes errors up to 2e-16 for x in the
24: interval [-8,8].
25: The bug is caused by an inappropriate order
26: of summation of the series. rhm will fix it
27: someday.
28:
29: Coefficients are from Hart & Cheney.
30: #5849 (19.22D)
31: #6549 (19.25D)
32: #6949 (19.41D)
33:
34: y0(x) returns the value of Y0(x)
35: for positive real values of x.
36: For x<=0, if on the VAX, error number EDOM is set and
37: the reserved operand fault is generated;
38: otherwise (an IEEE machine) an invalid operation is performed.
39:
40: Calls sin, cos, sqrt, log, j0.
41:
42: The values of Y0 have not been checked
43: to more than ten places.
44:
45: Coefficients are from Hart & Cheney.
46: #6245 (18.78D)
47: #6549 (19.25D)
48: #6949 (19.41D)
49: */
50:
51: #include "mathimpl.h"
52:
53: #if defined(vax)||defined(tahoe)
54: #include <errno.h>
55: #else /* defined(vax)||defined(tahoe) */
56: static const double zero = 0.e0;
57: #endif /* defined(vax)||defined(tahoe) */
58:
59: static double pzero, qzero;
60:
61: static const double tpi = .6366197723675813430755350535e0;
62: static const double pio4 = .7853981633974483096156608458e0;
63: static const double p1[] = {
64: 0.4933787251794133561816813446e21,
65: -.1179157629107610536038440800e21,
66: 0.6382059341072356562289432465e19,
67: -.1367620353088171386865416609e18,
68: 0.1434354939140344111664316553e16,
69: -.8085222034853793871199468171e13,
70: 0.2507158285536881945555156435e11,
71: -.4050412371833132706360663322e8,
72: 0.2685786856980014981415848441e5,
73: };
74: static const double q1[] = {
75: 0.4933787251794133562113278438e21,
76: 0.5428918384092285160200195092e19,
77: 0.3024635616709462698627330784e17,
78: 0.1127756739679798507056031594e15,
79: 0.3123043114941213172572469442e12,
80: 0.6699987672982239671814028660e9,
81: 0.1114636098462985378182402543e7,
82: 0.1363063652328970604442810507e4,
83: 1.0
84: };
85: static const double p2[] = {
86: 0.5393485083869438325262122897e7,
87: 0.1233238476817638145232406055e8,
88: 0.8413041456550439208464315611e7,
89: 0.2016135283049983642487182349e7,
90: 0.1539826532623911470917825993e6,
91: 0.2485271928957404011288128951e4,
92: 0.0,
93: };
94: static const double q2[] = {
95: 0.5393485083869438325560444960e7,
96: 0.1233831022786324960844856182e8,
97: 0.8426449050629797331554404810e7,
98: 0.2025066801570134013891035236e7,
99: 0.1560017276940030940592769933e6,
100: 0.2615700736920839685159081813e4,
101: 1.0,
102: };
103: static const double p3[] = {
104: -.3984617357595222463506790588e4,
105: -.1038141698748464093880530341e5,
106: -.8239066313485606568803548860e4,
107: -.2365956170779108192723612816e4,
108: -.2262630641933704113967255053e3,
109: -.4887199395841261531199129300e1,
110: 0.0,
111: };
112: static const double q3[] = {
113: 0.2550155108860942382983170882e6,
114: 0.6667454239319826986004038103e6,
115: 0.5332913634216897168722255057e6,
116: 0.1560213206679291652539287109e6,
117: 0.1570489191515395519392882766e5,
118: 0.4087714673983499223402830260e3,
119: 1.0,
120: };
121: static const double p4[] = {
122: -.2750286678629109583701933175e20,
123: 0.6587473275719554925999402049e20,
124: -.5247065581112764941297350814e19,
125: 0.1375624316399344078571335453e18,
126: -.1648605817185729473122082537e16,
127: 0.1025520859686394284509167421e14,
128: -.3436371222979040378171030138e11,
129: 0.5915213465686889654273830069e8,
130: -.4137035497933148554125235152e5,
131: };
132: static const double q4[] = {
133: 0.3726458838986165881989980e21,
134: 0.4192417043410839973904769661e19,
135: 0.2392883043499781857439356652e17,
136: 0.9162038034075185262489147968e14,
137: 0.2613065755041081249568482092e12,
138: 0.5795122640700729537480087915e9,
139: 0.1001702641288906265666651753e7,
140: 0.1282452772478993804176329391e4,
141: 1.0,
142: };
143:
144: static void asympt();
145:
146: double
147: j0(arg) double arg;{
148: double argsq, n, d;
149: int i;
150:
151: if(arg < 0.) arg = -arg;
152: if(arg > 8.){
153: asympt(arg);
154: n = arg - pio4;
155: return(sqrt(tpi/arg)*(pzero*cos(n) - qzero*sin(n)));
156: }
157: argsq = arg*arg;
158: for(n=0,d=0,i=8;i>=0;i--){
159: n = n*argsq + p1[i];
160: d = d*argsq + q1[i];
161: }
162: return(n/d);
163: }
164:
165: double
166: y0(arg) double arg;{
167: double argsq, n, d;
168: int i;
169:
170: if(arg <= 0.){
171: #if defined(vax)||defined(tahoe)
172: return(infnan(EDOM)); /* NaN */
173: #else /* defined(vax)||defined(tahoe) */
174: return(zero/zero); /* IEEE machines: invalid operation */
175: #endif /* defined(vax)||defined(tahoe) */
176: }
177: if(arg > 8.){
178: asympt(arg);
179: n = arg - pio4;
180: return(sqrt(tpi/arg)*(pzero*sin(n) + qzero*cos(n)));
181: }
182: argsq = arg*arg;
183: for(n=0,d=0,i=8;i>=0;i--){
184: n = n*argsq + p4[i];
185: d = d*argsq + q4[i];
186: }
187: return(n/d + tpi*j0(arg)*log(arg));
188: }
189:
190: static void
191: asympt(arg) double arg;{
192: double zsq, n, d;
193: int i;
194: zsq = 64./(arg*arg);
195: for(n=0,d=0,i=6;i>=0;i--){
196: n = n*zsq + p2[i];
197: d = d*zsq + q2[i];
198: }
199: pzero = n/d;
200: for(n=0,d=0,i=6;i>=0;i--){
201: n = n*zsq + p3[i];
202: d = d*zsq + q3[i];
203: }
204: qzero = (8./arg)*(n/d);
205: }
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.