OpenCores
URL https://opencores.org/ocsvn/openrisc/openrisc/trunk

Subversion Repositories openrisc

[/] [openrisc/] [trunk/] [gnu-dev/] [or1k-gcc/] [libquadmath/] [math/] [j1q.c] - Blame information for rev 740

Details | Compare with Previous | View Log

Line No. Rev Author Line
1 740 jeremybenn
/*                                                      j1l.c
2
 *
3
 *      Bessel function of order one
4
 *
5
 *
6
 *
7
 * SYNOPSIS:
8
 *
9
 * long double x, y, j1l();
10
 *
11
 * y = j1l( x );
12
 *
13
 *
14
 *
15
 * DESCRIPTION:
16
 *
17
 * Returns Bessel function of first kind, order one of the argument.
18
 *
19
 * The domain is divided into two major intervals [0, 2] and
20
 * (2, infinity). In the first interval the rational approximation is
21
 * J1(x) = .5x + x x^2 R(x^2)
22
 *
23
 * The second interval is further partitioned into eight equal segments
24
 * of 1/x.
25
 * J1(x) = sqrt(2/(pi x)) (P1(x) cos(X) - Q1(x) sin(X)),
26
 * X = x - 3 pi / 4,
27
 *
28
 * and the auxiliary functions are given by
29
 *
30
 * J1(x)cos(X) + Y1(x)sin(X) = sqrt( 2/(pi x)) P1(x),
31
 * P1(x) = 1 + 1/x^2 R(1/x^2)
32
 *
33
 * Y1(x)cos(X) - J1(x)sin(X) = sqrt( 2/(pi x)) Q1(x),
34
 * Q1(x) = 1/x (.375 + 1/x^2 R(1/x^2)).
35
 *
36
 *
37
 *
38
 * ACCURACY:
39
 *
40
 *                      Absolute error:
41
 * arithmetic   domain      # trials      peak         rms
42
 *    IEEE      0, 30       100000      2.8e-34      2.7e-35
43
 *
44
 *
45
 */
46
 
47
/*                                                      y1l.c
48
 *
49
 *      Bessel function of the second kind, order one
50
 *
51
 *
52
 *
53
 * SYNOPSIS:
54
 *
55
 * double x, y, y1l();
56
 *
57
 * y = y1l( x );
58
 *
59
 *
60
 *
61
 * DESCRIPTION:
62
 *
63
 * Returns Bessel function of the second kind, of order
64
 * one, of the argument.
65
 *
66
 * The domain is divided into two major intervals [0, 2] and
67
 * (2, infinity). In the first interval the rational approximation is
68
 * Y1(x) = 2/pi * (log(x) * J1(x) - 1/x) + x R(x^2) .
69
 * In the second interval the approximation is the same as for J1(x), and
70
 * Y1(x) = sqrt(2/(pi x)) (P1(x) sin(X) + Q1(x) cos(X)),
71
 * X = x - 3 pi / 4.
72
 *
73
 * ACCURACY:
74
 *
75
 *  Absolute error, when y0(x) < 1; else relative error:
76
 *
77
 * arithmetic   domain     # trials      peak         rms
78
 *    IEEE      0, 30       100000      2.7e-34     2.9e-35
79
 *
80
 */
81
 
82
/* Copyright 2001 by Stephen L. Moshier (moshier@na-net.onrl.gov).
83
 
84
    This library is free software; you can redistribute it and/or
85
    modify it under the terms of the GNU Lesser General Public
86
    License as published by the Free Software Foundation; either
87
    version 2.1 of the License, or (at your option) any later version.
88
 
89
    This library is distributed in the hope that it will be useful,
90
    but WITHOUT ANY WARRANTY; without even the implied warranty of
91
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
92
    Lesser General Public License for more details.
93
 
94
    You should have received a copy of the GNU Lesser General Public
95
    License along with this library; if not, write to the Free Software
96
    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307  USA */
97
 
98
#include "quadmath-imp.h"
99
 
100
/* 1 / sqrt(pi) */
101
static const __float128 ONEOSQPI = 5.6418958354775628694807945156077258584405E-1Q;
102
/* 2 / pi */
103
static const __float128 TWOOPI = 6.3661977236758134307553505349005744813784E-1Q;
104
static const __float128 zero = 0.0Q;
105
 
106
/* J1(x) = .5x + x x^2 R(x^2)
107
   Peak relative error 1.9e-35
108
 
109
#define NJ0_2N 6
110
static const __float128 J0_2N[NJ0_2N + 1] = {
111
 -5.943799577386942855938508697619735179660E16Q,
112
  1.812087021305009192259946997014044074711E15Q,
113
 -2.761698314264509665075127515729146460895E13Q,
114
  2.091089497823600978949389109350658815972E11Q,
115
 -8.546413231387036372945453565654130054307E8Q,
116
  1.797229225249742247475464052741320612261E6Q,
117
 -1.559552840946694171346552770008812083969E3Q
118
};
119
#define NJ0_2D 6
120
static const __float128 J0_2D[NJ0_2D + 1] = {
121
  9.510079323819108569501613916191477479397E17Q,
122
  1.063193817503280529676423936545854693915E16Q,
123
  5.934143516050192600795972192791775226920E13Q,
124
  2.168000911950620999091479265214368352883E11Q,
125
  5.673775894803172808323058205986256928794E8Q,
126
  1.080329960080981204840966206372671147224E6Q,
127
  1.411951256636576283942477881535283304912E3Q,
128
 /* 1.000000000000000000000000000000000000000E0Q */
129
};
130
 
131
/* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2),
132
 
133
   Peak relative error 3.6e-36  */
134
#define NP16_IN 9
135
static const __float128 P16_IN[NP16_IN + 1] = {
136
  5.143674369359646114999545149085139822905E-16Q,
137
  4.836645664124562546056389268546233577376E-13Q,
138
  1.730945562285804805325011561498453013673E-10Q,
139
  3.047976856147077889834905908605310585810E-8Q,
140
  2.855227609107969710407464739188141162386E-6Q,
141
  1.439362407936705484122143713643023998457E-4Q,
142
  3.774489768532936551500999699815873422073E-3Q,
143
  4.723962172984642566142399678920790598426E-2Q,
144
  2.359289678988743939925017240478818248735E-1Q,
145
  3.032580002220628812728954785118117124520E-1Q,
146
};
147
#define NP16_ID 9
148
static const __float128 P16_ID[NP16_ID + 1] = {
149
  4.389268795186898018132945193912677177553E-15Q,
150
  4.132671824807454334388868363256830961655E-12Q,
151
  1.482133328179508835835963635130894413136E-9Q,
152
  2.618941412861122118906353737117067376236E-7Q,
153
  2.467854246740858470815714426201888034270E-5Q,
154
  1.257192927368839847825938545925340230490E-3Q,
155
  3.362739031941574274949719324644120720341E-2Q,
156
  4.384458231338934105875343439265370178858E-1Q,
157
  2.412830809841095249170909628197264854651E0Q,
158
  4.176078204111348059102962617368214856874E0Q,
159
 /* 1.000000000000000000000000000000000000000E0 */
160
};
161
 
162
/* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2),
163
    0.0625 <= 1/x <= 0.125
164
    Peak relative error 1.9e-36  */
165
#define NP8_16N 11
166
static const __float128 P8_16N[NP8_16N + 1] = {
167
  2.984612480763362345647303274082071598135E-16Q,
168
  1.923651877544126103941232173085475682334E-13Q,
169
  4.881258879388869396043760693256024307743E-11Q,
170
  6.368866572475045408480898921866869811889E-9Q,
171
  4.684818344104910450523906967821090796737E-7Q,
172
  2.005177298271593587095982211091300382796E-5Q,
173
  4.979808067163957634120681477207147536182E-4Q,
174
  6.946005761642579085284689047091173581127E-3Q,
175
  5.074601112955765012750207555985299026204E-2Q,
176
  1.698599455896180893191766195194231825379E-1Q,
177
  1.957536905259237627737222775573623779638E-1Q,
178
  2.991314703282528370270179989044994319374E-2Q,
179
};
180
#define NP8_16D 10
181
static const __float128 P8_16D[NP8_16D + 1] = {
182
  2.546869316918069202079580939942463010937E-15Q,
183
  1.644650111942455804019788382157745229955E-12Q,
184
  4.185430770291694079925607420808011147173E-10Q,
185
  5.485331966975218025368698195861074143153E-8Q,
186
  4.062884421686912042335466327098932678905E-6Q,
187
  1.758139661060905948870523641319556816772E-4Q,
188
  4.445143889306356207566032244985607493096E-3Q,
189
  6.391901016293512632765621532571159071158E-2Q,
190
  4.933040207519900471177016015718145795434E-1Q,
191
  1.839144086168947712971630337250761842976E0Q,
192
  2.715120873995490920415616716916149586579E0Q,
193
 /* 1.000000000000000000000000000000000000000E0 */
194
};
195
 
196
/* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2),
197
  0.125 <= 1/x <= 0.1875
198
  Peak relative error 1.3e-36  */
199
#define NP5_8N 10
200
static const __float128 P5_8N[NP5_8N + 1] = {
201
  2.837678373978003452653763806968237227234E-12Q,
202
  9.726641165590364928442128579282742354806E-10Q,
203
  1.284408003604131382028112171490633956539E-7Q,
204
  8.524624695868291291250573339272194285008E-6Q,
205
  3.111516908953172249853673787748841282846E-4Q,
206
  6.423175156126364104172801983096596409176E-3Q,
207
  7.430220589989104581004416356260692450652E-2Q,
208
  4.608315409833682489016656279567605536619E-1Q,
209
  1.396870223510964882676225042258855977512E0Q,
210
  1.718500293904122365894630460672081526236E0Q,
211
  5.465927698800862172307352821870223855365E-1Q
212
};
213
#define NP5_8D 10
214
static const __float128 P5_8D[NP5_8D + 1] = {
215
  2.421485545794616609951168511612060482715E-11Q,
216
  8.329862750896452929030058039752327232310E-9Q,
217
  1.106137992233383429630592081375289010720E-6Q,
218
  7.405786153760681090127497796448503306939E-5Q,
219
  2.740364785433195322492093333127633465227E-3Q,
220
  5.781246470403095224872243564165254652198E-2Q,
221
  6.927711353039742469918754111511109983546E-1Q,
222
  4.558679283460430281188304515922826156690E0Q,
223
  1.534468499844879487013168065728837900009E1Q,
224
  2.313927430889218597919624843161569422745E1Q,
225
  1.194506341319498844336768473218382828637E1Q,
226
 /* 1.000000000000000000000000000000000000000E0 */
227
};
228
 
229
/* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2),
230
   Peak relative error 1.4e-36
231
   0.1875 <= 1/x <= 0.25  */
232
#define NP4_5N 10
233
static const __float128 P4_5N[NP4_5N + 1] = {
234
  1.846029078268368685834261260420933914621E-10Q,
235
  3.916295939611376119377869680335444207768E-8Q,
236
  3.122158792018920627984597530935323997312E-6Q,
237
  1.218073444893078303994045653603392272450E-4Q,
238
  2.536420827983485448140477159977981844883E-3Q,
239
  2.883011322006690823959367922241169171315E-2Q,
240
  1.755255190734902907438042414495469810830E-1Q,
241
  5.379317079922628599870898285488723736599E-1Q,
242
  7.284904050194300773890303361501726561938E-1Q,
243
  3.270110346613085348094396323925000362813E-1Q,
244
  1.804473805689725610052078464951722064757E-2Q,
245
};
246
#define NP4_5D 9
247
static const __float128 P4_5D[NP4_5D + 1] = {
248
  1.575278146806816970152174364308980863569E-9Q,
249
  3.361289173657099516191331123405675054321E-7Q,
250
  2.704692281550877810424745289838790693708E-5Q,
251
  1.070854930483999749316546199273521063543E-3Q,
252
  2.282373093495295842598097265627962125411E-2Q,
253
  2.692025460665354148328762368240343249830E-1Q,
254
  1.739892942593664447220951225734811133759E0Q,
255
  5.890727576752230385342377570386657229324E0Q,
256
  9.517442287057841500750256954117735128153E0Q,
257
  6.100616353935338240775363403030137736013E0Q,
258
 /* 1.000000000000000000000000000000000000000E0 */
259
};
260
 
261
/* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2),
262
   Peak relative error 3.0e-36
263
   0.25 <= 1/x <= 0.3125  */
264
#define NP3r2_4N 9
265
static const __float128 P3r2_4N[NP3r2_4N + 1] = {
266
  8.240803130988044478595580300846665863782E-8Q,
267
  1.179418958381961224222969866406483744580E-5Q,
268
  6.179787320956386624336959112503824397755E-4Q,
269
  1.540270833608687596420595830747166658383E-2Q,
270
  1.983904219491512618376375619598837355076E-1Q,
271
  1.341465722692038870390470651608301155565E0Q,
272
  4.617865326696612898792238245990854646057E0Q,
273
  7.435574801812346424460233180412308000587E0Q,
274
  4.671327027414635292514599201278557680420E0Q,
275
  7.299530852495776936690976966995187714739E-1Q,
276
};
277
#define NP3r2_4D 9
278
static const __float128 P3r2_4D[NP3r2_4D + 1] = {
279
  7.032152009675729604487575753279187576521E-7Q,
280
  1.015090352324577615777511269928856742848E-4Q,
281
  5.394262184808448484302067955186308730620E-3Q,
282
  1.375291438480256110455809354836988584325E-1Q,
283
  1.836247144461106304788160919310404376670E0Q,
284
  1.314378564254376655001094503090935880349E1Q,
285
  4.957184590465712006934452500894672343488E1Q,
286
  9.287394244300647738855415178790263465398E1Q,
287
  7.652563275535900609085229286020552768399E1Q,
288
  2.147042473003074533150718117770093209096E1Q,
289
 /* 1.000000000000000000000000000000000000000E0 */
290
};
291
 
292
/* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2),
293
   Peak relative error 1.0e-35
294
   0.3125 <= 1/x <= 0.375  */
295
#define NP2r7_3r2N 9
296
static const __float128 P2r7_3r2N[NP2r7_3r2N + 1] = {
297
  4.599033469240421554219816935160627085991E-7Q,
298
  4.665724440345003914596647144630893997284E-5Q,
299
  1.684348845667764271596142716944374892756E-3Q,
300
  2.802446446884455707845985913454440176223E-2Q,
301
  2.321937586453963310008279956042545173930E-1Q,
302
  9.640277413988055668692438709376437553804E-1Q,
303
  1.911021064710270904508663334033003246028E0Q,
304
  1.600811610164341450262992138893970224971E0Q,
305
  4.266299218652587901171386591543457861138E-1Q,
306
  1.316470424456061252962568223251247207325E-2Q,
307
};
308
#define NP2r7_3r2D 8
309
static const __float128 P2r7_3r2D[NP2r7_3r2D + 1] = {
310
  3.924508608545520758883457108453520099610E-6Q,
311
  4.029707889408829273226495756222078039823E-4Q,
312
  1.484629715787703260797886463307469600219E-2Q,
313
  2.553136379967180865331706538897231588685E-1Q,
314
  2.229457223891676394409880026887106228740E0Q,
315
  1.005708903856384091956550845198392117318E1Q,
316
  2.277082659664386953166629360352385889558E1Q,
317
  2.384726835193630788249826630376533988245E1Q,
318
  9.700989749041320895890113781610939632410E0Q,
319
 /* 1.000000000000000000000000000000000000000E0 */
320
};
321
 
322
/* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2),
323
   Peak relative error 1.7e-36
324
   0.3125 <= 1/x <= 0.4375  */
325
#define NP2r3_2r7N 9
326
static const __float128 P2r3_2r7N[NP2r3_2r7N + 1] = {
327
  3.916766777108274628543759603786857387402E-6Q,
328
  3.212176636756546217390661984304645137013E-4Q,
329
  9.255768488524816445220126081207248947118E-3Q,
330
  1.214853146369078277453080641911700735354E-1Q,
331
  7.855163309847214136198449861311404633665E-1Q,
332
  2.520058073282978403655488662066019816540E0Q,
333
  3.825136484837545257209234285382183711466E0Q,
334
  2.432569427554248006229715163865569506873E0Q,
335
  4.877934835018231178495030117729800489743E-1Q,
336
  1.109902737860249670981355149101343427885E-2Q,
337
};
338
#define NP2r3_2r7D 8
339
static const __float128 P2r3_2r7D[NP2r3_2r7D + 1] = {
340
  3.342307880794065640312646341190547184461E-5Q,
341
  2.782182891138893201544978009012096558265E-3Q,
342
  8.221304931614200702142049236141249929207E-2Q,
343
  1.123728246291165812392918571987858010949E0Q,
344
  7.740482453652715577233858317133423434590E0Q,
345
  2.737624677567945952953322566311201919139E1Q,
346
  4.837181477096062403118304137851260715475E1Q,
347
  3.941098643468580791437772701093795299274E1Q,
348
  1.245821247166544627558323920382547533630E1Q,
349
 /* 1.000000000000000000000000000000000000000E0 */
350
};
351
 
352
/* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2),
353
   Peak relative error 1.7e-35
354
   0.4375 <= 1/x <= 0.5  */
355
#define NP2_2r3N 8
356
static const __float128 P2_2r3N[NP2_2r3N + 1] = {
357
  3.397930802851248553545191160608731940751E-4Q,
358
  2.104020902735482418784312825637833698217E-2Q,
359
  4.442291771608095963935342749477836181939E-1Q,
360
  4.131797328716583282869183304291833754967E0Q,
361
  1.819920169779026500146134832455189917589E1Q,
362
  3.781779616522937565300309684282401791291E1Q,
363
  3.459605449728864218972931220783543410347E1Q,
364
  1.173594248397603882049066603238568316561E1Q,
365
  9.455702270242780642835086549285560316461E-1Q,
366
};
367
#define NP2_2r3D 8
368
static const __float128 P2_2r3D[NP2_2r3D + 1] = {
369
  2.899568897241432883079888249845707400614E-3Q,
370
  1.831107138190848460767699919531132426356E-1Q,
371
  3.999350044057883839080258832758908825165E0Q,
372
  3.929041535867957938340569419874195303712E1Q,
373
  1.884245613422523323068802689915538908291E2Q,
374
  4.461469948819229734353852978424629815929E2Q,
375
  5.004998753999796821224085972610636347903E2Q,
376
  2.386342520092608513170837883757163414100E2Q,
377
  3.791322528149347975999851588922424189957E1Q,
378
 /* 1.000000000000000000000000000000000000000E0 */
379
};
380
 
381
/* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
382
   Q1(x) = 1/x (.375 + 1/x^2 R(1/x^2)),
383
   Peak relative error 8.0e-36
384
 
385
#define NQ16_IN 10
386
static const __float128 Q16_IN[NQ16_IN + 1] = {
387
  -3.917420835712508001321875734030357393421E-18Q,
388
  -4.440311387483014485304387406538069930457E-15Q,
389
  -1.951635424076926487780929645954007139616E-12Q,
390
  -4.318256438421012555040546775651612810513E-10Q,
391
  -5.231244131926180765270446557146989238020E-8Q,
392
  -3.540072702902043752460711989234732357653E-6Q,
393
  -1.311017536555269966928228052917534882984E-4Q,
394
  -2.495184669674631806622008769674827575088E-3Q,
395
  -2.141868222987209028118086708697998506716E-2Q,
396
  -6.184031415202148901863605871197272650090E-2Q,
397
  -1.922298704033332356899546792898156493887E-2Q,
398
};
399
#define NQ16_ID 9
400
static const __float128 Q16_ID[NQ16_ID + 1] = {
401
  3.820418034066293517479619763498400162314E-17Q,
402
  4.340702810799239909648911373329149354911E-14Q,
403
  1.914985356383416140706179933075303538524E-11Q,
404
  4.262333682610888819476498617261895474330E-9Q,
405
  5.213481314722233980346462747902942182792E-7Q,
406
  3.585741697694069399299005316809954590558E-5Q,
407
  1.366513429642842006385029778105539457546E-3Q,
408
  2.745282599850704662726337474371355160594E-2Q,
409
  2.637644521611867647651200098449903330074E-1Q,
410
  1.006953426110765984590782655598680488746E0Q,
411
 /* 1.000000000000000000000000000000000000000E0 */
412
 };
413
 
414
/* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
415
   Q1(x) = 1/x (.375 + 1/x^2 R(1/x^2)),
416
   Peak relative error 1.9e-36
417
   0.0625 <= 1/x <= 0.125  */
418
#define NQ8_16N 11
419
static const __float128 Q8_16N[NQ8_16N + 1] = {
420
  -2.028630366670228670781362543615221542291E-17Q,
421
  -1.519634620380959966438130374006858864624E-14Q,
422
  -4.540596528116104986388796594639405114524E-12Q,
423
  -7.085151756671466559280490913558388648274E-10Q,
424
  -6.351062671323970823761883833531546885452E-8Q,
425
  -3.390817171111032905297982523519503522491E-6Q,
426
  -1.082340897018886970282138836861233213972E-4Q,
427
  -2.020120801187226444822977006648252379508E-3Q,
428
  -2.093169910981725694937457070649605557555E-2Q,
429
  -1.092176538874275712359269481414448063393E-1Q,
430
  -2.374790947854765809203590474789108718733E-1Q,
431
  -1.365364204556573800719985118029601401323E-1Q,
432
};
433
#define NQ8_16D 11
434
static const __float128 Q8_16D[NQ8_16D + 1] = {
435
  1.978397614733632533581207058069628242280E-16Q,
436
  1.487361156806202736877009608336766720560E-13Q,
437
  4.468041406888412086042576067133365913456E-11Q,
438
  7.027822074821007443672290507210594648877E-9Q,
439
  6.375740580686101224127290062867976007374E-7Q,
440
  3.466887658320002225888644977076410421940E-5Q,
441
  1.138625640905289601186353909213719596986E-3Q,
442
  2.224470799470414663443449818235008486439E-2Q,
443
  2.487052928527244907490589787691478482358E-1Q,
444
  1.483927406564349124649083853892380899217E0Q,
445
  4.182773513276056975777258788903489507705E0Q,
446
  4.419665392573449746043880892524360870944E0Q,
447
 /* 1.000000000000000000000000000000000000000E0 */
448
};
449
 
450
/* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
451
   Q1(x) = 1/x (.375 + 1/x^2 R(1/x^2)),
452
   Peak relative error 1.5e-35
453
   0.125 <= 1/x <= 0.1875  */
454
#define NQ5_8N 10
455
static const __float128 Q5_8N[NQ5_8N + 1] = {
456
  -3.656082407740970534915918390488336879763E-13Q,
457
  -1.344660308497244804752334556734121771023E-10Q,
458
  -1.909765035234071738548629788698150760791E-8Q,
459
  -1.366668038160120210269389551283666716453E-6Q,
460
  -5.392327355984269366895210704976314135683E-5Q,
461
  -1.206268245713024564674432357634540343884E-3Q,
462
  -1.515456784370354374066417703736088291287E-2Q,
463
  -1.022454301137286306933217746545237098518E-1Q,
464
  -3.373438906472495080504907858424251082240E-1Q,
465
  -4.510782522110845697262323973549178453405E-1Q,
466
  -1.549000892545288676809660828213589804884E-1Q,
467
};
468
#define NQ5_8D 10
469
static const __float128 Q5_8D[NQ5_8D + 1] = {
470
  3.565550843359501079050699598913828460036E-12Q,
471
  1.321016015556560621591847454285330528045E-9Q,
472
  1.897542728662346479999969679234270605975E-7Q,
473
  1.381720283068706710298734234287456219474E-5Q,
474
  5.599248147286524662305325795203422873725E-4Q,
475
  1.305442352653121436697064782499122164843E-2Q,
476
  1.750234079626943298160445750078631894985E-1Q,
477
  1.311420542073436520965439883806946678491E0Q,
478
  5.162757689856842406744504211089724926650E0Q,
479
  9.527760296384704425618556332087850581308E0Q,
480
  6.604648207463236667912921642545100248584E0Q,
481
 /* 1.000000000000000000000000000000000000000E0 */
482
};
483
 
484
/* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
485
   Q1(x) = 1/x (.375 + 1/x^2 R(1/x^2)),
486
   Peak relative error 1.3e-35
487
   0.1875 <= 1/x <= 0.25  */
488
#define NQ4_5N 10
489
static const __float128 Q4_5N[NQ4_5N + 1] = {
490
  -4.079513568708891749424783046520200903755E-11Q,
491
  -9.326548104106791766891812583019664893311E-9Q,
492
  -8.016795121318423066292906123815687003356E-7Q,
493
  -3.372350544043594415609295225664186750995E-5Q,
494
  -7.566238665947967882207277686375417983917E-4Q,
495
  -9.248861580055565402130441618521591282617E-3Q,
496
  -6.033106131055851432267702948850231270338E-2Q,
497
  -1.966908754799996793730369265431584303447E-1Q,
498
  -2.791062741179964150755788226623462207560E-1Q,
499
  -1.255478605849190549914610121863534191666E-1Q,
500
  -4.320429862021265463213168186061696944062E-3Q,
501
};
502
#define NQ4_5D 9
503
static const __float128 Q4_5D[NQ4_5D + 1] = {
504
  3.978497042580921479003851216297330701056E-10Q,
505
  9.203304163828145809278568906420772246666E-8Q,
506
  8.059685467088175644915010485174545743798E-6Q,
507
  3.490187375993956409171098277561669167446E-4Q,
508
  8.189109654456872150100501732073810028829E-3Q,
509
  1.072572867311023640958725265762483033769E-1Q,
510
  7.790606862409960053675717185714576937994E-1Q,
511
  3.016049768232011196434185423512777656328E0Q,
512
  5.722963851442769787733717162314477949360E0Q,
513
  4.510527838428473279647251350931380867663E0Q,
514
 /* 1.000000000000000000000000000000000000000E0 */
515
};
516
 
517
/* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
518
   Q1(x) = 1/x (.375 + 1/x^2 R(1/x^2)),
519
   Peak relative error 2.1e-35
520
   0.25 <= 1/x <= 0.3125  */
521
#define NQ3r2_4N 9
522
static const __float128 Q3r2_4N[NQ3r2_4N + 1] = {
523
  -1.087480809271383885936921889040388133627E-8Q,
524
  -1.690067828697463740906962973479310170932E-6Q,
525
  -9.608064416995105532790745641974762550982E-5Q,
526
  -2.594198839156517191858208513873961837410E-3Q,
527
  -3.610954144421543968160459863048062977822E-2Q,
528
  -2.629866798251843212210482269563961685666E-1Q,
529
  -9.709186825881775885917984975685752956660E-1Q,
530
  -1.667521829918185121727268867619982417317E0Q,
531
  -1.109255082925540057138766105229900943501E0Q,
532
  -1.812932453006641348145049323713469043328E-1Q,
533
};
534
#define NQ3r2_4D 9
535
static const __float128 Q3r2_4D[NQ3r2_4D + 1] = {
536
  1.060552717496912381388763753841473407026E-7Q,
537
  1.676928002024920520786883649102388708024E-5Q,
538
  9.803481712245420839301400601140812255737E-4Q,
539
  2.765559874262309494758505158089249012930E-2Q,
540
  4.117921827792571791298862613287549140706E-1Q,
541
  3.323769515244751267093378361930279161413E0Q,
542
  1.436602494405814164724810151689705353670E1Q,
543
  3.163087869617098638064881410646782408297E1Q,
544
  3.198181264977021649489103980298349589419E1Q,
545
  1.203649258862068431199471076202897823272E1Q,
546
 /* 1.000000000000000000000000000000000000000E0  */
547
};
548
 
549
/* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
550
   Q1(x) = 1/x (.375 + 1/x^2 R(1/x^2)),
551
   Peak relative error 1.6e-36
552
   0.3125 <= 1/x <= 0.375  */
553
#define NQ2r7_3r2N 9
554
static const __float128 Q2r7_3r2N[NQ2r7_3r2N + 1] = {
555
  -1.723405393982209853244278760171643219530E-7Q,
556
  -2.090508758514655456365709712333460087442E-5Q,
557
  -9.140104013370974823232873472192719263019E-4Q,
558
  -1.871349499990714843332742160292474780128E-2Q,
559
  -1.948930738119938669637865956162512983416E-1Q,
560
  -1.048764684978978127908439526343174139788E0Q,
561
  -2.827714929925679500237476105843643064698E0Q,
562
  -3.508761569156476114276988181329773987314E0Q,
563
  -1.669332202790211090973255098624488308989E0Q,
564
  -1.930796319299022954013840684651016077770E-1Q,
565
};
566
#define NQ2r7_3r2D 9
567
static const __float128 Q2r7_3r2D[NQ2r7_3r2D + 1] = {
568
  1.680730662300831976234547482334347983474E-6Q,
569
  2.084241442440551016475972218719621841120E-4Q,
570
  9.445316642108367479043541702688736295579E-3Q,
571
  2.044637889456631896650179477133252184672E-1Q,
572
  2.316091982244297350829522534435350078205E0Q,
573
  1.412031891783015085196708811890448488865E1Q,
574
  4.583830154673223384837091077279595496149E1Q,
575
  7.549520609270909439885998474045974122261E1Q,
576
  5.697605832808113367197494052388203310638E1Q,
577
  1.601496240876192444526383314589371686234E1Q,
578
  /* 1.000000000000000000000000000000000000000E0 */
579
};
580
 
581
/* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
582
   Q1(x) = 1/x (.375 + 1/x^2 R(1/x^2)),
583
   Peak relative error 9.5e-36
584
   0.375 <= 1/x <= 0.4375  */
585
#define NQ2r3_2r7N 9
586
static const __float128 Q2r3_2r7N[NQ2r3_2r7N + 1] = {
587
  -8.603042076329122085722385914954878953775E-7Q,
588
  -7.701746260451647874214968882605186675720E-5Q,
589
  -2.407932004380727587382493696877569654271E-3Q,
590
  -3.403434217607634279028110636919987224188E-2Q,
591
  -2.348707332185238159192422084985713102877E-1Q,
592
  -7.957498841538254916147095255700637463207E-1Q,
593
  -1.258469078442635106431098063707934348577E0Q,
594
  -8.162415474676345812459353639449971369890E-1Q,
595
  -1.581783890269379690141513949609572806898E-1Q,
596
  -1.890595651683552228232308756569450822905E-3Q,
597
};
598
#define NQ2r3_2r7D 8
599
static const __float128 Q2r3_2r7D[NQ2r3_2r7D + 1] = {
600
  8.390017524798316921170710533381568175665E-6Q,
601
  7.738148683730826286477254659973968763659E-4Q,
602
  2.541480810958665794368759558791634341779E-2Q,
603
  3.878879789711276799058486068562386244873E-1Q,
604
  3.003783779325811292142957336802456109333E0Q,
605
  1.206480374773322029883039064575464497400E1Q,
606
  2.458414064785315978408974662900438351782E1Q,
607
  2.367237826273668567199042088835448715228E1Q,
608
  9.231451197519171090875569102116321676763E0Q,
609
 /* 1.000000000000000000000000000000000000000E0 */
610
};
611
 
612
/* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
613
   Q1(x) = 1/x (.375 + 1/x^2 R(1/x^2)),
614
   Peak relative error 1.4e-36
615
   0.4375 <= 1/x <= 0.5  */
616
#define NQ2_2r3N 9
617
static const __float128 Q2_2r3N[NQ2_2r3N + 1] = {
618
  -5.552507516089087822166822364590806076174E-6Q,
619
  -4.135067659799500521040944087433752970297E-4Q,
620
  -1.059928728869218962607068840646564457980E-2Q,
621
  -1.212070036005832342565792241385459023801E-1Q,
622
  -6.688350110633603958684302153362735625156E-1Q,
623
  -1.793587878197360221340277951304429821582E0Q,
624
  -2.225407682237197485644647380483725045326E0Q,
625
  -1.123402135458940189438898496348239744403E0Q,
626
  -1.679187241566347077204805190763597299805E-1Q,
627
  -1.458550613639093752909985189067233504148E-3Q,
628
};
629
#define NQ2_2r3D 8
630
static const __float128 Q2_2r3D[NQ2_2r3D + 1] = {
631
  5.415024336507980465169023996403597916115E-5Q,
632
  4.179246497380453022046357404266022870788E-3Q,
633
  1.136306384261959483095442402929502368598E-1Q,
634
  1.422640343719842213484515445393284072830E0Q,
635
  8.968786703393158374728850922289204805764E0Q,
636
  2.914542473339246127533384118781216495934E1Q,
637
  4.781605421020380669870197378210457054685E1Q,
638
  3.693865837171883152382820584714795072937E1Q,
639
  1.153220502744204904763115556224395893076E1Q,
640
  /* 1.000000000000000000000000000000000000000E0 */
641
};
642
 
643
 
644
/* Evaluate P[n] x^n  +  P[n-1] x^(n-1)  +  ...  +  P[0] */
645
 
646
static __float128
647
neval (__float128 x, const __float128 *p, int n)
648
{
649
  __float128 y;
650
 
651
  p += n;
652
  y = *p--;
653
  do
654
    {
655
      y = y * x + *p--;
656
    }
657
  while (--n > 0);
658
  return y;
659
}
660
 
661
 
662
/* Evaluate x^n+1  +  P[n] x^(n)  +  P[n-1] x^(n-1)  +  ...  +  P[0] */
663
 
664
static __float128
665
deval (__float128 x, const __float128 *p, int n)
666
{
667
  __float128 y;
668
 
669
  p += n;
670
  y = x + *p--;
671
  do
672
    {
673
      y = y * x + *p--;
674
    }
675
  while (--n > 0);
676
  return y;
677
}
678
 
679
 
680
/* Bessel function of the first kind, order one.  */
681
 
682
__float128
683
j1q (__float128 x)
684
{
685
  __float128 xx, xinv, z, p, q, c, s, cc, ss;
686
 
687
  if (! finiteq (x))
688
    {
689
      if (x != x)
690
        return x;
691
      else
692
        return 0.0Q;
693
    }
694
  if (x == 0.0Q)
695
    return x;
696
  xx = fabsq (x);
697
  if (xx <= 2.0Q)
698
    {
699
      /* 0 <= x <= 2 */
700
      z = xx * xx;
701
      p = xx * z * neval (z, J0_2N, NJ0_2N) / deval (z, J0_2D, NJ0_2D);
702
      p += 0.5Q * xx;
703
      if (x < 0)
704
        p = -p;
705
      return p;
706
    }
707
 
708
  xinv = 1.0Q / xx;
709
  z = xinv * xinv;
710
  if (xinv <= 0.25)
711
    {
712
      if (xinv <= 0.125)
713
        {
714
          if (xinv <= 0.0625)
715
            {
716
              p = neval (z, P16_IN, NP16_IN) / deval (z, P16_ID, NP16_ID);
717
              q = neval (z, Q16_IN, NQ16_IN) / deval (z, Q16_ID, NQ16_ID);
718
            }
719
          else
720
            {
721
              p = neval (z, P8_16N, NP8_16N) / deval (z, P8_16D, NP8_16D);
722
              q = neval (z, Q8_16N, NQ8_16N) / deval (z, Q8_16D, NQ8_16D);
723
            }
724
        }
725
      else if (xinv <= 0.1875)
726
        {
727
          p = neval (z, P5_8N, NP5_8N) / deval (z, P5_8D, NP5_8D);
728
          q = neval (z, Q5_8N, NQ5_8N) / deval (z, Q5_8D, NQ5_8D);
729
        }
730
      else
731
        {
732
          p = neval (z, P4_5N, NP4_5N) / deval (z, P4_5D, NP4_5D);
733
          q = neval (z, Q4_5N, NQ4_5N) / deval (z, Q4_5D, NQ4_5D);
734
        }
735
    }                           /* .25 */
736
  else /* if (xinv <= 0.5) */
737
    {
738
      if (xinv <= 0.375)
739
        {
740
          if (xinv <= 0.3125)
741
            {
742
              p = neval (z, P3r2_4N, NP3r2_4N) / deval (z, P3r2_4D, NP3r2_4D);
743
              q = neval (z, Q3r2_4N, NQ3r2_4N) / deval (z, Q3r2_4D, NQ3r2_4D);
744
            }
745
          else
746
            {
747
              p = neval (z, P2r7_3r2N, NP2r7_3r2N)
748
                  / deval (z, P2r7_3r2D, NP2r7_3r2D);
749
              q = neval (z, Q2r7_3r2N, NQ2r7_3r2N)
750
                  / deval (z, Q2r7_3r2D, NQ2r7_3r2D);
751
            }
752
        }
753
      else if (xinv <= 0.4375)
754
        {
755
          p = neval (z, P2r3_2r7N, NP2r3_2r7N)
756
              / deval (z, P2r3_2r7D, NP2r3_2r7D);
757
          q = neval (z, Q2r3_2r7N, NQ2r3_2r7N)
758
              / deval (z, Q2r3_2r7D, NQ2r3_2r7D);
759
        }
760
      else
761
        {
762
          p = neval (z, P2_2r3N, NP2_2r3N) / deval (z, P2_2r3D, NP2_2r3D);
763
          q = neval (z, Q2_2r3N, NQ2_2r3N) / deval (z, Q2_2r3D, NQ2_2r3D);
764
        }
765
    }
766
  p = 1.0Q + z * p;
767
  q = z * q;
768
  q = q * xinv + 0.375Q * xinv;
769
  /* X = x - 3 pi/4
770
     cos(X) = cos(x) cos(3 pi/4) + sin(x) sin(3 pi/4)
771
     = 1/sqrt(2) * (-cos(x) + sin(x))
772
     sin(X) = sin(x) cos(3 pi/4) - cos(x) sin(3 pi/4)
773
     = -1/sqrt(2) * (sin(x) + cos(x))
774
     cf. Fdlibm.  */
775
  sincosq (xx, &s, &c);
776
  ss = -s - c;
777
  cc = s - c;
778
  z = cosq (xx + xx);
779
  if ((s * c) > 0)
780
    cc = z / ss;
781
  else
782
    ss = z / cc;
783
  z = ONEOSQPI * (p * cc - q * ss) / sqrtq (xx);
784
  if (x < 0)
785
    z = -z;
786
  return z;
787
}
788
 
789
 
790
/* Y1(x) = 2/pi * (log(x) * J1(x) - 1/x) + x R(x^2)
791
   Peak relative error 6.2e-38
792
 
793
#define NY0_2N 7
794
static __float128 Y0_2N[NY0_2N + 1] = {
795
  -6.804415404830253804408698161694720833249E19Q,
796
  1.805450517967019908027153056150465849237E19Q,
797
  -8.065747497063694098810419456383006737312E17Q,
798
  1.401336667383028259295830955439028236299E16Q,
799
  -1.171654432898137585000399489686629680230E14Q,
800
  5.061267920943853732895341125243428129150E11Q,
801
  -1.096677850566094204586208610960870217970E9Q,
802
  9.541172044989995856117187515882879304461E5Q,
803
};
804
#define NY0_2D 7
805
static __float128 Y0_2D[NY0_2D + 1] = {
806
  3.470629591820267059538637461549677594549E20Q,
807
  4.120796439009916326855848107545425217219E18Q,
808
  2.477653371652018249749350657387030814542E16Q,
809
  9.954678543353888958177169349272167762797E13Q,
810
  2.957927997613630118216218290262851197754E11Q,
811
  6.748421382188864486018861197614025972118E8Q,
812
  1.173453425218010888004562071020305709319E6Q,
813
  1.450335662961034949894009554536003377187E3Q,
814
  /* 1.000000000000000000000000000000000000000E0 */
815
};
816
 
817
 
818
/* Bessel function of the second kind, order one.  */
819
 
820
__float128
821
y1q (__float128 x)
822
{
823
  __float128 xx, xinv, z, p, q, c, s, cc, ss;
824
 
825
  if (! finiteq (x))
826
    {
827
      if (x != x)
828
        return x;
829
      else
830
        return 0.0Q;
831
    }
832
  if (x <= 0.0Q)
833
    {
834
      if (x < 0.0Q)
835
        return (zero / (zero * x));
836
      return -HUGE_VALQ + x;
837
    }
838
  xx = fabsq (x);
839
  if (xx <= 2.0Q)
840
    {
841
      /* 0 <= x <= 2 */
842
      z = xx * xx;
843
      p = xx * neval (z, Y0_2N, NY0_2N) / deval (z, Y0_2D, NY0_2D);
844
      p = -TWOOPI / xx + p;
845
      p = TWOOPI * logq (x) * j1q (x) + p;
846
      return p;
847
    }
848
 
849
  xinv = 1.0Q / xx;
850
  z = xinv * xinv;
851
  if (xinv <= 0.25)
852
    {
853
      if (xinv <= 0.125)
854
        {
855
          if (xinv <= 0.0625)
856
            {
857
              p = neval (z, P16_IN, NP16_IN) / deval (z, P16_ID, NP16_ID);
858
              q = neval (z, Q16_IN, NQ16_IN) / deval (z, Q16_ID, NQ16_ID);
859
            }
860
          else
861
            {
862
              p = neval (z, P8_16N, NP8_16N) / deval (z, P8_16D, NP8_16D);
863
              q = neval (z, Q8_16N, NQ8_16N) / deval (z, Q8_16D, NQ8_16D);
864
            }
865
        }
866
      else if (xinv <= 0.1875)
867
        {
868
          p = neval (z, P5_8N, NP5_8N) / deval (z, P5_8D, NP5_8D);
869
          q = neval (z, Q5_8N, NQ5_8N) / deval (z, Q5_8D, NQ5_8D);
870
        }
871
      else
872
        {
873
          p = neval (z, P4_5N, NP4_5N) / deval (z, P4_5D, NP4_5D);
874
          q = neval (z, Q4_5N, NQ4_5N) / deval (z, Q4_5D, NQ4_5D);
875
        }
876
    }                           /* .25 */
877
  else /* if (xinv <= 0.5) */
878
    {
879
      if (xinv <= 0.375)
880
        {
881
          if (xinv <= 0.3125)
882
            {
883
              p = neval (z, P3r2_4N, NP3r2_4N) / deval (z, P3r2_4D, NP3r2_4D);
884
              q = neval (z, Q3r2_4N, NQ3r2_4N) / deval (z, Q3r2_4D, NQ3r2_4D);
885
            }
886
          else
887
            {
888
              p = neval (z, P2r7_3r2N, NP2r7_3r2N)
889
                  / deval (z, P2r7_3r2D, NP2r7_3r2D);
890
              q = neval (z, Q2r7_3r2N, NQ2r7_3r2N)
891
                  / deval (z, Q2r7_3r2D, NQ2r7_3r2D);
892
            }
893
        }
894
      else if (xinv <= 0.4375)
895
        {
896
          p = neval (z, P2r3_2r7N, NP2r3_2r7N)
897
              / deval (z, P2r3_2r7D, NP2r3_2r7D);
898
          q = neval (z, Q2r3_2r7N, NQ2r3_2r7N)
899
              / deval (z, Q2r3_2r7D, NQ2r3_2r7D);
900
        }
901
      else
902
        {
903
          p = neval (z, P2_2r3N, NP2_2r3N) / deval (z, P2_2r3D, NP2_2r3D);
904
          q = neval (z, Q2_2r3N, NQ2_2r3N) / deval (z, Q2_2r3D, NQ2_2r3D);
905
        }
906
    }
907
  p = 1.0Q + z * p;
908
  q = z * q;
909
  q = q * xinv + 0.375Q * xinv;
910
  /* X = x - 3 pi/4
911
     cos(X) = cos(x) cos(3 pi/4) + sin(x) sin(3 pi/4)
912
     = 1/sqrt(2) * (-cos(x) + sin(x))
913
     sin(X) = sin(x) cos(3 pi/4) - cos(x) sin(3 pi/4)
914
     = -1/sqrt(2) * (sin(x) + cos(x))
915
     cf. Fdlibm.  */
916
  sincosq (xx, &s, &c);
917
  ss = -s - c;
918
  cc = s - c;
919
  z = cosq (xx + xx);
920
  if ((s * c) > 0)
921
    cc = z / ss;
922
  else
923
    ss = z / cc;
924
  z = ONEOSQPI * (p * ss + q * cc) / sqrtq (xx);
925
  return z;
926
}

powered by: WebSVN 2.1.0

© copyright 1999-2024 OpenCores.org, equivalent to Oliscience, all rights reserved. OpenCores®, registered trademark.