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