Print this page
5261 libm should stop using synonyms.h
5298 fabs is 0-sized, confuses dis(1) and others
Reviewed by: Josef 'Jeff' Sipek <jeffpc@josefsipek.net>
Approved by: Gordon Ross <gwr@nexenta.com>
| Split |
Close |
| Expand all |
| Collapse all |
--- old/usr/src/lib/libm/common/C/j0.c
+++ new/usr/src/lib/libm/common/C/j0.c
1 1 /*
2 2 * CDDL HEADER START
3 3 *
4 4 * The contents of this file are subject to the terms of the
5 5 * Common Development and Distribution License (the "License").
6 6 * You may not use this file except in compliance with the License.
7 7 *
8 8 * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
9 9 * or http://www.opensolaris.org/os/licensing.
10 10 * See the License for the specific language governing permissions
11 11 * and limitations under the License.
12 12 *
13 13 * When distributing Covered Code, include this CDDL HEADER in each
14 14 * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
15 15 * If applicable, add the following below this CDDL HEADER, with the
16 16 * fields enclosed by brackets "[]" replaced with your own identifying
17 17 * information: Portions Copyright [yyyy] [name of copyright owner]
18 18 *
19 19 * CDDL HEADER END
20 20 */
21 21
22 22 /*
23 23 * Copyright 2011 Nexenta Systems, Inc. All rights reserved.
24 24 */
25 25 /*
26 26 * Copyright 2006 Sun Microsystems, Inc. All rights reserved.
27 27 * Use is subject to license terms.
28 28 */
|
↓ open down ↓ |
28 lines elided |
↑ open up ↑ |
29 29
30 30 /*
31 31 * Floating point Bessel's function of the first and second kinds
32 32 * of order zero: j0(x),y0(x);
33 33 *
34 34 * Special cases:
35 35 * y0(0)=y1(0)=yn(n,0) = -inf with division by zero signal;
36 36 * y0(-ve)=y1(-ve)=yn(n,-ve) are NaN with invalid signal.
37 37 */
38 38
39 -#pragma weak j0 = __j0
40 -#pragma weak y0 = __y0
39 +#pragma weak __j0 = j0
40 +#pragma weak __y0 = y0
41 41
42 42 #include "libm.h"
43 -#include "libm_synonyms.h"
44 43 #include "libm_protos.h"
45 44 #include <math.h>
46 45 #include <values.h>
47 46
48 47 #define GENERIC double
49 48 static const GENERIC
50 49 zero = 0.0,
51 50 small = 1.0e-5,
52 51 tiny = 1.0e-18,
53 52 one = 1.0,
54 53 eight = 8.0,
55 54 invsqrtpi = 5.641895835477562869480794515607725858441e-0001,
56 55 tpi = 0.636619772367581343075535053490057448;
57 56
58 57 static GENERIC pzero(GENERIC), qzero(GENERIC);
59 58 static const GENERIC r0[4] = { /* [1.e-5, 1.28] */
60 59 -2.500000000000003622131880894830476755537e-0001,
61 60 1.095597547334830263234433855932375353303e-0002,
62 61 -1.819734750463320921799187258987098087697e-0004,
63 62 9.977001946806131657544212501069893930846e-0007,
64 63 };
65 64 static const GENERIC s0[4] = { /* [1.e-5, 1.28] */
66 65 1.0,
67 66 1.867609810662950169966782360588199673741e-0002,
68 67 1.590389206181565490878430827706972074208e-0004,
69 68 6.520867386742583632375520147714499522721e-0007,
70 69 };
71 70 static const GENERIC r1[9] = { /* [1.28,8] */
72 71 9.999999999999999942156495584397047660949e-0001,
73 72 -2.389887722731319130476839836908143731281e-0001,
74 73 1.293359476138939027791270393439493640570e-0002,
75 74 -2.770985642343140122168852400228563364082e-0004,
76 75 2.905241575772067678086738389169625218912e-0006,
77 76 -1.636846356264052597969042009265043251279e-0008,
78 77 5.072306160724884775085431059052611737827e-0011,
79 78 -8.187060730684066824228914775146536139112e-0014,
80 79 5.422219326959949863954297860723723423842e-0017,
81 80 };
82 81 static const GENERIC s1[9] = { /* [1.28,8] */
83 82 1.0,
84 83 1.101122772686807702762104741932076228349e-0002,
85 84 6.140169310641649223411427764669143978228e-0005,
86 85 2.292035877515152097976946119293215705250e-0007,
87 86 6.356910426504644334558832036362219583789e-0010,
88 87 1.366626326900219555045096999553948891401e-0012,
89 88 2.280399586866739522891837985560481180088e-0015,
90 89 2.801559820648939665270492520004836611187e-0018,
91 90 2.073101088320349159764410261466350732968e-0021,
92 91 };
93 92
94 93 GENERIC
95 94 j0(GENERIC x) {
96 95 GENERIC z, s, c, ss, cc, r, u, v, ox;
97 96 int i;
98 97
99 98 if (isnan(x))
100 99 return (x*x); /* + -> * for Cheetah */
101 100 ox = x;
102 101 x = fabs(x);
103 102 if (x > 8.0) {
104 103 if (!finite(x))
105 104 return (zero);
106 105 s = sin(x);
107 106 c = cos(x);
108 107 /*
109 108 * j0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x0)-q0(x)*sin(x0))
110 109 * where x0 = x-pi/4
111 110 * Better formula:
112 111 * cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4)
113 112 * = 1/sqrt(2) * (cos(x) + sin(x))
114 113 * sin(x0) = sin(x)cos(pi/4)-cos(x)sin(pi/4)
115 114 * = 1/sqrt(2) * (sin(x) - cos(x))
116 115 * To avoid cancellation, use
117 116 * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
118 117 * to compute the worse one.
119 118 */
120 119 if (x > 8.9e307) { /* x+x may overflow */
121 120 ss = s-c;
122 121 cc = s+c;
123 122 } else if (signbit(s) != signbit(c)) {
124 123 ss = s - c;
125 124 cc = -cos(x+x)/ss;
126 125 } else {
127 126 cc = s + c;
128 127 ss = -cos(x+x)/cc;
129 128 }
130 129 /*
131 130 * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x)
132 131 * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x)
133 132 */
134 133 if (x > 1.0e40) z = (invsqrtpi*cc)/sqrt(x);
135 134 else {
136 135 u = pzero(x); v = qzero(x);
137 136 z = invsqrtpi*(u*cc-v*ss)/sqrt(x);
138 137 }
139 138 /* force to pass SVR4 even the result is wrong (sign) */
140 139 if (x > X_TLOSS)
141 140 return (_SVID_libm_err(ox, z, 34));
142 141 else
143 142 return (z);
144 143 }
145 144 if (x <= small) {
146 145 if (x <= tiny)
147 146 return (one-x);
148 147 else
149 148 return (one-x*x*0.25);
150 149 }
151 150 z = x*x;
152 151 if (x <= 1.28) {
153 152 r = r0[0]+z*(r0[1]+z*(r0[2]+z*r0[3]));
154 153 s = s0[0]+z*(s0[1]+z*(s0[2]+z*s0[3]));
155 154 return (one + z*(r/s));
156 155 } else {
157 156 for (r = r1[8], s = s1[8], i = 7; i >= 0; i--) {
158 157 r = r*z + r1[i];
159 158 s = s*z + s1[i];
160 159 }
161 160 return (r/s);
162 161 }
163 162 }
164 163
165 164 static const GENERIC u0[13] = {
166 165 -7.380429510868722526754723020704317641941e-0002,
167 166 1.772607102684869924301459663049874294814e-0001,
168 167 -1.524370666542713828604078090970799356306e-0002,
169 168 4.650819100693891757143771557629924591915e-0004,
170 169 -7.125768872339528975036316108718239946022e-0006,
171 170 6.411017001656104598327565004771515257146e-0008,
172 171 -3.694275157433032553021246812379258781665e-0010,
173 172 1.434364544206266624252820889648445263842e-0012,
174 173 -3.852064731859936455895036286874139896861e-0015,
175 174 7.182052899726138381739945881914874579696e-0018,
176 175 -9.060556574619677567323741194079797987200e-0021,
177 176 7.124435467408860515265552217131230511455e-0024,
178 177 -2.709726774636397615328813121715432044771e-0027,
179 178 };
180 179 static const GENERIC v0[5] = {
181 180 1.0,
182 181 4.678678931512549002587702477349214886475e-0003,
183 182 9.486828955529948534822800829497565178985e-0006,
184 183 1.001495929158861646659010844136682454906e-0008,
185 184 4.725338116256021660204443235685358593611e-0012,
186 185 };
187 186
188 187 GENERIC
189 188 y0(GENERIC x) {
190 189 GENERIC z, /* d, */ s, c, ss, cc, u, v;
191 190 int i;
192 191
193 192 if (isnan(x))
194 193 return (x*x); /* + -> * for Cheetah */
195 194 if (x <= zero) {
196 195 if (x == zero)
197 196 /* d= -one/(x-x); */
198 197 return (_SVID_libm_err(x, x, 8));
199 198 else
200 199 /* d = zero/(x-x); */
201 200 return (_SVID_libm_err(x, x, 9));
202 201 }
203 202 if (x > 8.0) {
204 203 if (!finite(x))
205 204 return (zero);
206 205 s = sin(x);
207 206 c = cos(x);
208 207 /*
209 208 * j0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x0)-q0(x)*sin(x0))
210 209 * where x0 = x-pi/4
211 210 * Better formula:
212 211 * cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4)
213 212 * = 1/sqrt(2) * (cos(x) + sin(x))
214 213 * sin(x0) = sin(x)cos(pi/4)-cos(x)sin(pi/4)
215 214 * = 1/sqrt(2) * (sin(x) - cos(x))
216 215 * To avoid cancellation, use
217 216 * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
218 217 * to compute the worse one.
219 218 */
220 219 if (x > 8.9e307) { /* x+x may overflow */
221 220 ss = s-c;
222 221 cc = s+c;
223 222 } else if (signbit(s) != signbit(c)) {
224 223 ss = s - c;
225 224 cc = -cos(x+x)/ss;
226 225 } else {
227 226 cc = s + c;
228 227 ss = -cos(x+x)/cc;
229 228 }
230 229 /*
231 230 * j0(x) = 1/sqrt(pi*x) * (P(0,x)*cc - Q(0,x)*ss)
232 231 * y0(x) = 1/sqrt(pi*x) * (P(0,x)*ss + Q(0,x)*cc)
233 232 */
234 233 if (x > 1.0e40)
235 234 z = (invsqrtpi*ss)/sqrt(x);
236 235 else
237 236 z = invsqrtpi*(pzero(x)*ss+qzero(x)*cc)/sqrt(x);
238 237 if (x > X_TLOSS)
239 238 return (_SVID_libm_err(x, z, 35));
240 239 else
241 240 return (z);
242 241
243 242 }
244 243 if (x <= tiny) {
245 244 return (u0[0] + tpi*log(x));
246 245 }
247 246 z = x*x;
248 247 for (u = u0[12], i = 11; i >= 0; i--) u = u*z + u0[i];
249 248 v = v0[0]+z*(v0[1]+z*(v0[2]+z*(v0[3]+z*v0[4])));
250 249 return (u/v + tpi*(j0(x)*log(x)));
251 250 }
252 251
253 252 static const GENERIC pr[7] = { /* [8 -- inf] pzero 6550 */
254 253 .4861344183386052721391238447e5,
255 254 .1377662549407112278133438945e6,
256 255 .1222466364088289731869114004e6,
257 256 .4107070084315176135583353374e5,
258 257 .5026073801860637125889039915e4,
259 258 .1783193659125479654541542419e3,
260 259 .88010344055383421691677564e0,
261 260 };
262 261 static const GENERIC ps[7] = { /* [8 -- inf] pzero 6550 */
263 262 .4861344183386052721414037058e5,
264 263 .1378196632630384670477582699e6,
265 264 .1223967185341006542748936787e6,
266 265 .4120150243795353639995862617e5,
267 266 .5068271181053546392490184353e4,
268 267 .1829817905472769960535671664e3,
269 268 1.0,
270 269 };
271 270 static const GENERIC huge = 1.0e10;
272 271
273 272 static GENERIC
274 273 pzero(GENERIC x) {
275 274 GENERIC s, r, t, z;
276 275 int i;
277 276 if (x > huge)
278 277 return (one);
279 278 t = eight/x; z = t*t;
280 279 r = pr[5]+z*pr[6];
281 280 s = ps[5]+z;
282 281 for (i = 4; i >= 0; i--) {
283 282 r = r*z + pr[i];
284 283 s = s*z + ps[i];
285 284 }
286 285 return (r/s);
287 286 }
288 287
289 288 static const GENERIC qr[7] = { /* [8 -- inf] qzero 6950 */
290 289 -.1731210995701068539185611951e3,
291 290 -.5522559165936166961235240613e3,
292 291 -.5604935606637346590614529613e3,
293 292 -.2200430300226009379477365011e3,
294 293 -.323869355375648849771296746e2,
295 294 -.14294979207907956223499258e1,
296 295 -.834690374102384988158918e-2,
297 296 };
298 297 static const GENERIC qs[7] = { /* [8 -- inf] qzero 6950 */
299 298 .1107975037248683865326709645e5,
300 299 .3544581680627082674651471873e5,
301 300 .3619118937918394132179019059e5,
302 301 .1439895563565398007471485822e5,
303 302 .2190277023344363955930226234e4,
304 303 .106695157020407986137501682e3,
305 304 1.0,
306 305 };
307 306
308 307 static GENERIC
309 308 qzero(GENERIC x) {
310 309 GENERIC s, r, t, z;
311 310 int i;
312 311 if (x > huge)
313 312 return (-0.125/x);
314 313 t = eight/x; z = t*t;
315 314 r = qr[5]+z*qr[6];
316 315 s = qs[5]+z;
317 316 for (i = 4; i >= 0; i--) {
318 317 r = r*z + qr[i];
319 318 s = s*z + qs[i];
320 319 }
321 320 return (t*(r/s));
322 321 }
|
↓ open down ↓ |
269 lines elided |
↑ open up ↑ |
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX