5 * Common Development and Distribution License (the "License").
6 * You may not use this file except in compliance with the License.
7 *
8 * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
9 * or http://www.opensolaris.org/os/licensing.
10 * See the License for the specific language governing permissions
11 * and limitations under the License.
12 *
13 * When distributing Covered Code, include this CDDL HEADER in each
14 * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
15 * If applicable, add the following below this CDDL HEADER, with the
16 * fields enclosed by brackets "[]" replaced with your own identifying
17 * information: Portions Copyright [yyyy] [name of copyright owner]
18 *
19 * CDDL HEADER END
20 */
21
22 /*
23 * Copyright 2011 Nexenta Systems, Inc. All rights reserved.
24 */
25 /*
26 * Copyright 2006 Sun Microsystems, Inc. All rights reserved.
27 * Use is subject to license terms.
28 */
29
30 /*
31 * floating point Bessel's function of the first and second kinds
32 * of order zero: j1(x),y1(x);
33 *
34 * Special cases:
35 * y0(0)=y1(0)=yn(n,0) = -inf with division by zero signal;
36 * y0(-ve)=y1(-ve)=yn(n,-ve) are NaN with invalid signal.
37 */
38
39 #pragma weak __j1l = j1l
40 #pragma weak __y1l = y1l
41
42 #include "libm.h"
43 #include "longdouble.h"
44
45 #define GENERIC long double
46 static const GENERIC
47 zero = 0.0L,
48 small = 1.0e-9L,
49 tiny = 1.0e-38L,
50 one = 1.0L,
51 five = 5.0L,
52 invsqrtpi= 5.641895835477562869480794515607725858441e-0001L,
53 tpi = 0.636619772367581343075535053490057448L;
54
55 static GENERIC pone(GENERIC);
56 static GENERIC qone(GENERIC);
57
58 static const GENERIC r0[7] = {
59 -6.249999999999999999999999999999999627320e-0002L,
60 1.940606727194041716205384618494641565464e-0003L,
61 -3.005630423155733701856481469986459043883e-0005L,
62 2.345586219403918667468341047369572169358e-0007L,
63 -9.976809285885253587529010109133336669724e-0010L,
64 2.218743258363623946078958783775107473381e-0012L,
65 -2.071079656218700604767650924103578046280e-0015L,
66 };
67 static const GENERIC s0[7] = {
68 1.0e0L,
69 1.061695903156199920738051277075003059555e-0002L,
70 5.521860513111180371566951179398862692060e-0005L,
71 1.824214367413754193524107877084979441407e-0007L,
72 4.098957778439576834818838198039029353925e-0010L,
73 6.047735079699666389853240090925264056197e-0013L,
74 4.679044728878836197247923279512047035041e-0016L,
75 };
76
77 GENERIC
78 j1l(x) GENERIC x;{
79 GENERIC z, d, s,c,ss,cc,r;
80 int i, sgn;
81
82 if (!finitel(x)) return one/x;
83 sgn = signbitl(x);
84 x = fabsl(x);
85 if (x > 1.28L) {
86 s = sinl(x);
87 c = cosl(x);
88 /*
89 * j1(x) = sqrt(2/(pi*x))*(p1(x)*cos(x0)-q1(x)*sin(x0))
90 * where x0 = x-3pi/4
91 * Better formula:
92 * cos(x0) = cos(x)cos(3pi/4)+sin(x)sin(3pi/4)
93 * = 1/sqrt(2) * (sin(x) - cos(x))
94 * sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4)
95 * = -1/sqrt(2) * (cos(x) + sin(x))
96 * To avoid cancellation, use
97 * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
98 * to compute the worse one.
99 */
100 if (x>1.0e2450L) { /* x+x may overflow */
101 ss = -s-c;
102 cc = s-c;
103 } else if (signbitl(s) != signbitl(c)) {
104 cc = s - c;
105 ss = cosl(x+x)/cc;
106 } else {
107 ss = -s-c;
108 cc = cosl(x+x)/ss;
109 }
110 /*
111 * j1(x) = 1/sqrt(pi*x) * (P(1,x)*cc - Q(1,x)*ss)
112 * y1(x) = 1/sqrt(pi*x) * (P(1,x)*ss + Q(1,x)*cc)
113 */
114 if (x>1.0e120L) return (invsqrtpi*cc)/sqrtl(x);
115 d = invsqrtpi*(pone(x)*cc-qone(x)*ss)/sqrtl(x);
116 if (sgn == 0) return d; else return -d;
117 }
118 if (x<=small) {
119 if (x<=tiny) d = 0.5L*x;
120 else d = x*(0.5L-x*x*0.125L);
121 if (sgn == 0) return d; else return -d;
122 }
123 z = x*x;
124 r = r0[6];
125 s = s0[6];
126 for (i=5;i>=0;i--) {
127 r = r*z + r0[i];
128 s = s*z + s0[i];
129 }
130 d = x*0.5L+x*(z*(r/s));
131 if (sgn == 0) return d; else return -d;
132 }
133
134 static const GENERIC u0[7] = {
135 -1.960570906462389484060557273467558703503e-0001L,
136 5.166389353148318460304315890665450006495e-0002L,
137 -2.229699464105910913337190798743451115604e-0003L,
138 3.625437034548863342715657067759078267158e-0005L,
139 -2.689902826993117212255524537353883987171e-0007L,
140 9.304570592456930912969387719010256018466e-0010L,
141 -1.234878126794286643318321347997500346131e-0012L,
142 };
143 static const GENERIC v0[8] = {
144 1.0e0L,
145 1.369394302535807332517110204820556695644e-0002L,
146 9.508438148097659501433367062605935379588e-0005L,
147 4.399007309420092056052714797296467565655e-0007L,
148 1.488083087443756398305819693177715000787e-0009L,
149 3.751609832625793536245746965768587624922e-0012L,
150 6.680926434086257291872903276124244131448e-0015L,
151 6.676602383908906988160099057991121446058e-0018L,
152 };
153
154 GENERIC
155 y1l(x) GENERIC x;{
156 GENERIC z, s,c,ss,cc,u,v;
157 int i;
158
159 if (isnanl(x)) return x+x;
160 if (x <= zero) {
161 if (x == zero)
162 return -one/zero;
163 else
164 return zero/zero;
165 }
166 if (x > 1.28L) {
167 if (!finitel(x)) return zero;
168 s = sinl(x);
169 c = cosl(x);
170 /*
171 * j1(x) = sqrt(2/(pi*x))*(p1(x)*cos(x0)-q1(x)*sin(x0))
172 * where x0 = x-3pi/4
173 * Better formula:
174 * cos(x0) = cos(x)cos(3pi/4)+sin(x)sin(3pi/4)
175 * = 1/sqrt(2) * (sin(x) - cos(x))
176 * sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4)
177 * = -1/sqrt(2) * (cos(x) + sin(x))
178 * To avoid cancellation, use
179 * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
180 * to compute the worse one.
181 */
182 if (x>1.0e2450L) { /* x+x may overflow */
183 ss = -s-c;
184 cc = s-c;
185 } else if (signbitl(s) != signbitl(c)) {
186 cc = s - c;
187 ss = cosl(x+x)/cc;
188 } else {
189 ss = -s-c;
190 cc = cosl(x+x)/ss;
191 }
192 /*
193 * j1(x) = 1/sqrt(pi*x) * (P(1,x)*cc - Q(1,x)*ss)
194 * y1(x) = 1/sqrt(pi*x) * (P(1,x)*ss + Q(1,x)*cc)
195 */
196 if (x>1.0e91L) return (invsqrtpi*ss)/sqrtl(x);
197 return invsqrtpi*(pone(x)*ss+qone(x)*cc)/sqrtl(x);
198 }
199 if (x<=tiny) {
200 return (-tpi/x);
201 }
202 z = x*x;
203 u = u0[6]; v = v0[6]+z*v0[7];
204 for (i=5;i>=0;i--) {
205 u = u*z + u0[i];
206 v = v*z + v0[i];
207 }
208 return(x*(u/v) + tpi*(j1l(x)*logl(x)-one/x));
209 }
210
211 static const GENERIC pr0[12] = {
212 1.000000000000000000000000000000000000267e+0000L,
213 1.060717875045891455602180843276758003035e+0003L,
214 4.344347542892127024446687712181105852335e+0005L,
215 8.915680220724007016377924252717410457094e+0007L,
216 9.969502259938406062809873257569171272819e+0009L,
217 6.200290193138613035646510338707386316595e+0011L,
218 2.105978548788015119851815854422247330118e+0013L,
219 3.696635772784601239371730810311998368948e+0014L,
220 3.015913097920694682057958412534134515156e+0015L,
221 9.370298471339353098123277427328592725921e+0015L,
222 7.190349005196335967340799265074029443057e+0015L,
223 2.736097786240689996880391074927552517982e+0014L,
224 };
225 static const GENERIC ps0[11] = {
226 1.0e0L,
227 1.060600687545891455602180843276758095107e+0003L,
228 4.343106093416975589147153906505338900961e+0005L,
229 8.910605869002176566582072242244353399059e+0007L,
230 9.959122058635087888690713917622056540190e+0009L,
231 6.188744967234948231792482949171041843894e+0011L,
232 2.098863976953783506401759873801990304907e+0013L,
233 3.672870357018063196746729751479938908450e+0014L,
234 2.975538419246824921049011529574385888420e+0015L,
235 9.063657659995043205018686029284479837091e+0015L,
236 6.401953344314747916729366441508892711691e+0015L,
237 };
238 static const GENERIC pr1[12] = {
239 1.000000000000000000000023667524130660984e+0000L,
240 6.746154419979618754354803488126452971204e+0002L,
241 1.811210781083390154857018330296145970502e+0005L,
242 2.533098390379924268038005329095287842244e+0007L,
243 2.029683619805342145252338570875424600729e+0009L,
244 9.660859662192711465301069401598929980319e+0010L,
245 2.743396238644831519934098967716621316316e+0012L,
246 4.553097354140854377931023170263455246288e+0013L,
247 4.210245069852219757476169864974870720374e+0014L,
248 1.987334056229596485076645967176169801727e+0015L,
249 4.067120052787096893838970455751338930462e+0015L,
250 2.486539606380406398310845264910691398133e+0015L,
251 };
252 static const GENERIC ps1[14] = {
253 1.0e0L,
254 6.744982544979618754355808680196859521782e+0002L,
255 1.810421795396966762032155290441364740350e+0005L,
256 2.530986460644310651529583759699988435573e+0007L,
257 2.026743276048023121360249288818290224145e+0009L,
258 9.637461924407405935245269407052641341836e+0010L,
259 2.732378628423766417402292797028314160831e+0012L,
260 4.522345274960527124354844364012184278488e+0013L,
261 4.160650668341743132685335758415469856545e+0014L,
262 1.943730242988858208243492424892435901211e+0015L,
263 3.880228532692127989901131618598067450001e+0015L,
264 2.178020816161154615841000173683302999728e+0015L,
265 -8.994062666842225551554346698171600634173e+0013L,
266 1.368520368508851253495764806934619574990e+0013L,
267 };
268 static const GENERIC pr2[12] = {
269 1.000000000000000006938651621840396237282e+0000L,
270 3.658416291850404981407101077037948144698e+0002L,
271 5.267073772170356547709794670602812447537e+0004L,
272 3.912012101226837463014925210735894620442e+0006L,
273 1.651295648974103957193874928714180765625e+0008L,
274 4.114901144480797609972484998142146783499e+0009L,
275 6.092524309766036681542980572526335147672e+0010L,
276 5.263913178071282616719249969074134570577e+0011L,
277 2.538408581124324223367341020538081330994e+0012L,
278 6.288607929360291027895126983015365677648e+0012L,
279 6.848330048211148419047055075386525945280e+0012L,
280 2.290309646838867941423178163991423244690e+0012L,
281 };
282 static const GENERIC ps2[14] = {
283 1.0e0L,
284 3.657244416850405086459410165762319861856e+0002L,
285 5.262802358425023243992387075861237306312e+0004L,
286 3.905896813959919648136295861661483848364e+0006L,
287 1.646791907791461220742694842108202772763e+0008L,
288 4.096132803064256022224954120208201437344e+0009L,
289 6.046665195915950447544429445730680236759e+0010L,
290 5.198061739781991313414052212328653295168e+0011L,
291 2.484233851814333966401527626421254279796e+0012L,
292 6.047868806925315879339651539434315255940e+0012L,
293 6.333103831254091652501642567294101813354e+0012L,
294 1.875143098754284994467609936924685024968e+0012L,
295 -5.238330920563392692965412762508813601534e+0010L,
296 4.656888609439364725427789198383779259957e+0009L,
297 };
298 static const GENERIC pr3[13] = {
299 1.000000000000009336887318068056137842897e+0000L,
300 2.242719942728459588488051572002835729183e+0002L,
301 1.955450611382026550266257737331095691092e+0004L,
302 8.707143293993619899395400562409175590739e+0005L,
303 2.186267894487004565948324289010954505316e+0007L,
304 3.224328510541957792360691585667502864688e+0008L,
305 2.821057355151380597331792896882741364897e+0009L,
306 1.445371387295422404365584793796028979840e+0010L,
307 4.181743160669891357783011002656658107864e+0010L,
308 6.387371088767993119325536137794535513922e+0010L,
309 4.575619999412716078064070587767416436396e+0010L,
310 1.228415651211639160620284441690503550842e+0010L,
311 7.242170349875563053436050532153112882072e+0008L,
312 };
313 static const GENERIC ps3[13] = {
314 1.0e0L,
315 2.241548067728529551049804610486061401070e+0002L,
316 1.952838216795552145132137932931237181307e+0004L,
317 8.684574926493185744628127341069974575526e+0005L,
318 2.176357771067037962940853412819852189164e+0007L,
319 3.199958682356132977319258783167122100567e+0008L,
320 2.786218931525334687844675219914201872570e+0009L,
321 1.416283776951741549631417572317916039767e+0010L,
322 4.042962659271567948735676834609348842922e+0010L,
323 6.028168462646694510083847222968444402161e+0010L,
324 4.118410226794641413833887606580085281111e+0010L,
325 9.918735736297038430744161253338202230263e+0009L,
326 4.092967198238098023219124487437130332038e+0008L,
327 };
328 static const GENERIC pr4[13] = {
329 1.000000000001509220978157399042059553390e+0000L,
330 1.437551868378147851133499996323782607787e+0002L,
331 7.911335537418177296041518061404505428004e+0003L,
332 2.193710939115317214716518908935756104804e+0005L,
333 3.390662495136730962513489796538274984382e+0006L,
334 3.048655347929348891006070609293884274789e+0007L,
335 1.613781633489496606354045161527450975195e+0008L,
336 4.975089835037230277110156150038482159988e+0008L,
337 8.636047087015115403880904418339566323264e+0008L,
338 7.918202912328366140110671223076949101509e+0008L,
339 3.423294665798984733439650311722794853294e+0008L,
340 5.621904953441963961040503934782662613621e+0007L,
341 2.086303543310240260758670404509484499793e+0006L,
342 };
343 static const GENERIC ps4[13] = {
344 1.0e0L,
345 1.436379993384532371670493319591847362304e+0002L,
346 7.894647154785430678061053848847436659499e+0003L,
347 2.184659753392097529008981741550878586174e+0005L,
348 3.366109083305465176803513738147049499361e+0006L,
349 3.011911545968996817697665866587226343186e+0007L,
350 1.582262913779689851316760148459414895301e+0008L,
351 4.819268809494937919217938589530138201770e+0008L,
352 8.201355762990450679702837123432527154830e+0008L,
353 7.268232093982510937417446421282341425212e+0008L,
354 2.950911909015572933262131323934036480462e+0008L,
355 4.242839924305934423010858966540621219396e+0007L,
356 1.064387620445090779182117666330405186866e+0006L,
357 };
358 static const GENERIC pr5[13] = {
359 1.000000000102434805241171427253847353861e+0000L,
360 9.129332257083629259060502249025963234821e+0001L,
361 3.132238483586953037576119377504557191413e+0003L,
362 5.329782528269307971278943122454171107861e+0004L,
363 4.988460157184117790692873002103052944145e+0005L,
364 2.686602071615786816147010334256047469378e+0006L,
365 8.445418526028961197703799808701268301831e+0006L,
366 1.536575358646141157475725889907900827390e+0007L,
367 1.568405818236523821796862770586544811945e+0007L,
368 8.450876239888770102387618667362302173547e+0006L,
369 2.154414900139567328424026827163203446077e+0006L,
370 2.105656926565043898888460254808062352205e+0005L,
371 4.739165011023396507022134303736862812975e+0003L,
372 };
373 static const GENERIC ps5[13] = {
374 1.0e0L,
375 9.117613509595327476509152673394703847793e+0001L,
376 3.121697972484015639301279229281770795147e+0003L,
377 5.294447222735893568040911873834576440255e+0004L,
378 4.930368882192772335798256684110887882807e+0005L,
379 2.634854685641165298302167435798357437768e+0006L,
380 8.185462775400326393555896157031818280918e+0006L,
381 1.462417423080215192609668642663030667086e+0007L,
382 1.450624993985851675982860844153954896015e+0007L,
383 7.460467647561995283219086567162006113864e+0006L,
384 1.754210981405612478869227142579056338965e+0006L,
385 1.463286721155271971526264914524746699596e+0005L,
386 2.155894725796702015341211116579827039459e+0003L,
387 };
388 static const GENERIC pr6[13] = {
389 1.000000003564855546741735920315743157129e+0000L,
390 5.734003934862540458119423509909510288366e+0001L,
391 1.209572491935850486086559692291796887976e+0003L,
392 1.243398391422281247933674779163660286838e+0004L,
393 6.930996755181437937258220998601708278787e+0004L,
394 2.198067659532757598646722249966767620099e+0005L,
395 4.033659432712058633933179115820576858455e+0005L,
396 4.257759657219008027016047206574574358678e+0005L,
397 2.511917395876004349480721277445763916389e+0005L,
398 7.813756153070623654178731651381881953552e+0004L,
399 1.152069173381127881385588092905864352891e+0004L,
400 6.548580782804088553777816037551523398082e+0002L,
401 8.668725370116906132327542766127938496880e+0000L,
402 };
403 static const GENERIC ps6[13] = {
404 1.0e0L,
405 5.722285236357114566499221525736286205184e+0001L,
406 1.203010842878317935444582950620339570506e+0003L,
407 1.230058335378583550155825502172435371208e+0004L,
408 6.800998550607861288865300438648089894412e+0004L,
409 2.130767829599304262987769347536850885921e+0005L,
410 3.840483466643916681759936972992155310026e+0005L,
411 3.947432373459225542861819148108081160393e+0005L,
412 2.237816239393081111481588434457838526738e+0005L,
413 6.545820495124419723398946273790921540774e+0004L,
414 8.729563630320892741500726213278834737196e+0003L,
415 4.130762660291894753450174794196998813709e+0002L,
416 3.480368898672684645130335786015075595598e+0000L,
417 };
418 static const GENERIC sixteen = 16.0L;
419 static const GENERIC eight = 8.0L;
420 static const GENERIC huge = 1.0e30L;
421
422 static GENERIC pone(x)
423 GENERIC x;
424 {
425 GENERIC s,r,t,z;
426 int i;
427 if (x>huge) return one;
428 t = one/x; z = t*t;
429 if (x>sixteen) {
430 r = z*pr0[11]+pr0[10]; s = ps0[10];
431 for (i=9;i>=0;i--) {
432 r = z*r + pr0[i];
433 s = z*s + ps0[i];
434 }
435 } else if (x>eight) {
436 r = pr1[11]; s = ps1[11]+z*(ps1[12]+z*ps1[13]);
437 for (i=10;i>=0;i--) {
438 r = z*r + pr1[i];
439 s = z*s + ps1[i];
440 }
441 } else if (x>five) {
442 r = pr2[11]; s = ps2[11]+z*(ps2[12]+z*ps2[13]);
443 for (i=10;i>=0;i--) {
444 r = z*r + pr2[i];
445 s = z*s + ps2[i];
446 }
447 } else if (x>3.5L) {
448 r = pr3[12]; s = ps3[12];
449 for (i=11;i>=0;i--) {
450 r = z*r + pr3[i];
451 s = z*s + ps3[i];
452 }
453 } else if (x>2.5L) {
454 r = pr4[12]; s = ps4[12];
455 for (i=11;i>=0;i--) {
456 r = z*r + pr4[i];
457 s = z*s + ps4[i];
458 }
459 } else if (x> (1.0L/0.5625L)) {
460 r = pr5[12]; s = ps5[12];
461 for (i=11;i>=0;i--) {
462 r = z*r + pr5[i];
463 s = z*s + ps5[i];
464 }
465 } else { /* assume x > 1.28 */
466 r = pr6[12]; s = ps6[12];
467 for (i=11;i>=0;i--) {
468 r = z*r + pr6[i];
469 s = z*s + ps6[i];
470 }
471 }
472 return r/s;
473 }
474
475
476 static const GENERIC qr0[12] = {
477 3.749999999999999999999999999999999971033e-0001L,
478 4.256726035237050601607682277433094262226e+0002L,
479 1.875976490812878489192409978945401066066e+0005L,
480 4.170314268048041914273603680317745592790e+0007L,
481 5.092750132543855817293451118974555746551e+0009L,
482 3.494749676278488654103505795794139483404e+0011L,
483 1.327062148257437316997667817096694173709e+0013L,
484 2.648993138273427226907503742066551150490e+0014L,
485 2.511695665909547412222430494473998127684e+0015L,
486 9.274694506662289043224310499164702306096e+0015L,
487 8.150904170663663829331320302911792892002e+0015L,
488 -5.001918733707662355772037829620388765122e+0014L,
489 };
490 static const GENERIC qs0[11] = {
491 1.0e0L,
492 1.135400380229880160428715273982155760093e+0003L,
493 5.005701183877126164326765545516590744360e+0005L,
494 1.113444200113712167984337603933040102987e+0008L,
495 1.361074819925223062778717565699039471124e+0010L,
496 9.355750985802849484438933905325982809653e+0011L,
497 3.563462786008988825003965543857998084828e+0013L,
498 7.155145113900094163648726863803802910454e+0014L,
499 6.871266835834472758055559013851843654113e+0015L,
500 2.622030899226736712644974988157345234092e+0016L,
501 2.602912729172876330650077021706139707746e+0016L,
502 };
503 static const GENERIC qr1[12] = {
504 3.749999999999999999997762458207284405806e-0001L,
505 2.697883998881706839929255517498189980485e+0002L,
506 7.755195925781028489386938870473834411019e+0004L,
507 1.166777762104017777198211072895528968355e+0007L,
508 1.011504772984321168320010084520261069362e+0009L,
509 5.246007703574156853577754571720205550010e+0010L,
510 1.637692549885592683166116551691266537647e+0012L,
511 3.022303623698185669912990310925039382495e+0013L,
512 3.154769927290655684846107030265909987946e+0014L,
513 1.715819913441554770089730934808123360921e+0015L,
514 4.165044355759732622273534445131736188510e+0015L,
515 3.151381420874174705643100381708086287596e+0015L,
516 };
517 static const GENERIC qs1[14] = {
518 1.0e0L,
519 7.197091705351218239785633172408276982828e+0002L,
520 2.070012799599548685544883041297609861055e+0005L,
521 3.117014815317656221871840152778458754516e+0007L,
522 2.705719678902554974863325877025902971727e+0009L,
523 1.406113614727345726925060648750867264098e+0011L,
524 4.403777536067131320363005978631674817359e+0012L,
525 8.170725690209322283061499386703167242894e+0013L,
526 8.609458844975495289227794126964431210566e+0014L,
527 4.766766367015473481257280600694952920204e+0015L,
528 1.202286587943342194863557940888115641650e+0016L,
529 1.012474328306200909525063936061756024120e+0016L,
530 6.183552022678917858273222879615824070703e+0014L,
531 -9.756731548558226997573737400988225722740e+0013L,
532 };
533 static const GENERIC qr2[12] = {
534 3.749999999999999481245647262226994293189e-0001L,
535 1.471366807289771354491181140167359026735e+0002L,
536 2.279432486768448220142080962843526951250e+0004L,
537 1.828943048523771225163804043356958285893e+0006L,
538 8.379828388647823135832220596417725010837e+0007L,
539 2.279814029335044024585393671278378022053e+0009L,
540 3.711653952257118120832817785271466441420e+0010L,
541 3.557650914518554549916730572553105048068e+0011L,
542 1.924583483146095896259774329498934160650e+0012L,
543 5.424386256063736390759567088291887140278e+0012L,
544 6.839325621241776786206509704671746841737e+0012L,
545 2.702169563144001166291686452305436313971e+0012L,
546 };
547 static const GENERIC qs2[14] = {
548 1.0e0L,
549 3.926379194439388135703211933895203191089e+0002L,
550 6.089148804106598297488336063007609312276e+0004L,
551 4.893546162973278583711376356041614150645e+0006L,
552 2.247571119114497845046388801813832219404e+0008L,
553 6.137635663350177751290469334200757872645e+0009L,
554 1.005115019784102856424493519524998953678e+0011L,
555 9.725664462014503832860151384604677240620e+0011L,
556 5.345525100485511116148634192844434636072e+0012L,
557 1.549944007398946691720862738173956994779e+0013L,
558 2.067148441178952625710302124163264760362e+0013L,
559 9.401565402641963611295119487242595462301e+0012L,
560 3.548217088622398274748837287769709374385e+0011L,
561 -2.934470341719047120076509938432417352365e+0010L,
562 };
563 static const GENERIC qr3[13] = {
564 3.749999999999412724084579833297451472091e-0001L,
565 9.058478580291706212422978492938435582527e+0001L,
566 8.524056033161038750461083666711724381171e+0003L,
567 4.105967158629109427753434569223631014730e+0005L,
568 1.118326603378531348259783091972623333657e+0007L,
569 1.794636683403578918528064904714132329343e+0008L,
570 1.714314157463635959556133236004368896724e+0009L,
571 9.622092032236084846572067257267661456030e+0009L,
572 3.057759524485859159957762858780768355020e+0010L,
573 5.129306780754798531609621454415938890020e+0010L,
574 3.999122002794961070680636194346316041352e+0010L,
575 1.122298454643493485989721564358100345388e+0010L,
576 5.603981987645989709668830968522362582221e+0008L,
577 };
578 static const GENERIC qs3[13] = {
579 1.0e0L,
580 2.418328663076578169836155170053634419922e+0002L,
581 2.279620205900121042587523541281272875520e+0004L,
582 1.100984222585729521470129014992217092794e+0006L,
583 3.010743223679247091004262516286654516282e+0007L,
584 4.860925542827367817289619265215599433996e+0008L,
585 4.686668111035348691982715864307839581243e+0009L,
586 2.668701788405102017427214705946730894074e+0010L,
587 8.677395746106802640390580944836650584903e+0010L,
588 1.511936455574951790658498795945106643036e+0011L,
589 1.260845604432623478002018696873608353093e+0011L,
590 4.052692278419853853911440231600864589805e+0010L,
591 2.965516519212226064983267822243329694729e+0009L,
592 };
593 static const GENERIC qr4[13] = {
594 3.749999999919234164154669754440123072618e-0001L,
595 5.844218580776819864791168253485055101858e+0001L,
596 3.489273514092912982675669411371435670220e+0003L,
597 1.050523637774575684509663430018995479594e+0005L,
598 1.764549172059701565500717319792780115289e+0006L,
599 1.725532438844133795028063102681497371154e+0007L,
600 9.938114847359778539965140247590176334874e+0007L,
601 3.331710808184595545396883770200772842314e+0008L,
602 6.271970557641881511609560444872797282698e+0008L,
603 6.188529798677357075020774923903737913285e+0008L,
604 2.821905302742849974509982167877885011629e+0008L,
605 4.615467358646911976773290256984329814896e+0007L,
606 1.348140608731546467396685802693380693275e+0006L,
607 };
608 static const GENERIC qs4[13] = {
609 1.0e0L,
610 1.561192663112345185261418296389902133372e+0002L,
611 9.346678031144098270547225423124213083072e+0003L,
612 2.825851246482293547838023847601704751590e+0005L,
613 4.776572711622156091710902891124911556293e+0006L,
614 4.715106953717135402977938048006267859302e+0007L,
615 2.753962350894311316439652227611209035193e+0008L,
616 9.428501434615463207768964787500411575223e+0008L,
617 1.832650858775206787088236896454141572617e+0009L,
618 1.901697378939743226948920874296595242257e+0009L,
619 9.433322226854293780627188599226380812725e+0008L,
620 1.808520540608671608680284520798858587370e+0008L,
621 7.983342331736662753157217446919462398008e+0006L,
622 };
623 static const GENERIC qr5[13] = {
624 3.749999995331364437028988850515190446719e-0001L,
625 3.739356381766559882677514593041627547911e+0001L,
626 1.399562500629413529355265462912819802551e+0003L,
627 2.594154053098947925345332218062210111753e+0004L,
628 2.640149879297408640394163979394594318371e+0005L,
629 1.542471854873199142031889093591449397995e+0006L,
630 5.242272868972053374067572098992335425895e+0006L,
631 1.025834487769410221329633071426044839935e+0007L,
632 1.116553924239448940142230579060124209622e+0007L,
633 6.318076065595910176374916303525884653514e+0006L,
634 1.641218086168640408527639735915512881785e+0006L,
635 1.522369793529178644168813882912134706444e+0005L,
636 2.526530541062297200914180060208669584055e+0003L,
637 };
638 static const GENERIC qs5[13] = {
639 1.0e0L,
640 9.998960735935075380397545659016287506660e+0001L,
641 3.758767417842043742686475060540416737562e+0003L,
642 7.013652806952306520121959742519780781653e+0004L,
643 7.208949808818615099246529616211730446850e+0005L,
644 4.272753927109614455417836186072202009252e+0006L,
645 1.482524411356470699336129814111025434703e+0007L,
646 2.988750366665678233425279237627700803473e+0007L,
647 3.396957890261080492694709150553619185065e+0007L,
648 2.050652487738593004111578091156304540386e+0007L,
649 5.900504120811732547616511555946279451316e+0006L,
650 6.563391409260160897024498082273183468347e+0005L,
651 1.692629845012790205348966731477187041419e+0004L,
652 };
653 static const GENERIC qr6[13] = {
654 3.749999861516664133157566870858975421296e-0001L,
655 2.367863756747764863120797431599473468918e+0001L,
656 5.476715802114976248882067325630793143777e+0002L,
657 6.143190357869842894025012945444096170251e+0003L,
658 3.716250534677997850513733595140463851730e+0004L,
659 1.270883463823876752138326905022875657430e+0005L,
660 2.495301449636814481646371665429083801388e+0005L,
661 2.789578988212952248340486296254398601942e+0005L,
662 1.718247946911109055931819087137397324634e+0005L,
663 5.458973214011665714330326732204106364229e+0004L,
664 7.912102686687948786048943339759596652813e+0003L,
665 4.077961006160866935722030715149087938091e+0002L,
666 3.765206972770245085551057237882528510428e+0000L,
667 };
668 static const GENERIC qs6[13] = {
669 1.0e0L,
670 6.341646532940517305641893852673926809601e+0001L,
671 1.477058277414040790932597537920671025359e+0003L,
672 1.674406564031044491436044253393536487604e+0004L,
673 1.028516501369755949895050806908994650768e+0005L,
674 3.593620042532885295087463507733285434207e+0005L,
675 7.267924991381020915185873399453724799625e+0005L,
676 8.462277510768818399961191426205006083088e+0005L,
677 5.514399892230892163373611895645500250514e+0005L,
678 1.898084241009259353540620272932188102299e+0005L,
679 3.102941242117739015721984123081026253068e+0004L,
680 1.958971184431466907681440650181421086143e+0003L,
681 2.878853357310495087181721613889455121867e+0001L,
682 };
683 static GENERIC qone(x)
684 GENERIC x;
685 {
686 GENERIC s,r,t,z;
687 int i;
688 if (x>huge) return 0.375L/x;
689 t = one/x; z = t*t;
690 if (x>sixteen) {
691 r = z*qr0[11]+qr0[10]; s = qs0[10];
692 for (i=9;i>=0;i--) {
693 r = z*r + qr0[i];
694 s = z*s + qs0[i];
695 }
696 } else if (x>eight) {
697 r = qr1[11]; s = qs1[11]+z*(qs1[12]+z*qs1[13]);
698 for (i=10;i>=0;i--) {
699 r = z*r + qr1[i];
700 s = z*s + qs1[i];
701 }
702 } else if (x>five) { /* x > 5.0 */
703 r = qr2[11]; s = qs2[11]+z*(qs2[12]+z*qs2[13]);
704 for (i=10;i>=0;i--) {
705 r = z*r + qr2[i];
706 s = z*s + qs2[i];
707 }
708 } else if (x>3.5L) {
709 r = qr3[12]; s = qs3[12];
710 for (i=11;i>=0;i--) {
711 r = z*r + qr3[i];
712 s = z*s + qs3[i];
713 }
714 } else if (x>2.5L) {
715 r = qr4[12]; s = qs4[12];
716 for (i=11;i>=0;i--) {
717 r = z*r + qr4[i];
718 s = z*s + qs4[i];
719 }
720 } else if (x> (1.0L/0.5625L)) {
721 r = qr5[12]; s = qs5[12];
722 for (i=11;i>=0;i--) {
723 r = z*r + qr5[i];
724 s = z*s + qs5[i];
725 }
726 } else { /* assume x > 1.28 */
727 r = qr6[12]; s = qs6[12];
728 for (i=11;i>=0;i--) {
729 r = z*r + qr6[i];
730 s = z*s + qs6[i];
731 }
732 }
733 return t*(r/s);
734 }
|
5 * Common Development and Distribution License (the "License").
6 * You may not use this file except in compliance with the License.
7 *
8 * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
9 * or http://www.opensolaris.org/os/licensing.
10 * See the License for the specific language governing permissions
11 * and limitations under the License.
12 *
13 * When distributing Covered Code, include this CDDL HEADER in each
14 * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
15 * If applicable, add the following below this CDDL HEADER, with the
16 * fields enclosed by brackets "[]" replaced with your own identifying
17 * information: Portions Copyright [yyyy] [name of copyright owner]
18 *
19 * CDDL HEADER END
20 */
21
22 /*
23 * Copyright 2011 Nexenta Systems, Inc. All rights reserved.
24 */
25
26 /*
27 * Copyright 2006 Sun Microsystems, Inc. All rights reserved.
28 * Use is subject to license terms.
29 */
30
31 /*
32 * floating point Bessel's function of the first and second kinds
33 * of order zero: j1(x),y1(x);
34 *
35 * Special cases:
36 * y0(0)=y1(0)=yn(n,0) = -inf with division by zero signal;
37 * y0(-ve)=y1(-ve)=yn(n,-ve) are NaN with invalid signal.
38 */
39
40 #pragma weak __j1l = j1l
41 #pragma weak __y1l = y1l
42
43 #include "libm.h"
44 #include "longdouble.h"
45
46 #define GENERIC long double
47
48 static const GENERIC zero = 0.0L,
49 small = 1.0e-9L,
50 tiny = 1.0e-38L,
51 one = 1.0L,
52 five = 5.0L,
53 invsqrtpi = 5.641895835477562869480794515607725858441e-0001L,
54 tpi = 0.636619772367581343075535053490057448L;
55
56 static GENERIC pone(GENERIC);
57 static GENERIC qone(GENERIC);
58
59 static const GENERIC r0[7] = {
60 -6.249999999999999999999999999999999627320e-0002L,
61 1.940606727194041716205384618494641565464e-0003L,
62 -3.005630423155733701856481469986459043883e-0005L,
63 2.345586219403918667468341047369572169358e-0007L,
64 -9.976809285885253587529010109133336669724e-0010L,
65 2.218743258363623946078958783775107473381e-0012L,
66 -2.071079656218700604767650924103578046280e-0015L,
67 };
68
69 static const GENERIC s0[7] = {
70 1.0e0L,
71 1.061695903156199920738051277075003059555e-0002L,
72 5.521860513111180371566951179398862692060e-0005L,
73 1.824214367413754193524107877084979441407e-0007L,
74 4.098957778439576834818838198039029353925e-0010L,
75 6.047735079699666389853240090925264056197e-0013L,
76 4.679044728878836197247923279512047035041e-0016L,
77 };
78
79 GENERIC
80 j1l(GENERIC x)
81 {
82 GENERIC z, d, s, c, ss, cc, r;
83 int i, sgn;
84
85 if (!finitel(x))
86 return (one / x);
87
88 sgn = signbitl(x);
89 x = fabsl(x);
90
91 if (x > 1.28L) {
92 s = sinl(x);
93 c = cosl(x);
94
95 /* BEGIN CSTYLED */
96 /*
97 * j1(x) = sqrt(2/(pi*x))*(p1(x)*cos(x0)-q1(x)*sin(x0))
98 * where x0 = x-3pi/4
99 * Better formula:
100 * cos(x0) = cos(x)cos(3pi/4)+sin(x)sin(3pi/4)
101 * = 1/sqrt(2) * (sin(x) - cos(x))
102 * sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4)
103 * = -1/sqrt(2) * (cos(x) + sin(x))
104 * To avoid cancellation, use
105 * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
106 * to compute the worse one.
107 */
108 /* END CSTYLED */
109 if (x > 1.0e2450L) { /* x+x may overflow */
110 ss = -s - c;
111 cc = s - c;
112 } else if (signbitl(s) != signbitl(c)) {
113 cc = s - c;
114 ss = cosl(x + x) / cc;
115 } else {
116 ss = -s - c;
117 cc = cosl(x + x) / ss;
118 }
119
120 /*
121 * j1(x) = 1/sqrt(pi*x) * (P(1,x)*cc - Q(1,x)*ss)
122 * y1(x) = 1/sqrt(pi*x) * (P(1,x)*ss + Q(1,x)*cc)
123 */
124 if (x > 1.0e120L)
125 return ((invsqrtpi * cc) / sqrtl(x));
126
127 d = invsqrtpi * (pone(x) * cc - qone(x) * ss) / sqrtl(x);
128
129 if (sgn == 0)
130 return (d);
131 else
132 return (-d);
133 }
134
135 if (x <= small) {
136 if (x <= tiny)
137 d = 0.5L * x;
138 else
139 d = x * (0.5L - x * x * 0.125L);
140
141 if (sgn == 0)
142 return (d);
143 else
144 return (-d);
145 }
146
147 z = x * x;
148 r = r0[6];
149 s = s0[6];
150
151 for (i = 5; i >= 0; i--) {
152 r = r * z + r0[i];
153 s = s * z + s0[i];
154 }
155
156 d = x * 0.5L + x * (z * (r / s));
157
158 if (sgn == 0)
159 return (d);
160 else
161 return (-d);
162 }
163
164 static const GENERIC u0[7] = {
165 -1.960570906462389484060557273467558703503e-0001L,
166 5.166389353148318460304315890665450006495e-0002L,
167 -2.229699464105910913337190798743451115604e-0003L,
168 3.625437034548863342715657067759078267158e-0005L,
169 -2.689902826993117212255524537353883987171e-0007L,
170 9.304570592456930912969387719010256018466e-0010L,
171 -1.234878126794286643318321347997500346131e-0012L,
172 };
173
174 static const GENERIC v0[8] = {
175 1.0e0L,
176 1.369394302535807332517110204820556695644e-0002L,
177 9.508438148097659501433367062605935379588e-0005L,
178 4.399007309420092056052714797296467565655e-0007L,
179 1.488083087443756398305819693177715000787e-0009L,
180 3.751609832625793536245746965768587624922e-0012L,
181 6.680926434086257291872903276124244131448e-0015L,
182 6.676602383908906988160099057991121446058e-0018L,
183 };
184
185 GENERIC
186 y1l(GENERIC x)
187 {
188 GENERIC z, s, c, ss, cc, u, v;
189 int i;
190
191 if (isnanl(x))
192 return (x + x);
193
194 if (x <= zero) {
195 if (x == zero)
196 return (-one / zero);
197 else
198 return (zero / zero);
199 }
200
201 if (x > 1.28L) {
202 if (!finitel(x))
203 return (zero);
204
205 s = sinl(x);
206 c = cosl(x);
207
208 /* BEGIN CSTYLED */
209 /*
210 * j1(x) = sqrt(2/(pi*x))*(p1(x)*cos(x0)-q1(x)*sin(x0))
211 * where x0 = x-3pi/4
212 * Better formula:
213 * cos(x0) = cos(x)cos(3pi/4)+sin(x)sin(3pi/4)
214 * = 1/sqrt(2) * (sin(x) - cos(x))
215 * sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4)
216 * = -1/sqrt(2) * (cos(x) + sin(x))
217 * To avoid cancellation, use
218 * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
219 * to compute the worse one.
220 */
221 /* END CSTYLED */
222 if (x > 1.0e2450L) { /* x+x may overflow */
223 ss = -s - c;
224 cc = s - c;
225 } else if (signbitl(s) != signbitl(c)) {
226 cc = s - c;
227 ss = cosl(x + x) / cc;
228 } else {
229 ss = -s - c;
230 cc = cosl(x + x) / ss;
231 }
232
233 /*
234 * j1(x) = 1/sqrt(pi*x) * (P(1,x)*cc - Q(1,x)*ss)
235 * y1(x) = 1/sqrt(pi*x) * (P(1,x)*ss + Q(1,x)*cc)
236 */
237 if (x > 1.0e91L)
238 return ((invsqrtpi * ss) / sqrtl(x));
239
240 return (invsqrtpi * (pone(x) * ss + qone(x) * cc) / sqrtl(x));
241 }
242
243 if (x <= tiny)
244 return (-tpi / x);
245
246 z = x * x;
247 u = u0[6];
248 v = v0[6] + z * v0[7];
249
250 for (i = 5; i >= 0; i--) {
251 u = u * z + u0[i];
252 v = v * z + v0[i];
253 }
254
255 return (x * (u / v) + tpi * (j1l(x) * logl(x) - one / x));
256 }
257
258 static const GENERIC pr0[12] = {
259 1.000000000000000000000000000000000000267e+0000L,
260 1.060717875045891455602180843276758003035e+0003L,
261 4.344347542892127024446687712181105852335e+0005L,
262 8.915680220724007016377924252717410457094e+0007L,
263 9.969502259938406062809873257569171272819e+0009L,
264 6.200290193138613035646510338707386316595e+0011L,
265 2.105978548788015119851815854422247330118e+0013L,
266 3.696635772784601239371730810311998368948e+0014L,
267 3.015913097920694682057958412534134515156e+0015L,
268 9.370298471339353098123277427328592725921e+0015L,
269 7.190349005196335967340799265074029443057e+0015L,
270 2.736097786240689996880391074927552517982e+0014L,
271 };
272
273 static const GENERIC ps0[11] = {
274 1.0e0L,
275 1.060600687545891455602180843276758095107e+0003L,
276 4.343106093416975589147153906505338900961e+0005L,
277 8.910605869002176566582072242244353399059e+0007L,
278 9.959122058635087888690713917622056540190e+0009L,
279 6.188744967234948231792482949171041843894e+0011L,
280 2.098863976953783506401759873801990304907e+0013L,
281 3.672870357018063196746729751479938908450e+0014L,
282 2.975538419246824921049011529574385888420e+0015L,
283 9.063657659995043205018686029284479837091e+0015L,
284 6.401953344314747916729366441508892711691e+0015L,
285 };
286
287 static const GENERIC pr1[12] = {
288 1.000000000000000000000023667524130660984e+0000L,
289 6.746154419979618754354803488126452971204e+0002L,
290 1.811210781083390154857018330296145970502e+0005L,
291 2.533098390379924268038005329095287842244e+0007L,
292 2.029683619805342145252338570875424600729e+0009L,
293 9.660859662192711465301069401598929980319e+0010L,
294 2.743396238644831519934098967716621316316e+0012L,
295 4.553097354140854377931023170263455246288e+0013L,
296 4.210245069852219757476169864974870720374e+0014L,
297 1.987334056229596485076645967176169801727e+0015L,
298 4.067120052787096893838970455751338930462e+0015L,
299 2.486539606380406398310845264910691398133e+0015L,
300 };
301
302 static const GENERIC ps1[14] = {
303 1.0e0L,
304 6.744982544979618754355808680196859521782e+0002L,
305 1.810421795396966762032155290441364740350e+0005L,
306 2.530986460644310651529583759699988435573e+0007L,
307 2.026743276048023121360249288818290224145e+0009L,
308 9.637461924407405935245269407052641341836e+0010L,
309 2.732378628423766417402292797028314160831e+0012L,
310 4.522345274960527124354844364012184278488e+0013L,
311 4.160650668341743132685335758415469856545e+0014L,
312 1.943730242988858208243492424892435901211e+0015L,
313 3.880228532692127989901131618598067450001e+0015L,
314 2.178020816161154615841000173683302999728e+0015L,
315 -8.994062666842225551554346698171600634173e+0013L,
316 1.368520368508851253495764806934619574990e+0013L,
317 };
318
319 static const GENERIC pr2[12] = {
320 1.000000000000000006938651621840396237282e+0000L,
321 3.658416291850404981407101077037948144698e+0002L,
322 5.267073772170356547709794670602812447537e+0004L,
323 3.912012101226837463014925210735894620442e+0006L,
324 1.651295648974103957193874928714180765625e+0008L,
325 4.114901144480797609972484998142146783499e+0009L,
326 6.092524309766036681542980572526335147672e+0010L,
327 5.263913178071282616719249969074134570577e+0011L,
328 2.538408581124324223367341020538081330994e+0012L,
329 6.288607929360291027895126983015365677648e+0012L,
330 6.848330048211148419047055075386525945280e+0012L,
331 2.290309646838867941423178163991423244690e+0012L,
332 };
333
334 static const GENERIC ps2[14] = {
335 1.0e0L,
336 3.657244416850405086459410165762319861856e+0002L,
337 5.262802358425023243992387075861237306312e+0004L,
338 3.905896813959919648136295861661483848364e+0006L,
339 1.646791907791461220742694842108202772763e+0008L,
340 4.096132803064256022224954120208201437344e+0009L,
341 6.046665195915950447544429445730680236759e+0010L,
342 5.198061739781991313414052212328653295168e+0011L,
343 2.484233851814333966401527626421254279796e+0012L,
344 6.047868806925315879339651539434315255940e+0012L,
345 6.333103831254091652501642567294101813354e+0012L,
346 1.875143098754284994467609936924685024968e+0012L,
347 -5.238330920563392692965412762508813601534e+0010L,
348 4.656888609439364725427789198383779259957e+0009L,
349 };
350
351 static const GENERIC pr3[13] = {
352 1.000000000000009336887318068056137842897e+0000L,
353 2.242719942728459588488051572002835729183e+0002L,
354 1.955450611382026550266257737331095691092e+0004L,
355 8.707143293993619899395400562409175590739e+0005L,
356 2.186267894487004565948324289010954505316e+0007L,
357 3.224328510541957792360691585667502864688e+0008L,
358 2.821057355151380597331792896882741364897e+0009L,
359 1.445371387295422404365584793796028979840e+0010L,
360 4.181743160669891357783011002656658107864e+0010L,
361 6.387371088767993119325536137794535513922e+0010L,
362 4.575619999412716078064070587767416436396e+0010L,
363 1.228415651211639160620284441690503550842e+0010L,
364 7.242170349875563053436050532153112882072e+0008L,
365 };
366
367 static const GENERIC ps3[13] = {
368 1.0e0L,
369 2.241548067728529551049804610486061401070e+0002L,
370 1.952838216795552145132137932931237181307e+0004L,
371 8.684574926493185744628127341069974575526e+0005L,
372 2.176357771067037962940853412819852189164e+0007L,
373 3.199958682356132977319258783167122100567e+0008L,
374 2.786218931525334687844675219914201872570e+0009L,
375 1.416283776951741549631417572317916039767e+0010L,
376 4.042962659271567948735676834609348842922e+0010L,
377 6.028168462646694510083847222968444402161e+0010L,
378 4.118410226794641413833887606580085281111e+0010L,
379 9.918735736297038430744161253338202230263e+0009L,
380 4.092967198238098023219124487437130332038e+0008L,
381 };
382
383 static const GENERIC pr4[13] = {
384 1.000000000001509220978157399042059553390e+0000L,
385 1.437551868378147851133499996323782607787e+0002L,
386 7.911335537418177296041518061404505428004e+0003L,
387 2.193710939115317214716518908935756104804e+0005L,
388 3.390662495136730962513489796538274984382e+0006L,
389 3.048655347929348891006070609293884274789e+0007L,
390 1.613781633489496606354045161527450975195e+0008L,
391 4.975089835037230277110156150038482159988e+0008L,
392 8.636047087015115403880904418339566323264e+0008L,
393 7.918202912328366140110671223076949101509e+0008L,
394 3.423294665798984733439650311722794853294e+0008L,
395 5.621904953441963961040503934782662613621e+0007L,
396 2.086303543310240260758670404509484499793e+0006L,
397 };
398
399 static const GENERIC ps4[13] = {
400 1.0e0L,
401 1.436379993384532371670493319591847362304e+0002L,
402 7.894647154785430678061053848847436659499e+0003L,
403 2.184659753392097529008981741550878586174e+0005L,
404 3.366109083305465176803513738147049499361e+0006L,
405 3.011911545968996817697665866587226343186e+0007L,
406 1.582262913779689851316760148459414895301e+0008L,
407 4.819268809494937919217938589530138201770e+0008L,
408 8.201355762990450679702837123432527154830e+0008L,
409 7.268232093982510937417446421282341425212e+0008L,
410 2.950911909015572933262131323934036480462e+0008L,
411 4.242839924305934423010858966540621219396e+0007L,
412 1.064387620445090779182117666330405186866e+0006L,
413 };
414
415 static const GENERIC pr5[13] = {
416 1.000000000102434805241171427253847353861e+0000L,
417 9.129332257083629259060502249025963234821e+0001L,
418 3.132238483586953037576119377504557191413e+0003L,
419 5.329782528269307971278943122454171107861e+0004L,
420 4.988460157184117790692873002103052944145e+0005L,
421 2.686602071615786816147010334256047469378e+0006L,
422 8.445418526028961197703799808701268301831e+0006L,
423 1.536575358646141157475725889907900827390e+0007L,
424 1.568405818236523821796862770586544811945e+0007L,
425 8.450876239888770102387618667362302173547e+0006L,
426 2.154414900139567328424026827163203446077e+0006L,
427 2.105656926565043898888460254808062352205e+0005L,
428 4.739165011023396507022134303736862812975e+0003L,
429 };
430
431 static const GENERIC ps5[13] = {
432 1.0e0L,
433 9.117613509595327476509152673394703847793e+0001L,
434 3.121697972484015639301279229281770795147e+0003L,
435 5.294447222735893568040911873834576440255e+0004L,
436 4.930368882192772335798256684110887882807e+0005L,
437 2.634854685641165298302167435798357437768e+0006L,
438 8.185462775400326393555896157031818280918e+0006L,
439 1.462417423080215192609668642663030667086e+0007L,
440 1.450624993985851675982860844153954896015e+0007L,
441 7.460467647561995283219086567162006113864e+0006L,
442 1.754210981405612478869227142579056338965e+0006L,
443 1.463286721155271971526264914524746699596e+0005L,
444 2.155894725796702015341211116579827039459e+0003L,
445 };
446
447 static const GENERIC pr6[13] = {
448 1.000000003564855546741735920315743157129e+0000L,
449 5.734003934862540458119423509909510288366e+0001L,
450 1.209572491935850486086559692291796887976e+0003L,
451 1.243398391422281247933674779163660286838e+0004L,
452 6.930996755181437937258220998601708278787e+0004L,
453 2.198067659532757598646722249966767620099e+0005L,
454 4.033659432712058633933179115820576858455e+0005L,
455 4.257759657219008027016047206574574358678e+0005L,
456 2.511917395876004349480721277445763916389e+0005L,
457 7.813756153070623654178731651381881953552e+0004L,
458 1.152069173381127881385588092905864352891e+0004L,
459 6.548580782804088553777816037551523398082e+0002L,
460 8.668725370116906132327542766127938496880e+0000L,
461 };
462
463 static const GENERIC ps6[13] = {
464 1.0e0L,
465 5.722285236357114566499221525736286205184e+0001L,
466 1.203010842878317935444582950620339570506e+0003L,
467 1.230058335378583550155825502172435371208e+0004L,
468 6.800998550607861288865300438648089894412e+0004L,
469 2.130767829599304262987769347536850885921e+0005L,
470 3.840483466643916681759936972992155310026e+0005L,
471 3.947432373459225542861819148108081160393e+0005L,
472 2.237816239393081111481588434457838526738e+0005L,
473 6.545820495124419723398946273790921540774e+0004L,
474 8.729563630320892741500726213278834737196e+0003L,
475 4.130762660291894753450174794196998813709e+0002L,
476 3.480368898672684645130335786015075595598e+0000L,
477 };
478
479 static const GENERIC sixteen = 16.0L;
480 static const GENERIC eight = 8.0L;
481 static const GENERIC huge = 1.0e30L;
482
483 static GENERIC
484 pone(GENERIC x)
485 {
486 GENERIC s, r, t, z;
487 int i;
488
489 if (x > huge)
490 return (one);
491
492 t = one / x;
493 z = t * t;
494
495 if (x > sixteen) {
496 r = z * pr0[11] + pr0[10];
497 s = ps0[10];
498
499 for (i = 9; i >= 0; i--) {
500 r = z * r + pr0[i];
501 s = z * s + ps0[i];
502 }
503 } else if (x > eight) {
504 r = pr1[11];
505 s = ps1[11] + z * (ps1[12] + z * ps1[13]);
506
507 for (i = 10; i >= 0; i--) {
508 r = z * r + pr1[i];
509 s = z * s + ps1[i];
510 }
511 } else if (x > five) {
512 r = pr2[11];
513 s = ps2[11] + z * (ps2[12] + z * ps2[13]);
514
515 for (i = 10; i >= 0; i--) {
516 r = z * r + pr2[i];
517 s = z * s + ps2[i];
518 }
519 } else if (x > 3.5L) {
520 r = pr3[12];
521 s = ps3[12];
522
523 for (i = 11; i >= 0; i--) {
524 r = z * r + pr3[i];
525 s = z * s + ps3[i];
526 }
527 } else if (x > 2.5L) {
528 r = pr4[12];
529 s = ps4[12];
530
531 for (i = 11; i >= 0; i--) {
532 r = z * r + pr4[i];
533 s = z * s + ps4[i];
534 }
535 } else if (x > (1.0L / 0.5625L)) {
536 r = pr5[12];
537 s = ps5[12];
538
539 for (i = 11; i >= 0; i--) {
540 r = z * r + pr5[i];
541 s = z * s + ps5[i];
542 }
543 } else { /* assume x > 1.28 */
544 r = pr6[12];
545 s = ps6[12];
546
547 for (i = 11; i >= 0; i--) {
548 r = z * r + pr6[i];
549 s = z * s + ps6[i];
550 }
551 }
552
553 return (r / s);
554 }
555 static const GENERIC qr0[12] = {
556 3.749999999999999999999999999999999971033e-0001L,
557 4.256726035237050601607682277433094262226e+0002L,
558 1.875976490812878489192409978945401066066e+0005L,
559 4.170314268048041914273603680317745592790e+0007L,
560 5.092750132543855817293451118974555746551e+0009L,
561 3.494749676278488654103505795794139483404e+0011L,
562 1.327062148257437316997667817096694173709e+0013L,
563 2.648993138273427226907503742066551150490e+0014L,
564 2.511695665909547412222430494473998127684e+0015L,
565 9.274694506662289043224310499164702306096e+0015L,
566 8.150904170663663829331320302911792892002e+0015L,
567 -5.001918733707662355772037829620388765122e+0014L,
568 };
569
570 static const GENERIC qs0[11] = {
571 1.0e0L,
572 1.135400380229880160428715273982155760093e+0003L,
573 5.005701183877126164326765545516590744360e+0005L,
574 1.113444200113712167984337603933040102987e+0008L,
575 1.361074819925223062778717565699039471124e+0010L,
576 9.355750985802849484438933905325982809653e+0011L,
577 3.563462786008988825003965543857998084828e+0013L,
578 7.155145113900094163648726863803802910454e+0014L,
579 6.871266835834472758055559013851843654113e+0015L,
580 2.622030899226736712644974988157345234092e+0016L,
581 2.602912729172876330650077021706139707746e+0016L,
582 };
583
584 static const GENERIC qr1[12] = {
585 3.749999999999999999997762458207284405806e-0001L,
586 2.697883998881706839929255517498189980485e+0002L,
587 7.755195925781028489386938870473834411019e+0004L,
588 1.166777762104017777198211072895528968355e+0007L,
589 1.011504772984321168320010084520261069362e+0009L,
590 5.246007703574156853577754571720205550010e+0010L,
591 1.637692549885592683166116551691266537647e+0012L,
592 3.022303623698185669912990310925039382495e+0013L,
593 3.154769927290655684846107030265909987946e+0014L,
594 1.715819913441554770089730934808123360921e+0015L,
595 4.165044355759732622273534445131736188510e+0015L,
596 3.151381420874174705643100381708086287596e+0015L,
597 };
598
599 static const GENERIC qs1[14] = {
600 1.0e0L,
601 7.197091705351218239785633172408276982828e+0002L,
602 2.070012799599548685544883041297609861055e+0005L,
603 3.117014815317656221871840152778458754516e+0007L,
604 2.705719678902554974863325877025902971727e+0009L,
605 1.406113614727345726925060648750867264098e+0011L,
606 4.403777536067131320363005978631674817359e+0012L,
607 8.170725690209322283061499386703167242894e+0013L,
608 8.609458844975495289227794126964431210566e+0014L,
609 4.766766367015473481257280600694952920204e+0015L,
610 1.202286587943342194863557940888115641650e+0016L,
611 1.012474328306200909525063936061756024120e+0016L,
612 6.183552022678917858273222879615824070703e+0014L,
613 -9.756731548558226997573737400988225722740e+0013L,
614 };
615
616 static const GENERIC qr2[12] = {
617 3.749999999999999481245647262226994293189e-0001L,
618 1.471366807289771354491181140167359026735e+0002L,
619 2.279432486768448220142080962843526951250e+0004L,
620 1.828943048523771225163804043356958285893e+0006L,
621 8.379828388647823135832220596417725010837e+0007L,
622 2.279814029335044024585393671278378022053e+0009L,
623 3.711653952257118120832817785271466441420e+0010L,
624 3.557650914518554549916730572553105048068e+0011L,
625 1.924583483146095896259774329498934160650e+0012L,
626 5.424386256063736390759567088291887140278e+0012L,
627 6.839325621241776786206509704671746841737e+0012L,
628 2.702169563144001166291686452305436313971e+0012L,
629 };
630
631 static const GENERIC qs2[14] = {
632 1.0e0L,
633 3.926379194439388135703211933895203191089e+0002L,
634 6.089148804106598297488336063007609312276e+0004L,
635 4.893546162973278583711376356041614150645e+0006L,
636 2.247571119114497845046388801813832219404e+0008L,
637 6.137635663350177751290469334200757872645e+0009L,
638 1.005115019784102856424493519524998953678e+0011L,
639 9.725664462014503832860151384604677240620e+0011L,
640 5.345525100485511116148634192844434636072e+0012L,
641 1.549944007398946691720862738173956994779e+0013L,
642 2.067148441178952625710302124163264760362e+0013L,
643 9.401565402641963611295119487242595462301e+0012L,
644 3.548217088622398274748837287769709374385e+0011L,
645 -2.934470341719047120076509938432417352365e+0010L,
646 };
647
648 static const GENERIC qr3[13] = {
649 3.749999999999412724084579833297451472091e-0001L,
650 9.058478580291706212422978492938435582527e+0001L,
651 8.524056033161038750461083666711724381171e+0003L,
652 4.105967158629109427753434569223631014730e+0005L,
653 1.118326603378531348259783091972623333657e+0007L,
654 1.794636683403578918528064904714132329343e+0008L,
655 1.714314157463635959556133236004368896724e+0009L,
656 9.622092032236084846572067257267661456030e+0009L,
657 3.057759524485859159957762858780768355020e+0010L,
658 5.129306780754798531609621454415938890020e+0010L,
659 3.999122002794961070680636194346316041352e+0010L,
660 1.122298454643493485989721564358100345388e+0010L,
661 5.603981987645989709668830968522362582221e+0008L,
662 };
663
664 static const GENERIC qs3[13] = {
665 1.0e0L,
666 2.418328663076578169836155170053634419922e+0002L,
667 2.279620205900121042587523541281272875520e+0004L,
668 1.100984222585729521470129014992217092794e+0006L,
669 3.010743223679247091004262516286654516282e+0007L,
670 4.860925542827367817289619265215599433996e+0008L,
671 4.686668111035348691982715864307839581243e+0009L,
672 2.668701788405102017427214705946730894074e+0010L,
673 8.677395746106802640390580944836650584903e+0010L,
674 1.511936455574951790658498795945106643036e+0011L,
675 1.260845604432623478002018696873608353093e+0011L,
676 4.052692278419853853911440231600864589805e+0010L,
677 2.965516519212226064983267822243329694729e+0009L,
678 };
679
680 static const GENERIC qr4[13] = {
681 3.749999999919234164154669754440123072618e-0001L,
682 5.844218580776819864791168253485055101858e+0001L,
683 3.489273514092912982675669411371435670220e+0003L,
684 1.050523637774575684509663430018995479594e+0005L,
685 1.764549172059701565500717319792780115289e+0006L,
686 1.725532438844133795028063102681497371154e+0007L,
687 9.938114847359778539965140247590176334874e+0007L,
688 3.331710808184595545396883770200772842314e+0008L,
689 6.271970557641881511609560444872797282698e+0008L,
690 6.188529798677357075020774923903737913285e+0008L,
691 2.821905302742849974509982167877885011629e+0008L,
692 4.615467358646911976773290256984329814896e+0007L,
693 1.348140608731546467396685802693380693275e+0006L,
694 };
695
696 static const GENERIC qs4[13] = {
697 1.0e0L,
698 1.561192663112345185261418296389902133372e+0002L,
699 9.346678031144098270547225423124213083072e+0003L,
700 2.825851246482293547838023847601704751590e+0005L,
701 4.776572711622156091710902891124911556293e+0006L,
702 4.715106953717135402977938048006267859302e+0007L,
703 2.753962350894311316439652227611209035193e+0008L,
704 9.428501434615463207768964787500411575223e+0008L,
705 1.832650858775206787088236896454141572617e+0009L,
706 1.901697378939743226948920874296595242257e+0009L,
707 9.433322226854293780627188599226380812725e+0008L,
708 1.808520540608671608680284520798858587370e+0008L,
709 7.983342331736662753157217446919462398008e+0006L,
710 };
711
712 static const GENERIC qr5[13] = {
713 3.749999995331364437028988850515190446719e-0001L,
714 3.739356381766559882677514593041627547911e+0001L,
715 1.399562500629413529355265462912819802551e+0003L,
716 2.594154053098947925345332218062210111753e+0004L,
717 2.640149879297408640394163979394594318371e+0005L,
718 1.542471854873199142031889093591449397995e+0006L,
719 5.242272868972053374067572098992335425895e+0006L,
720 1.025834487769410221329633071426044839935e+0007L,
721 1.116553924239448940142230579060124209622e+0007L,
722 6.318076065595910176374916303525884653514e+0006L,
723 1.641218086168640408527639735915512881785e+0006L,
724 1.522369793529178644168813882912134706444e+0005L,
725 2.526530541062297200914180060208669584055e+0003L,
726 };
727
728 static const GENERIC qs5[13] = {
729 1.0e0L,
730 9.998960735935075380397545659016287506660e+0001L,
731 3.758767417842043742686475060540416737562e+0003L,
732 7.013652806952306520121959742519780781653e+0004L,
733 7.208949808818615099246529616211730446850e+0005L,
734 4.272753927109614455417836186072202009252e+0006L,
735 1.482524411356470699336129814111025434703e+0007L,
736 2.988750366665678233425279237627700803473e+0007L,
737 3.396957890261080492694709150553619185065e+0007L,
738 2.050652487738593004111578091156304540386e+0007L,
739 5.900504120811732547616511555946279451316e+0006L,
740 6.563391409260160897024498082273183468347e+0005L,
741 1.692629845012790205348966731477187041419e+0004L,
742 };
743
744 static const GENERIC qr6[13] = {
745 3.749999861516664133157566870858975421296e-0001L,
746 2.367863756747764863120797431599473468918e+0001L,
747 5.476715802114976248882067325630793143777e+0002L,
748 6.143190357869842894025012945444096170251e+0003L,
749 3.716250534677997850513733595140463851730e+0004L,
750 1.270883463823876752138326905022875657430e+0005L,
751 2.495301449636814481646371665429083801388e+0005L,
752 2.789578988212952248340486296254398601942e+0005L,
753 1.718247946911109055931819087137397324634e+0005L,
754 5.458973214011665714330326732204106364229e+0004L,
755 7.912102686687948786048943339759596652813e+0003L,
756 4.077961006160866935722030715149087938091e+0002L,
757 3.765206972770245085551057237882528510428e+0000L,
758 };
759
760 static const GENERIC qs6[13] = {
761 1.0e0L, 6.341646532940517305641893852673926809601e+0001L,
762 1.477058277414040790932597537920671025359e+0003L,
763 1.674406564031044491436044253393536487604e+0004L,
764 1.028516501369755949895050806908994650768e+0005L,
765 3.593620042532885295087463507733285434207e+0005L,
766 7.267924991381020915185873399453724799625e+0005L,
767 8.462277510768818399961191426205006083088e+0005L,
768 5.514399892230892163373611895645500250514e+0005L,
769 1.898084241009259353540620272932188102299e+0005L,
770 3.102941242117739015721984123081026253068e+0004L,
771 1.958971184431466907681440650181421086143e+0003L,
772 2.878853357310495087181721613889455121867e+0001L,
773 };
774
775 static GENERIC
776 qone(GENERIC x)
777 {
778 GENERIC s, r, t, z;
779 int i;
780
781 if (x > huge)
782 return (0.375L / x);
783
784 t = one / x;
785 z = t * t;
786
787 if (x > sixteen) {
788 r = z * qr0[11] + qr0[10];
789 s = qs0[10];
790
791 for (i = 9; i >= 0; i--) {
792 r = z * r + qr0[i];
793 s = z * s + qs0[i];
794 }
795 } else if (x > eight) {
796 r = qr1[11];
797 s = qs1[11] + z * (qs1[12] + z * qs1[13]);
798
799 for (i = 10; i >= 0; i--) {
800 r = z * r + qr1[i];
801 s = z * s + qs1[i];
802 }
803 } else if (x > five) { /* x > 5.0 */
804 r = qr2[11];
805 s = qs2[11] + z * (qs2[12] + z * qs2[13]);
806
807 for (i = 10; i >= 0; i--) {
808 r = z * r + qr2[i];
809 s = z * s + qs2[i];
810 }
811 } else if (x > 3.5L) {
812 r = qr3[12];
813 s = qs3[12];
814
815 for (i = 11; i >= 0; i--) {
816 r = z * r + qr3[i];
817 s = z * s + qs3[i];
818 }
819 } else if (x > 2.5L) {
820 r = qr4[12];
821 s = qs4[12];
822
823 for (i = 11; i >= 0; i--) {
824 r = z * r + qr4[i];
825 s = z * s + qs4[i];
826 }
827 } else if (x > (1.0L / 0.5625L)) {
828 r = qr5[12];
829 s = qs5[12];
830
831 for (i = 11; i >= 0; i--) {
832 r = z * r + qr5[i];
833 s = z * s + qs5[i];
834 }
835 } else { /* assume x > 1.28 */
836 r = qr6[12];
837 s = qs6[12];
838
839 for (i = 11; i >= 0; i--) {
840 r = z * r + qr6[i];
841 s = z * s + qs6[i];
842 }
843 }
844
845 return (t * (r / s));
846 }
|