![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to j1.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 j1.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[] = "@(#)j1.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 one
15:
16: j1(x) returns the value of J1(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 J1 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: #6050 (20.98D)
31: #6750 (19.19D)
32: #7150 (19.35D)
33:
34: y1(x) returns the value of Y1(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, j1.
41:
42: The values of Y1 have not been checked
43: to more than ten places.
44:
45: Coefficients are from Hart & Cheney.
46: #6447 (22.18D)
47: #6750 (19.19D)
48: #7150 (19.35D)
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.581199354001606143928050809e21,
65: -.6672106568924916298020941484e20,
66: 0.2316433580634002297931815435e19,
67: -.3588817569910106050743641413e17,
68: 0.2908795263834775409737601689e15,
69: -.1322983480332126453125473247e13,
70: 0.3413234182301700539091292655e10,
71: -.4695753530642995859767162166e7,
72: 0.2701122710892323414856790990e4,
73: };
74: static const double q1[] = {
75: 0.1162398708003212287858529400e22,
76: 0.1185770712190320999837113348e20,
77: 0.6092061398917521746105196863e17,
78: 0.2081661221307607351240184229e15,
79: 0.5243710262167649715406728642e12,
80: 0.1013863514358673989967045588e10,
81: 0.1501793594998585505921097578e7,
82: 0.1606931573481487801970916749e4,
83: 1.0,
84: };
85: static const double p2[] = {
86: -.4435757816794127857114720794e7,
87: -.9942246505077641195658377899e7,
88: -.6603373248364939109255245434e7,
89: -.1523529351181137383255105722e7,
90: -.1098240554345934672737413139e6,
91: -.1611616644324610116477412898e4,
92: 0.0,
93: };
94: static const double q2[] = {
95: -.4435757816794127856828016962e7,
96: -.9934124389934585658967556309e7,
97: -.6585339479723087072826915069e7,
98: -.1511809506634160881644546358e7,
99: -.1072638599110382011903063867e6,
100: -.1455009440190496182453565068e4,
101: 1.0,
102: };
103: static const double p3[] = {
104: 0.3322091340985722351859704442e5,
105: 0.8514516067533570196555001171e5,
106: 0.6617883658127083517939992166e5,
107: 0.1849426287322386679652009819e5,
108: 0.1706375429020768002061283546e4,
109: 0.3526513384663603218592175580e2,
110: 0.0,
111: };
112: static const double q3[] = {
113: 0.7087128194102874357377502472e6,
114: 0.1819458042243997298924553839e7,
115: 0.1419460669603720892855755253e7,
116: 0.4002944358226697511708610813e6,
117: 0.3789022974577220264142952256e5,
118: 0.8638367769604990967475517183e3,
119: 1.0,
120: };
121: static const double p4[] = {
122: -.9963753424306922225996744354e23,
123: 0.2655473831434854326894248968e23,
124: -.1212297555414509577913561535e22,
125: 0.2193107339917797592111427556e20,
126: -.1965887462722140658820322248e18,
127: 0.9569930239921683481121552788e15,
128: -.2580681702194450950541426399e13,
129: 0.3639488548124002058278999428e10,
130: -.2108847540133123652824139923e7,
131: 0.0,
132: };
133: static const double q4[] = {
134: 0.5082067366941243245314424152e24,
135: 0.5435310377188854170800653097e22,
136: 0.2954987935897148674290758119e20,
137: 0.1082258259408819552553850180e18,
138: 0.2976632125647276729292742282e15,
139: 0.6465340881265275571961681500e12,
140: 0.1128686837169442121732366891e10,
141: 0.1563282754899580604737366452e7,
142: 0.1612361029677000859332072312e4,
143: 1.0,
144: };
145:
146: static void asympt();
147:
148: double
149: j1(arg) double arg;{
150: double xsq, n, d, x;
151: int i;
152:
153: x = arg;
154: if(x < 0.) x = -x;
155: if(x > 8.){
156: asympt(x);
157: n = x - 3.*pio4;
158: n = sqrt(tpi/x)*(pzero*cos(n) - qzero*sin(n));
159: if(arg <0.) n = -n;
160: return(n);
161: }
162: xsq = x*x;
163: for(n=0,d=0,i=8;i>=0;i--){
164: n = n*xsq + p1[i];
165: d = d*xsq + q1[i];
166: }
167: return(arg*n/d);
168: }
169:
170: double
171: y1(arg) double arg;{
172: double xsq, n, d, x;
173: int i;
174:
175: x = arg;
176: if(x <= 0.){
177: #if defined(vax)||defined(tahoe)
178: return(infnan(EDOM)); /* NaN */
179: #else /* defined(vax)||defined(tahoe) */
180: return(zero/zero); /* IEEE machines: invalid operation */
181: #endif /* defined(vax)||defined(tahoe) */
182: }
183: if(x > 8.){
184: asympt(x);
185: n = x - 3*pio4;
186: return(sqrt(tpi/x)*(pzero*sin(n) + qzero*cos(n)));
187: }
188: xsq = x*x;
189: for(n=0,d=0,i=9;i>=0;i--){
190: n = n*xsq + p4[i];
191: d = d*xsq + q4[i];
192: }
193: return(x*n/d + tpi*(j1(x)*log(x)-1./x));
194: }
195:
196: static void
197: asympt(arg) double arg;{
198: double zsq, n, d;
199: int i;
200: zsq = 64./(arg*arg);
201: for(n=0,d=0,i=6;i>=0;i--){
202: n = n*zsq + p2[i];
203: d = d*zsq + q2[i];
204: }
205: pzero = n/d;
206: for(n=0,d=0,i=6;i>=0;i--){
207: n = n*zsq + p3[i];
208: d = d*zsq + q3[i];
209: }
210: qzero = (8./arg)*(n/d);
211: }
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.