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/j1.c
+++ new/usr/src/lib/libm/common/C/j1.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: j1(x),y1(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 j1 = __j1
40 -#pragma weak y1 = __y1
39 +#pragma weak __j1 = j1
40 +#pragma weak __y1 = y1
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-20,
53 52 one = 1.0,
54 53 invsqrtpi = 5.641895835477562869480794515607725858441e-0001,
55 54 tpi = 0.636619772367581343075535053490057448;
56 55
57 56 static GENERIC pone(GENERIC), qone(GENERIC);
58 57 static const GENERIC r0[4] = {
59 58 -6.250000000000002203053200981413218949548e-0002,
60 59 1.600998455640072901321605101981501263762e-0003,
61 60 -1.963888815948313758552511884390162864930e-0005,
62 61 8.263917341093549759781339713418201620998e-0008,
63 62 };
64 63 static const GENERIC s0[7] = {
65 64 1.0e0,
66 65 1.605069137643004242395356851797873766927e-0002,
67 66 1.149454623251299996428500249509098499383e-0004,
68 67 3.849701673735260970379681807910852327825e-0007,
69 68 };
70 69 static const GENERIC r1[12] = {
71 70 4.999999999999999995517408894340485471724e-0001,
72 71 -6.003825028120475684835384519945468075423e-0002,
73 72 2.301719899263321828388344461995355419832e-0003,
74 73 -4.208494869238892934859525221654040304068e-0005,
75 74 4.377745135188837783031540029700282443388e-0007,
76 75 -2.854106755678624335145364226735677754179e-0009,
77 76 1.234002865443952024332943901323798413689e-0011,
78 77 -3.645498437039791058951273508838177134310e-0014,
79 78 7.404320596071797459925377103787837414422e-0017,
80 79 -1.009457448277522275262808398517024439084e-0019,
81 80 8.520158355824819796968771418801019930585e-0023,
82 81 -3.458159926081163274483854614601091361424e-0026,
83 82 };
84 83 static const GENERIC s1[5] = {
85 84 1.0e0,
86 85 4.923499437590484879081138588998986303306e-0003,
87 86 1.054389489212184156499666953501976688452e-0005,
88 87 1.180768373106166527048240364872043816050e-0008,
89 88 5.942665743476099355323245707680648588540e-0012,
90 89 };
91 90
92 91 GENERIC
93 92 j1(GENERIC x) {
94 93 GENERIC z, d, s, c, ss, cc, r;
95 94 int i, sgn;
96 95
97 96 if (!finite(x))
98 97 return (one/x);
99 98 sgn = signbit(x);
100 99 x = fabs(x);
101 100 if (x > 8.00) {
102 101 s = sin(x);
103 102 c = cos(x);
104 103 /*
105 104 * j1(x) = sqrt(2/(pi*x))*(p1(x)*cos(x0)-q1(x)*sin(x0))
106 105 * where x0 = x-3pi/4
107 106 * Better formula:
108 107 * cos(x0) = cos(x)cos(3pi/4)+sin(x)sin(3pi/4)
109 108 * = 1/sqrt(2) * (sin(x) - cos(x))
110 109 * sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4)
111 110 * = -1/sqrt(2) * (cos(x) + sin(x))
112 111 * To avoid cancellation, use
113 112 * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
114 113 * to compute the worse one.
115 114 */
116 115 if (x > 8.9e307) { /* x+x may overflow */
117 116 ss = -s-c;
118 117 cc = s-c;
119 118 } else if (signbit(s) != signbit(c)) {
120 119 cc = s - c;
121 120 ss = cos(x+x)/cc;
122 121 } else {
123 122 ss = -s-c;
124 123 cc = cos(x+x)/ss;
125 124 }
126 125 /*
127 126 * j1(x) = 1/sqrt(pi*x) * (P(1,x)*cc - Q(1,x)*ss)
128 127 * y1(x) = 1/sqrt(pi*x) * (P(1,x)*ss + Q(1,x)*cc)
129 128 */
130 129 if (x > 1.0e40)
131 130 d = (invsqrtpi*cc)/sqrt(x);
132 131 else
133 132 d = invsqrtpi*(pone(x)*cc-qone(x)*ss)/sqrt(x);
134 133
135 134 if (x > X_TLOSS) {
136 135 if (sgn != 0) { d = -d; x = -x; }
137 136 return (_SVID_libm_err(x, d, 36));
138 137 } else
139 138 if (sgn == 0)
140 139 return (d);
141 140 else
142 141 return (-d);
143 142 }
144 143 if (x <= small) {
145 144 if (x <= tiny)
146 145 d = 0.5*x;
147 146 else
148 147 d = x*(0.5-x*x*0.125);
149 148 if (sgn == 0)
150 149 return (d);
151 150 else
152 151 return (-d);
153 152 }
154 153 z = x*x;
155 154 if (x < 1.28) {
156 155 r = r0[3];
157 156 s = s0[3];
158 157 for (i = 2; i >= 0; i--) {
159 158 r = r*z + r0[i];
160 159 s = s*z + s0[i];
161 160 }
162 161 d = x*0.5+x*(z*(r/s));
163 162 } else {
164 163 r = r1[11];
165 164 for (i = 10; i >= 0; i--) r = r*z + r1[i];
166 165 s = s1[0]+z*(s1[1]+z*(s1[2]+z*(s1[3]+z*s1[4])));
167 166 d = x*(r/s);
168 167 }
169 168 if (sgn == 0)
170 169 return (d);
171 170 else
172 171 return (-d);
173 172 }
174 173
175 174 static const GENERIC u0[4] = {
176 175 -1.960570906462389461018983259589655961560e-0001,
177 176 4.931824118350661953459180060007970291139e-0002,
178 177 -1.626975871565393656845930125424683008677e-0003,
179 178 1.359657517926394132692884168082224258360e-0005,
180 179 };
181 180 static const GENERIC v0[5] = {
182 181 1.0e0,
183 182 2.565807214838390835108224713630901653793e-0002,
184 183 3.374175208978404268650522752520906231508e-0004,
185 184 2.840368571306070719539936935220728843177e-0006,
186 185 1.396387402048998277638900944415752207592e-0008,
187 186 };
188 187 static const GENERIC u1[12] = {
189 188 -1.960570906462389473336339614647555351626e-0001,
190 189 5.336268030335074494231369159933012844735e-0002,
191 190 -2.684137504382748094149184541866332033280e-0003,
192 191 5.737671618979185736981543498580051903060e-0005,
193 192 -6.642696350686335339171171785557663224892e-0007,
194 193 4.692417922568160354012347591960362101664e-0009,
195 194 -2.161728635907789319335231338621412258355e-0011,
196 195 6.727353419738316107197644431844194668702e-0014,
197 196 -1.427502986803861372125234355906790573422e-0016,
198 197 2.020392498726806769468143219616642940371e-0019,
199 198 -1.761371948595104156753045457888272716340e-0022,
200 199 7.352828391941157905175042420249225115816e-0026,
201 200 };
202 201 static const GENERIC v1[5] = {
203 202 1.0e0,
204 203 5.029187436727947764916247076102283399442e-0003,
205 204 1.102693095808242775074856548927801750627e-0005,
206 205 1.268035774543174837829534603830227216291e-0008,
207 206 6.579416271766610825192542295821308730206e-0012,
208 207 };
209 208
210 209
211 210 GENERIC
212 211 y1(GENERIC x) {
213 212 GENERIC z, d, s, c, ss, cc, u, v;
214 213 int i;
215 214
216 215 if (isnan(x))
217 216 return (x*x); /* + -> * for Cheetah */
218 217 if (x <= zero) {
219 218 if (x == zero)
220 219 /* return -one/zero; */
221 220 return (_SVID_libm_err(x, x, 10));
222 221 else
223 222 /* return zero/zero; */
224 223 return (_SVID_libm_err(x, x, 11));
225 224 }
226 225 if (x > 8.0) {
227 226 if (!finite(x))
228 227 return (zero);
229 228 s = sin(x);
230 229 c = cos(x);
231 230 /*
232 231 * j1(x) = sqrt(2/(pi*x))*(p1(x)*cos(x0)-q1(x)*sin(x0))
233 232 * where x0 = x-3pi/4
234 233 * Better formula:
235 234 * cos(x0) = cos(x)cos(3pi/4)+sin(x)sin(3pi/4)
236 235 * = 1/sqrt(2) * (sin(x) - cos(x))
237 236 * sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4)
238 237 * = -1/sqrt(2) * (cos(x) + sin(x))
239 238 * To avoid cancellation, use
240 239 * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
241 240 * to compute the worse one.
242 241 */
243 242 if (x > 8.9e307) { /* x+x may overflow */
244 243 ss = -s-c;
245 244 cc = s-c;
246 245 } else if (signbit(s) != signbit(c)) {
247 246 cc = s - c;
248 247 ss = cos(x+x)/cc;
249 248 } else {
250 249 ss = -s-c;
251 250 cc = cos(x+x)/ss;
252 251 }
253 252 /*
254 253 * j1(x) = 1/sqrt(pi*x) * (P(1,x)*cc - Q(1,x)*ss)
255 254 * y1(x) = 1/sqrt(pi*x) * (P(1,x)*ss + Q(1,x)*cc)
256 255 */
257 256 if (x > 1.0e91)
258 257 d = (invsqrtpi*ss)/sqrt(x);
259 258 else
260 259 d = invsqrtpi*(pone(x)*ss+qone(x)*cc)/sqrt(x);
261 260
262 261 if (x > X_TLOSS)
263 262 return (_SVID_libm_err(x, d, 37));
264 263 else
265 264 return (d);
266 265 }
267 266 if (x <= tiny) {
268 267 return (-tpi/x);
269 268 }
270 269 z = x*x;
271 270 if (x < 1.28) {
272 271 u = u0[3]; v = v0[3]+z*v0[4];
273 272 for (i = 2; i >= 0; i--) {
274 273 u = u*z + u0[i];
275 274 v = v*z + v0[i];
276 275 }
277 276 } else {
278 277 for (u = u1[11], i = 10; i >= 0; i--) u = u*z+u1[i];
279 278 v = v1[0]+z*(v1[1]+z*(v1[2]+z*(v1[3]+z*v1[4])));
280 279 }
281 280 return (x*(u/v) + tpi*(j1(x)*log(x)-one/x));
282 281 }
283 282
284 283 static const GENERIC pr0[6] = {
285 284 -.4435757816794127857114720794e7,
286 285 -.9942246505077641195658377899e7,
287 286 -.6603373248364939109255245434e7,
288 287 -.1523529351181137383255105722e7,
289 288 -.1098240554345934672737413139e6,
290 289 -.1611616644324610116477412898e4,
291 290 };
292 291 static const GENERIC ps0[6] = {
293 292 -.4435757816794127856828016962e7,
294 293 -.9934124389934585658967556309e7,
295 294 -.6585339479723087072826915069e7,
296 295 -.1511809506634160881644546358e7,
297 296 -.1072638599110382011903063867e6,
298 297 -.1455009440190496182453565068e4,
299 298 };
300 299 static const GENERIC huge = 1.0e10;
301 300
302 301 static GENERIC
303 302 pone(GENERIC x) {
304 303 GENERIC s, r, t, z;
305 304 int i;
306 305 /* assume x > 8 */
307 306 if (x > huge)
308 307 return (one);
309 308
310 309 t = 8.0/x; z = t*t;
311 310 r = pr0[5]; s = ps0[5]+z;
312 311 for (i = 4; i >= 0; i--) {
313 312 r = z*r + pr0[i];
314 313 s = z*s + ps0[i];
315 314 }
316 315 return (r/s);
317 316 }
318 317
319 318
320 319 static const GENERIC qr0[6] = {
321 320 0.3322091340985722351859704442e5,
322 321 0.8514516067533570196555001171e5,
323 322 0.6617883658127083517939992166e5,
324 323 0.1849426287322386679652009819e5,
325 324 0.1706375429020768002061283546e4,
326 325 0.3526513384663603218592175580e2,
327 326 };
328 327 static const GENERIC qs0[6] = {
329 328 0.7087128194102874357377502472e6,
330 329 0.1819458042243997298924553839e7,
331 330 0.1419460669603720892855755253e7,
332 331 0.4002944358226697511708610813e6,
333 332 0.3789022974577220264142952256e5,
334 333 0.8638367769604990967475517183e3,
335 334 };
336 335
337 336 static GENERIC
338 337 qone(GENERIC x) {
339 338 GENERIC s, r, t, z;
340 339 int i;
341 340 if (x > huge)
342 341 return (0.375/x);
343 342
344 343 t = 8.0/x; z = t*t;
345 344 /* assume x > 8 */
346 345 r = qr0[5]; s = qs0[5]+z;
347 346 for (i = 4; i >= 0; i--) {
348 347 r = z*r + qr0[i];
349 348 s = z*s + qs0[i];
350 349 }
351 350 return (t*(r/s));
352 351 }
|
↓ open down ↓ |
299 lines elided |
↑ open up ↑ |
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX