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

Subversion Repositories openrisc

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

Go to most recent revision | Details | Compare with Previous | View Log

Line No. Rev Author Line
1 740 jeremybenn 
/*                                                      j0l.c
2 
 *
3 
 *      Bessel function of order zero
4 
 *
5 
 *
6 
 *
7 
 * SYNOPSIS:
8 
 *
9 
 * long double x, y, j0l();
10 
 *
11 
 * y = j0l( x );
12 
 *
13 
 *
14 
 *
15 
 * DESCRIPTION:
16 
 *
17 
 * Returns Bessel function of first kind, order zero 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
21 
 * is J0(x) = 1 - x^2 / 4 + x^4 R(x^2)
22 
 * The second interval is further partitioned into eight equal segments
23 
 * of 1/x.
24 
 *
25 
 * J0(x) = sqrt(2/(pi x)) (P0(x) cos(X) - Q0(x) sin(X)),
26 
 * X = x - pi/4,
27 
 *
28 
 * and the auxiliary functions are given by
29 
 *
30 
 * J0(x)cos(X) + Y0(x)sin(X) = sqrt( 2/(pi x)) P0(x),
31 
 * P0(x) = 1 + 1/x^2 R(1/x^2)
32 
 *
33 
 * Y0(x)cos(X) - J0(x)sin(X) = sqrt( 2/(pi x)) Q0(x),
34 
 * Q0(x) = 1/x (-.125 + 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      1.7e-34      2.4e-35
43 
 *
44 
 *
45 
 */
46 
 
47 
/*                                                      y0l.c
48 
 *
49 
 *      Bessel function of the second kind, order zero
50 
 *
51 
 *
52 
 *
53 
 * SYNOPSIS:
54 
 *
55 
 * double x, y, y0l();
56 
 *
57 
 * y = y0l( x );
58 
 *
59 
 *
60 
 *
61 
 * DESCRIPTION:
62 
 *
63 
 * Returns Bessel function of the second kind, of order
64 
 * zero, of the argument.
65 
 *
66 
 * The approximation is the same as for J0(x), and
67 
 * Y0(x) = sqrt(2/(pi x)) (P0(x) sin(X) + Q0(x) cos(X)).
68 
 *
69 
 * ACCURACY:
70 
 *
71 
 *  Absolute error, when y0(x) < 1; else relative error:
72 
 *
73 
 * arithmetic   domain     # trials      peak         rms
74 
 *    IEEE      0, 30       100000      3.0e-34     2.7e-35
75 
 *
76 
 */
77 
 
78 
/* Copyright 2001 by Stephen L. Moshier (moshier@na-net.ornl.gov).
79 
 
80 
    This library is free software; you can redistribute it and/or
81 
    modify it under the terms of the GNU Lesser General Public
82 
    License as published by the Free Software Foundation; either
83 
    version 2.1 of the License, or (at your option) any later version.
84 
 
85 
    This library is distributed in the hope that it will be useful,
86 
    but WITHOUT ANY WARRANTY; without even the implied warranty of
87 
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
88 
    Lesser General Public License for more details.
89 
 
90 
    You should have received a copy of the GNU Lesser General Public
91 
    License along with this library; if not, write to the Free Software
92 
    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307  USA */
93 
 
94 
#include "quadmath-imp.h"
95 
 
96 
/* 1 / sqrt(pi) */
97 
static const __float128 ONEOSQPI = 5.6418958354775628694807945156077258584405E-1Q;
98 
/* 2 / pi */
99 
static const __float128 TWOOPI = 6.3661977236758134307553505349005744813784E-1Q;
100 
static const __float128 zero = 0.0Q;
101 
 
102 
/* J0(x) = 1 - x^2/4 + x^2 x^2 R(x^2)
103 
   Peak relative error 3.4e-37
104 
 
105 
#define NJ0_2N 6
106 
static const __float128 J0_2N[NJ0_2N + 1] = {
107 
  3.133239376997663645548490085151484674892E16Q,
108 
 -5.479944965767990821079467311839107722107E14Q,
109 
  6.290828903904724265980249871997551894090E12Q,
110 
 -3.633750176832769659849028554429106299915E10Q,
111 
  1.207743757532429576399485415069244807022E8Q,
112 
 -2.107485999925074577174305650549367415465E5Q,
113 
  1.562826808020631846245296572935547005859E2Q,
114 
};
115 
#define NJ0_2D 6
116 
static const __float128 J0_2D[NJ0_2D + 1] = {
117 
  2.005273201278504733151033654496928968261E18Q,
118 
  2.063038558793221244373123294054149790864E16Q,
119 
  1.053350447931127971406896594022010524994E14Q,
120 
  3.496556557558702583143527876385508882310E11Q,
121 
  8.249114511878616075860654484367133976306E8Q,
122 
  1.402965782449571800199759247964242790589E6Q,
123 
  1.619910762853439600957801751815074787351E3Q,
124 
 /* 1.000000000000000000000000000000000000000E0 */
125 
};
126 
 
127 
/* J0(x)cosX + Y0(x)sinX = sqrt( 2/(pi x)) P0(x), P0(x) = 1 + 1/x^2 R(1/x^2),
128 
 
129 
   Peak relative error 3.3e-36  */
130 
#define NP16_IN 9
131 
static const __float128 P16_IN[NP16_IN + 1] = {
132 
  -1.901689868258117463979611259731176301065E-16Q,
133 
  -1.798743043824071514483008340803573980931E-13Q,
134 
  -6.481746687115262291873324132944647438959E-11Q,
135 
  -1.150651553745409037257197798528294248012E-8Q,
136 
  -1.088408467297401082271185599507222695995E-6Q,
137 
  -5.551996725183495852661022587879817546508E-5Q,
138 
  -1.477286941214245433866838787454880214736E-3Q,
139 
  -1.882877976157714592017345347609200402472E-2Q,
140 
  -9.620983176855405325086530374317855880515E-2Q,
141 
  -1.271468546258855781530458854476627766233E-1Q,
142 
};
143 
#define NP16_ID 9
144 
static const __float128 P16_ID[NP16_ID + 1] = {
145 
  2.704625590411544837659891569420764475007E-15Q,
146 
  2.562526347676857624104306349421985403573E-12Q,
147 
  9.259137589952741054108665570122085036246E-10Q,
148 
  1.651044705794378365237454962653430805272E-7Q,
149 
  1.573561544138733044977714063100859136660E-5Q,
150 
  8.134482112334882274688298469629884804056E-4Q,
151 
  2.219259239404080863919375103673593571689E-2Q,
152 
  2.976990606226596289580242451096393862792E-1Q,
153 
  1.713895630454693931742734911930937246254E0Q,
154 
  3.231552290717904041465898249160757368855E0Q,
155 
  /* 1.000000000000000000000000000000000000000E0 */
156 
};
157 
 
158 
/* J0(x)cosX + Y0(x)sinX = sqrt( 2/(pi x)) P0(x), P0(x) = 1 + 1/x^2 R(1/x^2)
159 
    0.0625 <= 1/x <= 0.125
160 
    Peak relative error 2.4e-35  */
161 
#define NP8_16N 10
162 
static const __float128 P8_16N[NP8_16N + 1] = {
163 
  -2.335166846111159458466553806683579003632E-15Q,
164 
  -1.382763674252402720401020004169367089975E-12Q,
165 
  -3.192160804534716696058987967592784857907E-10Q,
166 
  -3.744199606283752333686144670572632116899E-8Q,
167 
  -2.439161236879511162078619292571922772224E-6Q,
168 
  -9.068436986859420951664151060267045346549E-5Q,
169 
  -1.905407090637058116299757292660002697359E-3Q,
170 
  -2.164456143936718388053842376884252978872E-2Q,
171 
  -1.212178415116411222341491717748696499966E-1Q,
172 
  -2.782433626588541494473277445959593334494E-1Q,
173 
  -1.670703190068873186016102289227646035035E-1Q,
174 
};
175 
#define NP8_16D 10
176 
static const __float128 P8_16D[NP8_16D + 1] = {
177 
  3.321126181135871232648331450082662856743E-14Q,
178 
  1.971894594837650840586859228510007703641E-11Q,
179 
  4.571144364787008285981633719513897281690E-9Q,
180 
  5.396419143536287457142904742849052402103E-7Q,
181 
  3.551548222385845912370226756036899901549E-5Q,
182 
  1.342353874566932014705609788054598013516E-3Q,
183 
  2.899133293006771317589357444614157734385E-2Q,
184 
  3.455374978185770197704507681491574261545E-1Q,
185 
  2.116616964297512311314454834712634820514E0Q,
186 
  5.850768316827915470087758636881584174432E0Q,
187 
  5.655273858938766830855753983631132928968E0Q,
188 
  /* 1.000000000000000000000000000000000000000E0 */
189 
};
190 
 
191 
/* J0(x)cosX + Y0(x)sinX = sqrt( 2/(pi x)) P0(x), P0(x) = 1 + 1/x^2 R(1/x^2)
192 
  0.125 <= 1/x <= 0.1875
193 
  Peak relative error 2.7e-35  */
194 
#define NP5_8N 10
195 
static const __float128 P5_8N[NP5_8N + 1] = {
196 
  -1.270478335089770355749591358934012019596E-12Q,
197 
  -4.007588712145412921057254992155810347245E-10Q,
198 
  -4.815187822989597568124520080486652009281E-8Q,
199 
  -2.867070063972764880024598300408284868021E-6Q,
200 
  -9.218742195161302204046454768106063638006E-5Q,
201 
  -1.635746821447052827526320629828043529997E-3Q,
202 
  -1.570376886640308408247709616497261011707E-2Q,
203 
  -7.656484795303305596941813361786219477807E-2Q,
204 
  -1.659371030767513274944805479908858628053E-1Q,
205 
  -1.185340550030955660015841796219919804915E-1Q,
206 
  -8.920026499909994671248893388013790366712E-3Q,
207 
};
208 
#define NP5_8D 9
209 
static const __float128 P5_8D[NP5_8D + 1] = {
210 
  1.806902521016705225778045904631543990314E-11Q,
211 
  5.728502760243502431663549179135868966031E-9Q,
212 
  6.938168504826004255287618819550667978450E-7Q,
213 
  4.183769964807453250763325026573037785902E-5Q,
214 
  1.372660678476925468014882230851637878587E-3Q,
215 
  2.516452105242920335873286419212708961771E-2Q,
216 
  2.550502712902647803796267951846557316182E-1Q,
217 
  1.365861559418983216913629123778747617072E0Q,
218 
  3.523825618308783966723472468855042541407E0Q,
219 
  3.656365803506136165615111349150536282434E0Q,
220 
  /* 1.000000000000000000000000000000000000000E0 */
221 
};
222 
 
223 
/* J0(x)cosX + Y0(x)sinX = sqrt( 2/(pi x)) P0(x), P0(x) = 1 + 1/x^2 R(1/x^2)
224 
   Peak relative error 3.5e-35
225 
   0.1875 <= 1/x <= 0.25  */
226 
#define NP4_5N 9
227 
static const __float128 P4_5N[NP4_5N + 1] = {
228 
  -9.791405771694098960254468859195175708252E-10Q,
229 
  -1.917193059944531970421626610188102836352E-7Q,
230 
  -1.393597539508855262243816152893982002084E-5Q,
231 
  -4.881863490846771259880606911667479860077E-4Q,
232 
  -8.946571245022470127331892085881699269853E-3Q,
233 
  -8.707474232568097513415336886103899434251E-2Q,
234 
  -4.362042697474650737898551272505525973766E-1Q,
235 
  -1.032712171267523975431451359962375617386E0Q,
236 
  -9.630502683169895107062182070514713702346E-1Q,
237 
  -2.251804386252969656586810309252357233320E-1Q,
238 
};
239 
#define NP4_5D 9
240 
static const __float128 P4_5D[NP4_5D + 1] = {
241 
  1.392555487577717669739688337895791213139E-8Q,
242 
  2.748886559120659027172816051276451376854E-6Q,
243 
  2.024717710644378047477189849678576659290E-4Q,
244 
  7.244868609350416002930624752604670292469E-3Q,
245 
  1.373631762292244371102989739300382152416E-1Q,
246 
  1.412298581400224267910294815260613240668E0Q,
247 
  7.742495637843445079276397723849017617210E0Q,
248 
  2.138429269198406512028307045259503811861E1Q,
249 
  2.651547684548423476506826951831712762610E1Q,
250 
  1.167499382465291931571685222882909166935E1Q,
251 
  /* 1.000000000000000000000000000000000000000E0 */
252 
};
253 
 
254 
/* J0(x)cosX + Y0(x)sinX = sqrt( 2/(pi x)) P0(x), P0(x) = 1 + 1/x^2 R(1/x^2)
255 
   Peak relative error 2.3e-36
256 
   0.25 <= 1/x <= 0.3125  */
257 
#define NP3r2_4N 9
258 
static const __float128 P3r2_4N[NP3r2_4N + 1] = {
259 
  -2.589155123706348361249809342508270121788E-8Q,
260 
  -3.746254369796115441118148490849195516593E-6Q,
261 
  -1.985595497390808544622893738135529701062E-4Q,
262 
  -5.008253705202932091290132760394976551426E-3Q,
263 
  -6.529469780539591572179155511840853077232E-2Q,
264 
  -4.468736064761814602927408833818990271514E-1Q,
265 
  -1.556391252586395038089729428444444823380E0Q,
266 
  -2.533135309840530224072920725976994981638E0Q,
267 
  -1.605509621731068453869408718565392869560E0Q,
268 
  -2.518966692256192789269859830255724429375E-1Q,
269 
};
270 
#define NP3r2_4D 9
271 
static const __float128 P3r2_4D[NP3r2_4D + 1] = {
272 
  3.682353957237979993646169732962573930237E-7Q,
273 
  5.386741661883067824698973455566332102029E-5Q,
274 
  2.906881154171822780345134853794241037053E-3Q,
275 
  7.545832595801289519475806339863492074126E-2Q,
276 
  1.029405357245594877344360389469584526654E0Q,
277 
  7.565706120589873131187989560509757626725E0Q,
278 
  2.951172890699569545357692207898667665796E1Q,
279 
  5.785723537170311456298467310529815457536E1Q,
280 
  5.095621464598267889126015412522773474467E1Q,
281 
  1.602958484169953109437547474953308401442E1Q,
282 
  /* 1.000000000000000000000000000000000000000E0 */
283 
};
284 
 
285 
/* J0(x)cosX + Y0(x)sinX = sqrt( 2/(pi x)) P0(x), P0(x) = 1 + 1/x^2 R(1/x^2)
286 
   Peak relative error 1.0e-35
287 
   0.3125 <= 1/x <= 0.375  */
288 
#define NP2r7_3r2N 9
289 
static const __float128 P2r7_3r2N[NP2r7_3r2N + 1] = {
290 
  -1.917322340814391131073820537027234322550E-7Q,
291 
  -1.966595744473227183846019639723259011906E-5Q,
292 
  -7.177081163619679403212623526632690465290E-4Q,
293 
  -1.206467373860974695661544653741899755695E-2Q,
294 
  -1.008656452188539812154551482286328107316E-1Q,
295 
  -4.216016116408810856620947307438823892707E-1Q,
296 
  -8.378631013025721741744285026537009814161E-1Q,
297 
  -6.973895635309960850033762745957946272579E-1Q,
298 
  -1.797864718878320770670740413285763554812E-1Q,
299 
  -4.098025357743657347681137871388402849581E-3Q,
300 
};
301 
#define NP2r7_3r2D 8
302 
static const __float128 P2r7_3r2D[NP2r7_3r2D + 1] = {
303 
  2.726858489303036441686496086962545034018E-6Q,
304 
  2.840430827557109238386808968234848081424E-4Q,
305 
  1.063826772041781947891481054529454088832E-2Q,
306 
  1.864775537138364773178044431045514405468E-1Q,
307 
  1.665660052857205170440952607701728254211E0Q,
308 
  7.723745889544331153080842168958348568395E0Q,
309 
  1.810726427571829798856428548102077799835E1Q,
310 
  1.986460672157794440666187503833545388527E1Q,
311 
  8.645503204552282306364296517220055815488E0Q,
312 
  /* 1.000000000000000000000000000000000000000E0 */
313 
};
314 
 
315 
/* J0(x)cosX + Y0(x)sinX = sqrt( 2/(pi x)) P0(x), P0(x) = 1 + 1/x^2 R(1/x^2)
316 
   Peak relative error 1.3e-36
317 
   0.3125 <= 1/x <= 0.4375  */
318 
#define NP2r3_2r7N 9
319 
static const __float128 P2r3_2r7N[NP2r3_2r7N + 1] = {
320 
  -1.594642785584856746358609622003310312622E-6Q,
321 
  -1.323238196302221554194031733595194539794E-4Q,
322 
  -3.856087818696874802689922536987100372345E-3Q,
323 
  -5.113241710697777193011470733601522047399E-2Q,
324 
  -3.334229537209911914449990372942022350558E-1Q,
325 
  -1.075703518198127096179198549659283422832E0Q,
326 
  -1.634174803414062725476343124267110981807E0Q,
327 
  -1.030133247434119595616826842367268304880E0Q,
328 
  -1.989811539080358501229347481000707289391E-1Q,
329 
  -3.246859189246653459359775001466924610236E-3Q,
330 
};
331 
#define NP2r3_2r7D 8
332 
static const __float128 P2r3_2r7D[NP2r3_2r7D + 1] = {
333 
  2.267936634217251403663034189684284173018E-5Q,
334 
  1.918112982168673386858072491437971732237E-3Q,
335 
  5.771704085468423159125856786653868219522E-2Q,
336 
  8.056124451167969333717642810661498890507E-1Q,
337 
  5.687897967531010276788680634413789328776E0Q,
338 
  2.072596760717695491085444438270778394421E1Q,
339 
  3.801722099819929988585197088613160496684E1Q,
340 
  3.254620235902912339534998592085115836829E1Q,
341 
  1.104847772130720331801884344645060675036E1Q,
342 
  /* 1.000000000000000000000000000000000000000E0 */
343 
};
344 
 
345 
/* J0(x)cosX + Y0(x)sinX = sqrt( 2/(pi x)) P0(x), P0(x) = 1 + 1/x^2 R(1/x^2)
346 
   Peak relative error 1.2e-35
347 
   0.4375 <= 1/x <= 0.5  */
348 
#define NP2_2r3N 8
349 
static const __float128 P2_2r3N[NP2_2r3N + 1] = {
350 
  -1.001042324337684297465071506097365389123E-4Q,
351 
  -6.289034524673365824853547252689991418981E-3Q,
352 
  -1.346527918018624234373664526930736205806E-1Q,
353 
  -1.268808313614288355444506172560463315102E0Q,
354 
  -5.654126123607146048354132115649177406163E0Q,
355 
  -1.186649511267312652171775803270911971693E1Q,
356 
  -1.094032424931998612551588246779200724257E1Q,
357 
  -3.728792136814520055025256353193674625267E0Q,
358 
  -3.000348318524471807839934764596331810608E-1Q,
359 
};
360 
#define NP2_2r3D 8
361 
static const __float128 P2_2r3D[NP2_2r3D + 1] = {
362 
  1.423705538269770974803901422532055612980E-3Q,
363 
  9.171476630091439978533535167485230575894E-2Q,
364 
  2.049776318166637248868444600215942828537E0Q,
365 
  2.068970329743769804547326701946144899583E1Q,
366 
  1.025103500560831035592731539565060347709E2Q,
367 
  2.528088049697570728252145557167066708284E2Q,
368 
  2.992160327587558573740271294804830114205E2Q,
369 
  1.540193761146551025832707739468679973036E2Q,
370 
  2.779516701986912132637672140709452502650E1Q,
371 
  /* 1.000000000000000000000000000000000000000E0 */
372 
};
373 
 
374 
/* Y0(x)cosX - J0(x)sinX = sqrt( 2/(pi x)) Q0(x),
375 
   Q0(x) = 1/x (-.125 + 1/x^2 R(1/x^2))
376 
   Peak relative error 2.2e-35
377 
 
378 
#define NQ16_IN 10
379 
static const __float128 Q16_IN[NQ16_IN + 1] = {
380 
  2.343640834407975740545326632205999437469E-18Q,
381 
  2.667978112927811452221176781536278257448E-15Q,
382 
  1.178415018484555397390098879501969116536E-12Q,
383 
  2.622049767502719728905924701288614016597E-10Q,
384 
  3.196908059607618864801313380896308968673E-8Q,
385 
  2.179466154171673958770030655199434798494E-6Q,
386 
  8.139959091628545225221976413795645177291E-5Q,
387 
  1.563900725721039825236927137885747138654E-3Q,
388 
  1.355172364265825167113562519307194840307E-2Q,
389 
  3.928058355906967977269780046844768588532E-2Q,
390 
  1.107891967702173292405380993183694932208E-2Q,
391 
};
392 
#define NQ16_ID 9
393 
static const __float128 Q16_ID[NQ16_ID + 1] = {
394 
  3.199850952578356211091219295199301766718E-17Q,
395 
  3.652601488020654842194486058637953363918E-14Q,
396 
  1.620179741394865258354608590461839031281E-11Q,
397 
  3.629359209474609630056463248923684371426E-9Q,
398 
  4.473680923894354600193264347733477363305E-7Q,
399 
  3.106368086644715743265603656011050476736E-5Q,
400 
  1.198239259946770604954664925153424252622E-3Q,
401 
  2.446041004004283102372887804475767568272E-2Q,
402 
  2.403235525011860603014707768815113698768E-1Q,
403 
  9.491006790682158612266270665136910927149E-1Q,
404 
 /* 1.000000000000000000000000000000000000000E0 */
405 
 };
406 
 
407 
/* Y0(x)cosX - J0(x)sinX = sqrt( 2/(pi x)) Q0(x),
408 
   Q0(x) = 1/x (-.125 + 1/x^2 R(1/x^2))
409 
   Peak relative error 5.1e-36
410 
   0.0625 <= 1/x <= 0.125  */
411 
#define NQ8_16N 11
412 
static const __float128 Q8_16N[NQ8_16N + 1] = {
413 
  1.001954266485599464105669390693597125904E-17Q,
414 
  7.545499865295034556206475956620160007849E-15Q,
415 
  2.267838684785673931024792538193202559922E-12Q,
416 
  3.561909705814420373609574999542459912419E-10Q,
417 
  3.216201422768092505214730633842924944671E-8Q,
418 
  1.731194793857907454569364622452058554314E-6Q,
419 
  5.576944613034537050396518509871004586039E-5Q,
420 
  1.051787760316848982655967052985391418146E-3Q,
421 
  1.102852974036687441600678598019883746959E-2Q,
422 
  5.834647019292460494254225988766702933571E-2Q,
423 
  1.290281921604364618912425380717127576529E-1Q,
424 
  7.598886310387075708640370806458926458301E-2Q,
425 
};
426 
#define NQ8_16D 11
427 
static const __float128 Q8_16D[NQ8_16D + 1] = {
428 
  1.368001558508338469503329967729951830843E-16Q,
429 
  1.034454121857542147020549303317348297289E-13Q,
430 
  3.128109209247090744354764050629381674436E-11Q,
431 
  4.957795214328501986562102573522064468671E-9Q,
432 
  4.537872468606711261992676606899273588899E-7Q,
433 
  2.493639207101727713192687060517509774182E-5Q,
434 
  8.294957278145328349785532236663051405805E-4Q,
435 
  1.646471258966713577374948205279380115839E-2Q,
436 
  1.878910092770966718491814497982191447073E-1Q,
437 
  1.152641605706170353727903052525652504075E0Q,
438 
  3.383550240669773485412333679367792932235E0Q,
439 
  3.823875252882035706910024716609908473970E0Q,
440 
 /* 1.000000000000000000000000000000000000000E0 */
441 
};
442 
 
443 
/* Y0(x)cosX - J0(x)sinX = sqrt( 2/(pi x)) Q0(x),
444 
   Q0(x) = 1/x (-.125 + 1/x^2 R(1/x^2))
445 
   Peak relative error 3.9e-35
446 
   0.125 <= 1/x <= 0.1875  */
447 
#define NQ5_8N 10
448 
static const __float128 Q5_8N[NQ5_8N + 1] = {
449 
  1.750399094021293722243426623211733898747E-13Q,
450 
  6.483426211748008735242909236490115050294E-11Q,
451 
  9.279430665656575457141747875716899958373E-9Q,
452 
  6.696634968526907231258534757736576340266E-7Q,
453 
  2.666560823798895649685231292142838188061E-5Q,
454 
  6.025087697259436271271562769707550594540E-4Q,
455 
  7.652807734168613251901945778921336353485E-3Q,
456 
  5.226269002589406461622551452343519078905E-2Q,
457 
  1.748390159751117658969324896330142895079E-1Q,
458 
  2.378188719097006494782174902213083589660E-1Q,
459 
  8.383984859679804095463699702165659216831E-2Q,
460 
};
461 
#define NQ5_8D 10
462 
static const __float128 Q5_8D[NQ5_8D + 1] = {
463 
  2.389878229704327939008104855942987615715E-12Q,
464 
  8.926142817142546018703814194987786425099E-10Q,
465 
  1.294065862406745901206588525833274399038E-7Q,
466 
  9.524139899457666250828752185212769682191E-6Q,
467 
  3.908332488377770886091936221573123353489E-4Q,
468 
  9.250427033957236609624199884089916836748E-3Q,
469 
  1.263420066165922645975830877751588421451E-1Q,
470 
  9.692527053860420229711317379861733180654E-1Q,
471 
  3.937813834630430172221329298841520707954E0Q,
472 
  7.603126427436356534498908111445191312181E0Q,
473 
  5.670677653334105479259958485084550934305E0Q,
474 
 /* 1.000000000000000000000000000000000000000E0 */
475 
};
476 
 
477 
/* Y0(x)cosX - J0(x)sinX = sqrt( 2/(pi x)) Q0(x),
478 
   Q0(x) = 1/x (-.125 + 1/x^2 R(1/x^2))
479 
   Peak relative error 3.2e-35
480 
   0.1875 <= 1/x <= 0.25  */
481 
#define NQ4_5N 10
482 
static const __float128 Q4_5N[NQ4_5N + 1] = {
483 
  2.233870042925895644234072357400122854086E-11Q,
484 
  5.146223225761993222808463878999151699792E-9Q,
485 
  4.459114531468296461688753521109797474523E-7Q,
486 
  1.891397692931537975547242165291668056276E-5Q,
487 
  4.279519145911541776938964806470674565504E-4Q,
488 
  5.275239415656560634702073291768904783989E-3Q,
489 
  3.468698403240744801278238473898432608887E-2Q,
490 
  1.138773146337708415188856882915457888274E-1Q,
491 
  1.622717518946443013587108598334636458955E-1Q,
492 
  7.249040006390586123760992346453034628227E-2Q,
493 
  1.941595365256460232175236758506411486667E-3Q,
494 
};
495 
#define NQ4_5D 9
496 
static const __float128 Q4_5D[NQ4_5D + 1] = {
497 
  3.049977232266999249626430127217988047453E-10Q,
498 
  7.120883230531035857746096928889676144099E-8Q,
499 
  6.301786064753734446784637919554359588859E-6Q,
500 
  2.762010530095069598480766869426308077192E-4Q,
501 
  6.572163250572867859316828886203406361251E-3Q,
502 
  8.752566114841221958200215255461843397776E-2Q,
503 
  6.487654992874805093499285311075289932664E-1Q,
504 
  2.576550017826654579451615283022812801435E0Q,
505 
  5.056392229924022835364779562707348096036E0Q,
506 
  4.179770081068251464907531367859072157773E0Q,
507 
 /* 1.000000000000000000000000000000000000000E0 */
508 
};
509 
 
510 
/* Y0(x)cosX - J0(x)sinX = sqrt( 2/(pi x)) Q0(x),
511 
   Q0(x) = 1/x (-.125 + 1/x^2 R(1/x^2))
512 
   Peak relative error 1.4e-36
513 
   0.25 <= 1/x <= 0.3125  */
514 
#define NQ3r2_4N 10
515 
static const __float128 Q3r2_4N[NQ3r2_4N + 1] = {
516 
  6.126167301024815034423262653066023684411E-10Q,
517 
  1.043969327113173261820028225053598975128E-7Q,
518 
  6.592927270288697027757438170153763220190E-6Q,
519 
  2.009103660938497963095652951912071336730E-4Q,
520 
  3.220543385492643525985862356352195896964E-3Q,
521 
  2.774405975730545157543417650436941650990E-2Q,
522 
  1.258114008023826384487378016636555041129E-1Q,
523 
  2.811724258266902502344701449984698323860E-1Q,
524 
  2.691837665193548059322831687432415014067E-1Q,
525 
  7.949087384900985370683770525312735605034E-2Q,
526 
  1.229509543620976530030153018986910810747E-3Q,
527 
};
528 
#define NQ3r2_4D 9
529 
static const __float128 Q3r2_4D[NQ3r2_4D + 1] = {
530 
  8.364260446128475461539941389210166156568E-9Q,
531 
  1.451301850638956578622154585560759862764E-6Q,
532 
  9.431830010924603664244578867057141839463E-5Q,
533 
  3.004105101667433434196388593004526182741E-3Q,
534 
  5.148157397848271739710011717102773780221E-2Q,
535 
  4.901089301726939576055285374953887874895E-1Q,
536 
  2.581760991981709901216967665934142240346E0Q,
537 
  7.257105880775059281391729708630912791847E0Q,
538 
  1.006014717326362868007913423810737369312E1Q,
539 
  5.879416600465399514404064187445293212470E0Q,
540 
 /* 1.000000000000000000000000000000000000000E0*/
541 
};
542 
 
543 
/* Y0(x)cosX - J0(x)sinX = sqrt( 2/(pi x)) Q0(x),
544 
   Q0(x) = 1/x (-.125 + 1/x^2 R(1/x^2))
545 
   Peak relative error 3.8e-36
546 
   0.3125 <= 1/x <= 0.375  */
547 
#define NQ2r7_3r2N 9
548 
static const __float128 Q2r7_3r2N[NQ2r7_3r2N + 1] = {
549 
  7.584861620402450302063691901886141875454E-8Q,
550 
  9.300939338814216296064659459966041794591E-6Q,
551 
  4.112108906197521696032158235392604947895E-4Q,
552 
  8.515168851578898791897038357239630654431E-3Q,
553 
  8.971286321017307400142720556749573229058E-2Q,
554 
  4.885856732902956303343015636331874194498E-1Q,
555 
  1.334506268733103291656253500506406045846E0Q,
556 
  1.681207956863028164179042145803851824654E0Q,
557 
  8.165042692571721959157677701625853772271E-1Q,
558 
  9.805848115375053300608712721986235900715E-2Q,
559 
};
560 
#define NQ2r7_3r2D 9
561 
static const __float128 Q2r7_3r2D[NQ2r7_3r2D + 1] = {
562 
  1.035586492113036586458163971239438078160E-6Q,
563 
  1.301999337731768381683593636500979713689E-4Q,
564 
  5.993695702564527062553071126719088859654E-3Q,
565 
  1.321184892887881883489141186815457808785E-1Q,
566 
  1.528766555485015021144963194165165083312E0Q,
567 
  9.561463309176490874525827051566494939295E0Q,
568 
  3.203719484883967351729513662089163356911E1Q,
569 
  5.497294687660930446641539152123568668447E1Q,
570 
  4.391158169390578768508675452986948391118E1Q,
571 
  1.347836630730048077907818943625789418378E1Q,
572 
 /* 1.000000000000000000000000000000000000000E0 */
573 
};
574 
 
575 
/* Y0(x)cosX - J0(x)sinX = sqrt( 2/(pi x)) Q0(x),
576 
   Q0(x) = 1/x (-.125 + 1/x^2 R(1/x^2))
577 
   Peak relative error 2.2e-35
578 
   0.375 <= 1/x <= 0.4375  */
579 
#define NQ2r3_2r7N 9
580 
static const __float128 Q2r3_2r7N[NQ2r3_2r7N + 1] = {
581 
  4.455027774980750211349941766420190722088E-7Q,
582 
  4.031998274578520170631601850866780366466E-5Q,
583 
  1.273987274325947007856695677491340636339E-3Q,
584 
  1.818754543377448509897226554179659122873E-2Q,
585 
  1.266748858326568264126353051352269875352E-1Q,
586 
  4.327578594728723821137731555139472880414E-1Q,
587 
  6.892532471436503074928194969154192615359E-1Q,
588 
  4.490775818438716873422163588640262036506E-1Q,
589 
  8.649615949297322440032000346117031581572E-2Q,
590 
  7.261345286655345047417257611469066147561E-4Q,
591 
};
592 
#define NQ2r3_2r7D 8
593 
static const __float128 Q2r3_2r7D[NQ2r3_2r7D + 1] = {
594 
  6.082600739680555266312417978064954793142E-6Q,
595 
  5.693622538165494742945717226571441747567E-4Q,
596 
  1.901625907009092204458328768129666975975E-2Q,
597 
  2.958689532697857335456896889409923371570E-1Q,
598 
  2.343124711045660081603809437993368799568E0Q,
599 
  9.665894032187458293568704885528192804376E0Q,
600 
  2.035273104990617136065743426322454881353E1Q,
601 
  2.044102010478792896815088858740075165531E1Q,
602 
  8.445937177863155827844146643468706599304E0Q,
603 
 /* 1.000000000000000000000000000000000000000E0 */
604 
};
605 
 
606 
/* Y0(x)cosX - J0(x)sinX = sqrt( 2/(pi x)) Q0(x),
607 
   Q0(x) = 1/x (-.125 + 1/x^2 R(1/x^2))
608 
   Peak relative error 3.1e-36
609 
   0.4375 <= 1/x <= 0.5  */
610 
#define NQ2_2r3N 9
611 
static const __float128 Q2_2r3N[NQ2_2r3N + 1] = {
612 
  2.817566786579768804844367382809101929314E-6Q,
613 
  2.122772176396691634147024348373539744935E-4Q,
614 
  5.501378031780457828919593905395747517585E-3Q,
615 
  6.355374424341762686099147452020466524659E-2Q,
616 
  3.539652320122661637429658698954748337223E-1Q,
617 
  9.571721066119617436343740541777014319695E-1Q,
618 
  1.196258777828426399432550698612171955305E0Q,
619 
  6.069388659458926158392384709893753793967E-1Q,
620 
  9.026746127269713176512359976978248763621E-2Q,
621 
  5.317668723070450235320878117210807236375E-4Q,
622 
};
623 
#define NQ2_2r3D 8
624 
static const __float128 Q2_2r3D[NQ2_2r3D + 1] = {
625 
  3.846924354014260866793741072933159380158E-5Q,
626 
  3.017562820057704325510067178327449946763E-3Q,
627 
  8.356305620686867949798885808540444210935E-2Q,
628 
  1.068314930499906838814019619594424586273E0Q,
629 
  6.900279623894821067017966573640732685233E0Q,
630 
  2.307667390886377924509090271780839563141E1Q,
631 
  3.921043465412723970791036825401273528513E1Q,
632 
  3.167569478939719383241775717095729233436E1Q,
633 
  1.051023841699200920276198346301543665909E1Q,
634 
 /* 1.000000000000000000000000000000000000000E0*/
635 
};
636 
 
637 
 
638 
/* Evaluate P[n] x^n  +  P[n-1] x^(n-1)  +  ...  +  P[0] */
639 
 
640 
static __float128
641 
neval (__float128 x, const __float128 *p, int n)
642 
{
643 
  __float128 y;
644 
 
645 
  p += n;
646 
  y = *p--;
647 
  do
648 
    {
649 
      y = y * x + *p--;
650 
    }
651 
  while (--n > 0);
652 
  return y;
653 
}
654 
 
655 
 
656 
/* Evaluate x^n+1  +  P[n] x^(n)  +  P[n-1] x^(n-1)  +  ...  +  P[0] */
657 
 
658 
static __float128
659 
deval (__float128 x, const __float128 *p, int n)
660 
{
661 
  __float128 y;
662 
 
663 
  p += n;
664 
  y = x + *p--;
665 
  do
666 
    {
667 
      y = y * x + *p--;
668 
    }
669 
  while (--n > 0);
670 
  return y;
671 
}
672 
 
673 
 
674 
/* Bessel function of the first kind, order zero.  */
675 
 
676 
__float128
677 
j0q (__float128 x)
678 
{
679 
  __float128 xx, xinv, z, p, q, c, s, cc, ss;
680 
 
681 
  if (! finiteq (x))
682 
    {
683 
      if (x != x)
684 
        return x;
685 
      else
686 
        return 0.0Q;
687 
    }
688 
  if (x == 0.0Q)
689 
    return 1.0Q;
690 
 
691 
  xx = fabsq (x);
692 
  if (xx <= 2.0Q)
693 
    {
694 
      /* 0 <= x <= 2 */
695 
      z = xx * xx;
696 
      p = z * z * neval (z, J0_2N, NJ0_2N) / deval (z, J0_2D, NJ0_2D);
697 
      p -= 0.25Q * z;
698 
      p += 1.0Q;
699 
      return p;
700 
    }
701 
 
702 
  xinv = 1.0Q / xx;
703 
  z = xinv * xinv;
704 
  if (xinv <= 0.25)
705 
    {
706 
      if (xinv <= 0.125)
707 
        {
708 
          if (xinv <= 0.0625)
709 
            {
710 
              p = neval (z, P16_IN, NP16_IN) / deval (z, P16_ID, NP16_ID);
711 
              q = neval (z, Q16_IN, NQ16_IN) / deval (z, Q16_ID, NQ16_ID);
712 
            }
713 
          else
714 
            {
715 
              p = neval (z, P8_16N, NP8_16N) / deval (z, P8_16D, NP8_16D);
716 
              q = neval (z, Q8_16N, NQ8_16N) / deval (z, Q8_16D, NQ8_16D);
717 
            }
718 
        }
719 
      else if (xinv <= 0.1875)
720 
        {
721 
          p = neval (z, P5_8N, NP5_8N) / deval (z, P5_8D, NP5_8D);
722 
          q = neval (z, Q5_8N, NQ5_8N) / deval (z, Q5_8D, NQ5_8D);
723 
        }
724 
      else
725 
        {
726 
          p = neval (z, P4_5N, NP4_5N) / deval (z, P4_5D, NP4_5D);
727 
          q = neval (z, Q4_5N, NQ4_5N) / deval (z, Q4_5D, NQ4_5D);
728 
        }
729 
    }                           /* .25 */
730 
  else /* if (xinv <= 0.5) */
731 
    {
732 
      if (xinv <= 0.375)
733 
        {
734 
          if (xinv <= 0.3125)
735 
            {
736 
              p = neval (z, P3r2_4N, NP3r2_4N) / deval (z, P3r2_4D, NP3r2_4D);
737 
              q = neval (z, Q3r2_4N, NQ3r2_4N) / deval (z, Q3r2_4D, NQ3r2_4D);
738 
            }
739 
          else
740 
            {
741 
              p = neval (z, P2r7_3r2N, NP2r7_3r2N)
742 
                  / deval (z, P2r7_3r2D, NP2r7_3r2D);
743 
              q = neval (z, Q2r7_3r2N, NQ2r7_3r2N)
744 
                  / deval (z, Q2r7_3r2D, NQ2r7_3r2D);
745 
            }
746 
        }
747 
      else if (xinv <= 0.4375)
748 
        {
749 
          p = neval (z, P2r3_2r7N, NP2r3_2r7N)
750 
              / deval (z, P2r3_2r7D, NP2r3_2r7D);
751 
          q = neval (z, Q2r3_2r7N, NQ2r3_2r7N)
752 
              / deval (z, Q2r3_2r7D, NQ2r3_2r7D);
753 
        }
754 
      else
755 
        {
756 
          p = neval (z, P2_2r3N, NP2_2r3N) / deval (z, P2_2r3D, NP2_2r3D);
757 
          q = neval (z, Q2_2r3N, NQ2_2r3N) / deval (z, Q2_2r3D, NQ2_2r3D);
758 
        }
759 
    }
760 
  p = 1.0Q + z * p;
761 
  q = z * xinv * q;
762 
  q = q - 0.125Q * xinv;
763 
  /* X = x - pi/4
764 
     cos(X) = cos(x) cos(pi/4) + sin(x) sin(pi/4)
765 
     = 1/sqrt(2) * (cos(x) + sin(x))
766 
     sin(X) = sin(x) cos(pi/4) - cos(x) sin(pi/4)
767 
     = 1/sqrt(2) * (sin(x) - cos(x))
768 
     sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
769 
     cf. Fdlibm.  */
770 
  sincosq (xx, &s, &c);
771 
  ss = s - c;
772 
  cc = s + c;
773 
  z = - cosq (xx + xx);
774 
  if ((s * c) < 0)
775 
    cc = z / ss;
776 
  else
777 
    ss = z / cc;
778 
  z = ONEOSQPI * (p * cc - q * ss) / sqrtq (xx);
779 
  return z;
780 
}
781 
 
782 
 
783 
/* Y0(x) = 2/pi * log(x) * J0(x) + R(x^2)
784 
   Peak absolute error 1.7e-36 (relative where Y0 > 1)
785 
 
786 
#define NY0_2N 7
787 
static __float128 Y0_2N[NY0_2N + 1] = {
788 
 -1.062023609591350692692296993537002558155E19Q,
789 
  2.542000883190248639104127452714966858866E19Q,
790 
 -1.984190771278515324281415820316054696545E18Q,
791 
  4.982586044371592942465373274440222033891E16Q,
792 
 -5.529326354780295177243773419090123407550E14Q,
793 
  3.013431465522152289279088265336861140391E12Q,
794 
 -7.959436160727126750732203098982718347785E9Q,
795 
  8.230845651379566339707130644134372793322E6Q,
796 
};
797 
#define NY0_2D 7
798 
static __float128 Y0_2D[NY0_2D + 1] = {
799 
  1.438972634353286978700329883122253752192E20Q,
800 
  1.856409101981569254247700169486907405500E18Q,
801 
  1.219693352678218589553725579802986255614E16Q,
802 
  5.389428943282838648918475915779958097958E13Q,
803 
  1.774125762108874864433872173544743051653E11Q,
804 
  4.522104832545149534808218252434693007036E8Q,
805 
  8.872187401232943927082914504125234454930E5Q,
806 
  1.251945613186787532055610876304669413955E3Q,
807 
 /* 1.000000000000000000000000000000000000000E0 */
808 
};
809 
 
810 
 
811 
/* Bessel function of the second kind, order zero.  */
812 
 
813 
__float128
814 
y0q (__float128 x)
815 
{
816 
  __float128 xx, xinv, z, p, q, c, s, cc, ss;
817 
 
818 
  if (! finiteq (x))
819 
    {
820 
      if (x != x)
821 
        return x;
822 
      else
823 
        return 0.0Q;
824 
    }
825 
  if (x <= 0.0Q)
826 
    {
827 
      if (x < 0.0Q)
828 
        return (zero / (zero * x));
829 
      return -HUGE_VALQ + x;
830 
    }
831 
  xx = fabsq (x);
832 
  if (xx <= 2.0Q)
833 
    {
834 
      /* 0 <= x <= 2 */
835 
      z = xx * xx;
836 
      p = neval (z, Y0_2N, NY0_2N) / deval (z, Y0_2D, NY0_2D);
837 
      p = TWOOPI * logq (x) * j0q (x) + p;
838 
      return p;
839 
    }
840 
 
841 
  xinv = 1.0Q / xx;
842 
  z = xinv * xinv;
843 
  if (xinv <= 0.25)
844 
    {
845 
      if (xinv <= 0.125)
846 
        {
847 
          if (xinv <= 0.0625)
848 
            {
849 
              p = neval (z, P16_IN, NP16_IN) / deval (z, P16_ID, NP16_ID);
850 
              q = neval (z, Q16_IN, NQ16_IN) / deval (z, Q16_ID, NQ16_ID);
851 
            }
852 
          else
853 
            {
854 
              p = neval (z, P8_16N, NP8_16N) / deval (z, P8_16D, NP8_16D);
855 
              q = neval (z, Q8_16N, NQ8_16N) / deval (z, Q8_16D, NQ8_16D);
856 
            }
857 
        }
858 
      else if (xinv <= 0.1875)
859 
        {
860 
          p = neval (z, P5_8N, NP5_8N) / deval (z, P5_8D, NP5_8D);
861 
          q = neval (z, Q5_8N, NQ5_8N) / deval (z, Q5_8D, NQ5_8D);
862 
        }
863 
      else
864 
        {
865 
          p = neval (z, P4_5N, NP4_5N) / deval (z, P4_5D, NP4_5D);
866 
          q = neval (z, Q4_5N, NQ4_5N) / deval (z, Q4_5D, NQ4_5D);
867 
        }
868 
    }                           /* .25 */
869 
  else /* if (xinv <= 0.5) */
870 
    {
871 
      if (xinv <= 0.375)
872 
        {
873 
          if (xinv <= 0.3125)
874 
            {
875 
              p = neval (z, P3r2_4N, NP3r2_4N) / deval (z, P3r2_4D, NP3r2_4D);
876 
              q = neval (z, Q3r2_4N, NQ3r2_4N) / deval (z, Q3r2_4D, NQ3r2_4D);
877 
            }
878 
          else
879 
            {
880 
              p = neval (z, P2r7_3r2N, NP2r7_3r2N)
881 
                  / deval (z, P2r7_3r2D, NP2r7_3r2D);
882 
              q = neval (z, Q2r7_3r2N, NQ2r7_3r2N)
883 
                  / deval (z, Q2r7_3r2D, NQ2r7_3r2D);
884 
            }
885 
        }
886 
      else if (xinv <= 0.4375)
887 
        {
888 
          p = neval (z, P2r3_2r7N, NP2r3_2r7N)
889 
              / deval (z, P2r3_2r7D, NP2r3_2r7D);
890 
          q = neval (z, Q2r3_2r7N, NQ2r3_2r7N)
891 
              / deval (z, Q2r3_2r7D, NQ2r3_2r7D);
892 
        }
893 
      else
894 
        {
895 
          p = neval (z, P2_2r3N, NP2_2r3N) / deval (z, P2_2r3D, NP2_2r3D);
896 
          q = neval (z, Q2_2r3N, NQ2_2r3N) / deval (z, Q2_2r3D, NQ2_2r3D);
897 
        }
898 
    }
899 
  p = 1.0Q + z * p;
900 
  q = z * xinv * q;
901 
  q = q - 0.125Q * xinv;
902 
  /* X = x - pi/4
903 
     cos(X) = cos(x) cos(pi/4) + sin(x) sin(pi/4)
904 
     = 1/sqrt(2) * (cos(x) + sin(x))
905 
     sin(X) = sin(x) cos(pi/4) - cos(x) sin(pi/4)
906 
     = 1/sqrt(2) * (sin(x) - cos(x))
907 
     sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
908 
     cf. Fdlibm.  */
909 
  sincosq (x, &s, &c);
910 
  ss = s - c;
911 
  cc = s + c;
912 
  z = - cosq (x + x);
913 
  if ((s * c) < 0)
914 
    cc = z / ss;
915 
  else
916 
    ss = z / cc;
917 
  z = ONEOSQPI * (p * ss + q * cc) / sqrtq (x);
918 
  return z;
919 
}

powered by: WebSVN 2.1.0

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