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

Subversion Repositories openrisc

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

Details | Compare with Previous | View Log

Line No. Rev Author Line
1 740 jeremybenn
/*
2
 * ====================================================
3
 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
4
 *
5
 * Developed at SunPro, a Sun Microsystems, Inc. business.
6
 * Permission to use, copy, modify, and distribute this
7
 * software is freely granted, provided that this notice
8
 * is preserved.
9
 * ====================================================
10
 */
11
 
12
/* Modifications and expansions for 128-bit long double are
13
   Copyright (C) 2001 Stephen L. Moshier <moshier@na-net.ornl.gov>
14
   and are incorporated herein by permission of the author.  The author
15
   reserves the right to distribute this material elsewhere under different
16
   copying permissions.  These modifications are distributed here under
17
   the following terms:
18
 
19
    This library is free software; you can redistribute it and/or
20
    modify it under the terms of the GNU Lesser General Public
21
    License as published by the Free Software Foundation; either
22
    version 2.1 of the License, or (at your option) any later version.
23
 
24
    This library is distributed in the hope that it will be useful,
25
    but WITHOUT ANY WARRANTY; without even the implied warranty of
26
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
27
    Lesser General Public License for more details.
28
 
29
    You should have received a copy of the GNU Lesser General Public
30
    License along with this library; if not, write to the Free Software
31
    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307  USA */
32
 
33
/* double erf(double x)
34
 * double erfc(double x)
35
 *                           x
36
 *                    2      |\
37
 *     erf(x)  =  ---------  | exp(-t*t)dt
38
 *                 sqrt(pi) \|
39
 *                           0
40
 *
41
 *     erfc(x) =  1-erf(x)
42
 *  Note that
43
 *              erf(-x) = -erf(x)
44
 *              erfc(-x) = 2 - erfc(x)
45
 *
46
 * Method:
47
 *      1.  erf(x)  = x + x*R(x^2) for |x| in [0, 7/8]
48
 *         Remark. The formula is derived by noting
49
 *          erf(x) = (2/sqrt(pi))*(x - x^3/3 + x^5/10 - x^7/42 + ....)
50
 *         and that
51
 *          2/sqrt(pi) = 1.128379167095512573896158903121545171688
52
 *         is close to one.
53
 *
54
 *      1a. erf(x)  = 1 - erfc(x), for |x| > 1.0
55
 *          erfc(x) = 1 - erf(x)  if |x| < 1/4
56
 *
57
 *      2. For |x| in [7/8, 1], let s = |x| - 1, and
58
 *         c = 0.84506291151 rounded to single (24 bits)
59
 *      erf(s + c)  = sign(x) * (c  + P1(s)/Q1(s))
60
 *         Remark: here we use the taylor series expansion at x=1.
61
 *              erf(1+s) = erf(1) + s*Poly(s)
62
 *                       = 0.845.. + P1(s)/Q1(s)
63
 *         Note that |P1/Q1|< 0.078 for x in [0.84375,1.25]
64
 *
65
 *      3. For x in [1/4, 5/4],
66
 *      erfc(s + const)  = erfc(const)  + s P1(s)/Q1(s)
67
 *              for const = 1/4, 3/8, ..., 9/8
68
 *              and 0 <= s <= 1/8 .
69
 *
70
 *      4. For x in [5/4, 107],
71
 *      erfc(x) = (1/x)*exp(-x*x-0.5625 + R(z))
72
 *              z=1/x^2
73
 *         The interval is partitioned into several segments
74
 *         of width 1/8 in 1/x.
75
 *
76
 *      Note1:
77
 *         To compute exp(-x*x-0.5625+R/S), let s be a single
78
 *         precision number and s := x; then
79
 *              -x*x = -s*s + (s-x)*(s+x)
80
 *              exp(-x*x-0.5626+R/S) =
81
 *                      exp(-s*s-0.5625)*exp((s-x)*(s+x)+R/S);
82
 *      Note2:
83
 *         Here 4 and 5 make use of the asymptotic series
84
 *                        exp(-x*x)
85
 *              erfc(x) ~ ---------- * ( 1 + Poly(1/x^2) )
86
 *                        x*sqrt(pi)
87
 *
88
 *      5. For inf > x >= 107
89
 *      erf(x)  = sign(x) *(1 - tiny)  (raise inexact)
90
 *      erfc(x) = tiny*tiny (raise underflow) if x > 0
91
 *                      = 2 - tiny if x<0
92
 *
93
 *      7. Special case:
94
 *      erf(0)  = 0, erf(inf)  = 1, erf(-inf) = -1,
95
 *      erfc(0) = 1, erfc(inf) = 0, erfc(-inf) = 2,
96
 *              erfc/erf(NaN) is NaN
97
 */
98
 
99
#include "quadmath-imp.h"
100
 
101
 
102
 
103
__float128 erfcq (__float128);
104
 
105
 
106
/* Evaluate P[n] x^n  +  P[n-1] x^(n-1)  +  ...  +  P[0] */
107
 
108
static __float128
109
neval (__float128 x, const __float128 *p, int n)
110
{
111
  __float128 y;
112
 
113
  p += n;
114
  y = *p--;
115
  do
116
    {
117
      y = y * x + *p--;
118
    }
119
  while (--n > 0);
120
  return y;
121
}
122
 
123
 
124
/* Evaluate x^n+1  +  P[n] x^(n)  +  P[n-1] x^(n-1)  +  ...  +  P[0] */
125
 
126
static __float128
127
deval (__float128 x, const __float128 *p, int n)
128
{
129
  __float128 y;
130
 
131
  p += n;
132
  y = x + *p--;
133
  do
134
    {
135
      y = y * x + *p--;
136
    }
137
  while (--n > 0);
138
  return y;
139
}
140
 
141
 
142
 
143
static const __float128
144
tiny = 1e-4931Q,
145
  half = 0.5Q,
146
  one = 1.0Q,
147
  two = 2.0Q,
148
  /* 2/sqrt(pi) - 1 */
149
  efx = 1.2837916709551257389615890312154517168810E-1Q,
150
  /* 8 * (2/sqrt(pi) - 1) */
151
  efx8 = 1.0270333367641005911692712249723613735048E0Q;
152
 
153
 
154
/* erf(x)  = x  + x R(x^2)
155
 
156
   Peak relative error 1.8e-35  */
157
#define NTN1 8
158
static const __float128 TN1[NTN1 + 1] =
159
{
160
 -3.858252324254637124543172907442106422373E10Q,
161
  9.580319248590464682316366876952214879858E10Q,
162
  1.302170519734879977595901236693040544854E10Q,
163
  2.922956950426397417800321486727032845006E9Q,
164
  1.764317520783319397868923218385468729799E8Q,
165
  1.573436014601118630105796794840834145120E7Q,
166
  4.028077380105721388745632295157816229289E5Q,
167
  1.644056806467289066852135096352853491530E4Q,
168
  3.390868480059991640235675479463287886081E1Q
169
};
170
#define NTD1 8
171
static const __float128 TD1[NTD1 + 1] =
172
{
173
  -3.005357030696532927149885530689529032152E11Q,
174
  -1.342602283126282827411658673839982164042E11Q,
175
  -2.777153893355340961288511024443668743399E10Q,
176
  -3.483826391033531996955620074072768276974E9Q,
177
  -2.906321047071299585682722511260895227921E8Q,
178
  -1.653347985722154162439387878512427542691E7Q,
179
  -6.245520581562848778466500301865173123136E5Q,
180
  -1.402124304177498828590239373389110545142E4Q,
181
  -1.209368072473510674493129989468348633579E2Q
182
/* 1.0E0 */
183
};
184
 
185
 
186
/* erf(z+1)  = erf_const + P(z)/Q(z)
187
   -.125 <= z <= 0
188
   Peak relative error 7.3e-36  */
189
static const __float128 erf_const = 0.845062911510467529296875Q;
190
#define NTN2 8
191
static const __float128 TN2[NTN2 + 1] =
192
{
193
 -4.088889697077485301010486931817357000235E1Q,
194
  7.157046430681808553842307502826960051036E3Q,
195
 -2.191561912574409865550015485451373731780E3Q,
196
  2.180174916555316874988981177654057337219E3Q,
197
  2.848578658049670668231333682379720943455E2Q,
198
  1.630362490952512836762810462174798925274E2Q,
199
  6.317712353961866974143739396865293596895E0Q,
200
  2.450441034183492434655586496522857578066E1Q,
201
  5.127662277706787664956025545897050896203E-1Q
202
};
203
#define NTD2 8
204
static const __float128 TD2[NTD2 + 1] =
205
{
206
  1.731026445926834008273768924015161048885E4Q,
207
  1.209682239007990370796112604286048173750E4Q,
208
  1.160950290217993641320602282462976163857E4Q,
209
  5.394294645127126577825507169061355698157E3Q,
210
  2.791239340533632669442158497532521776093E3Q,
211
  8.989365571337319032943005387378993827684E2Q,
212
  2.974016493766349409725385710897298069677E2Q,
213
  6.148192754590376378740261072533527271947E1Q,
214
  1.178502892490738445655468927408440847480E1Q
215
 /* 1.0E0 */
216
};
217
 
218
 
219
/* erfc(x + 0.25) = erfc(0.25) + x R(x)
220
 
221
   Peak relative error 1.4e-35  */
222
#define NRNr13 8
223
static const __float128 RNr13[NRNr13 + 1] =
224
{
225
 -2.353707097641280550282633036456457014829E3Q,
226
  3.871159656228743599994116143079870279866E2Q,
227
 -3.888105134258266192210485617504098426679E2Q,
228
 -2.129998539120061668038806696199343094971E1Q,
229
 -8.125462263594034672468446317145384108734E1Q,
230
  8.151549093983505810118308635926270319660E0Q,
231
 -5.033362032729207310462422357772568553670E0Q,
232
 -4.253956621135136090295893547735851168471E-2Q,
233
 -8.098602878463854789780108161581050357814E-2Q
234
};
235
#define NRDr13 7
236
static const __float128 RDr13[NRDr13 + 1] =
237
{
238
  2.220448796306693503549505450626652881752E3Q,
239
  1.899133258779578688791041599040951431383E2Q,
240
  1.061906712284961110196427571557149268454E3Q,
241
  7.497086072306967965180978101974566760042E1Q,
242
  2.146796115662672795876463568170441327274E2Q,
243
  1.120156008362573736664338015952284925592E1Q,
244
  2.211014952075052616409845051695042741074E1Q,
245
  6.469655675326150785692908453094054988938E-1Q
246
 /* 1.0E0 */
247
};
248
/* erfc(0.25) = C13a + C13b to extra precision.  */
249
static const __float128 C13a = 0.723663330078125Q;
250
static const __float128 C13b = 1.0279753638067014931732235184287934646022E-5Q;
251
 
252
 
253
/* erfc(x + 0.375) = erfc(0.375) + x R(x)
254
 
255
   Peak relative error 1.2e-35  */
256
#define NRNr14 8
257
static const __float128 RNr14[NRNr14 + 1] =
258
{
259
 -2.446164016404426277577283038988918202456E3Q,
260
  6.718753324496563913392217011618096698140E2Q,
261
 -4.581631138049836157425391886957389240794E2Q,
262
 -2.382844088987092233033215402335026078208E1Q,
263
 -7.119237852400600507927038680970936336458E1Q,
264
  1.313609646108420136332418282286454287146E1Q,
265
 -6.188608702082264389155862490056401365834E0Q,
266
 -2.787116601106678287277373011101132659279E-2Q,
267
 -2.230395570574153963203348263549700967918E-2Q
268
};
269
#define NRDr14 7
270
static const __float128 RDr14[NRDr14 + 1] =
271
{
272
  2.495187439241869732696223349840963702875E3Q,
273
  2.503549449872925580011284635695738412162E2Q,
274
  1.159033560988895481698051531263861842461E3Q,
275
  9.493751466542304491261487998684383688622E1Q,
276
  2.276214929562354328261422263078480321204E2Q,
277
  1.367697521219069280358984081407807931847E1Q,
278
  2.276988395995528495055594829206582732682E1Q,
279
  7.647745753648996559837591812375456641163E-1Q
280
 /* 1.0E0 */
281
};
282
/* erfc(0.375) = C14a + C14b to extra precision.  */
283
static const __float128 C14a = 0.5958709716796875Q;
284
static const __float128 C14b = 1.2118885490201676174914080878232469565953E-5Q;
285
 
286
/* erfc(x + 0.5) = erfc(0.5) + x R(x)
287
 
288
   Peak relative error 4.7e-36  */
289
#define NRNr15 8
290
static const __float128 RNr15[NRNr15 + 1] =
291
{
292
 -2.624212418011181487924855581955853461925E3Q,
293
  8.473828904647825181073831556439301342756E2Q,
294
 -5.286207458628380765099405359607331669027E2Q,
295
 -3.895781234155315729088407259045269652318E1Q,
296
 -6.200857908065163618041240848728398496256E1Q,
297
  1.469324610346924001393137895116129204737E1Q,
298
 -6.961356525370658572800674953305625578903E0Q,
299
  5.145724386641163809595512876629030548495E-3Q,
300
  1.990253655948179713415957791776180406812E-2Q
301
};
302
#define NRDr15 7
303
static const __float128 RDr15[NRDr15 + 1] =
304
{
305
  2.986190760847974943034021764693341524962E3Q,
306
  5.288262758961073066335410218650047725985E2Q,
307
  1.363649178071006978355113026427856008978E3Q,
308
  1.921707975649915894241864988942255320833E2Q,
309
  2.588651100651029023069013885900085533226E2Q,
310
  2.628752920321455606558942309396855629459E1Q,
311
  2.455649035885114308978333741080991380610E1Q,
312
  1.378826653595128464383127836412100939126E0Q
313
  /* 1.0E0 */
314
};
315
/* erfc(0.5) = C15a + C15b to extra precision.  */
316
static const __float128 C15a = 0.4794921875Q;
317
static const __float128 C15b = 7.9346869534623172533461080354712635484242E-6Q;
318
 
319
/* erfc(x + 0.625) = erfc(0.625) + x R(x)
320
 
321
   Peak relative error 5.1e-36  */
322
#define NRNr16 8
323
static const __float128 RNr16[NRNr16 + 1] =
324
{
325
 -2.347887943200680563784690094002722906820E3Q,
326
  8.008590660692105004780722726421020136482E2Q,
327
 -5.257363310384119728760181252132311447963E2Q,
328
 -4.471737717857801230450290232600243795637E1Q,
329
 -4.849540386452573306708795324759300320304E1Q,
330
  1.140885264677134679275986782978655952843E1Q,
331
 -6.731591085460269447926746876983786152300E0Q,
332
  1.370831653033047440345050025876085121231E-1Q,
333
  2.022958279982138755020825717073966576670E-2Q,
334
};
335
#define NRDr16 7
336
static const __float128 RDr16[NRDr16 + 1] =
337
{
338
  3.075166170024837215399323264868308087281E3Q,
339
  8.730468942160798031608053127270430036627E2Q,
340
  1.458472799166340479742581949088453244767E3Q,
341
  3.230423687568019709453130785873540386217E2Q,
342
  2.804009872719893612081109617983169474655E2Q,
343
  4.465334221323222943418085830026979293091E1Q,
344
  2.612723259683205928103787842214809134746E1Q,
345
  2.341526751185244109722204018543276124997E0Q,
346
  /* 1.0E0 */
347
};
348
/* erfc(0.625) = C16a + C16b to extra precision.  */
349
static const __float128 C16a = 0.3767547607421875Q;
350
static const __float128 C16b = 4.3570693945275513594941232097252997287766E-6Q;
351
 
352
/* erfc(x + 0.75) = erfc(0.75) + x R(x)
353
 
354
   Peak relative error 1.7e-35  */
355
#define NRNr17 8
356
static const __float128 RNr17[NRNr17 + 1] =
357
{
358
  -1.767068734220277728233364375724380366826E3Q,
359
  6.693746645665242832426891888805363898707E2Q,
360
  -4.746224241837275958126060307406616817753E2Q,
361
  -2.274160637728782675145666064841883803196E1Q,
362
  -3.541232266140939050094370552538987982637E1Q,
363
  6.988950514747052676394491563585179503865E0Q,
364
  -5.807687216836540830881352383529281215100E0Q,
365
  3.631915988567346438830283503729569443642E-1Q,
366
  -1.488945487149634820537348176770282391202E-2Q
367
};
368
#define NRDr17 7
369
static const __float128 RDr17[NRDr17 + 1] =
370
{
371
  2.748457523498150741964464942246913394647E3Q,
372
  1.020213390713477686776037331757871252652E3Q,
373
  1.388857635935432621972601695296561952738E3Q,
374
  3.903363681143817750895999579637315491087E2Q,
375
  2.784568344378139499217928969529219886578E2Q,
376
  5.555800830216764702779238020065345401144E1Q,
377
  2.646215470959050279430447295801291168941E1Q,
378
  2.984905282103517497081766758550112011265E0Q,
379
  /* 1.0E0 */
380
};
381
/* erfc(0.75) = C17a + C17b to extra precision.  */
382
static const __float128 C17a = 0.2888336181640625Q;
383
static const __float128 C17b = 1.0748182422368401062165408589222625794046E-5Q;
384
 
385
 
386
/* erfc(x + 0.875) = erfc(0.875) + x R(x)
387
 
388
   Peak relative error 2.2e-35  */
389
#define NRNr18 8
390
static const __float128 RNr18[NRNr18 + 1] =
391
{
392
 -1.342044899087593397419622771847219619588E3Q,
393
  6.127221294229172997509252330961641850598E2Q,
394
 -4.519821356522291185621206350470820610727E2Q,
395
  1.223275177825128732497510264197915160235E1Q,
396
 -2.730789571382971355625020710543532867692E1Q,
397
  4.045181204921538886880171727755445395862E0Q,
398
 -4.925146477876592723401384464691452700539E0Q,
399
  5.933878036611279244654299924101068088582E-1Q,
400
 -5.557645435858916025452563379795159124753E-2Q
401
};
402
#define NRDr18 7
403
static const __float128 RDr18[NRDr18 + 1] =
404
{
405
  2.557518000661700588758505116291983092951E3Q,
406
  1.070171433382888994954602511991940418588E3Q,
407
  1.344842834423493081054489613250688918709E3Q,
408
  4.161144478449381901208660598266288188426E2Q,
409
  2.763670252219855198052378138756906980422E2Q,
410
  5.998153487868943708236273854747564557632E1Q,
411
  2.657695108438628847733050476209037025318E1Q,
412
  3.252140524394421868923289114410336976512E0Q,
413
  /* 1.0E0 */
414
};
415
/* erfc(0.875) = C18a + C18b to extra precision.  */
416
static const __float128 C18a = 0.215911865234375Q;
417
static const __float128 C18b = 1.3073705765341685464282101150637224028267E-5Q;
418
 
419
/* erfc(x + 1.0) = erfc(1.0) + x R(x)
420
 
421
   Peak relative error 1.6e-35  */
422
#define NRNr19 8
423
static const __float128 RNr19[NRNr19 + 1] =
424
{
425
 -1.139180936454157193495882956565663294826E3Q,
426
  6.134903129086899737514712477207945973616E2Q,
427
 -4.628909024715329562325555164720732868263E2Q,
428
  4.165702387210732352564932347500364010833E1Q,
429
 -2.286979913515229747204101330405771801610E1Q,
430
  1.870695256449872743066783202326943667722E0Q,
431
 -4.177486601273105752879868187237000032364E0Q,
432
  7.533980372789646140112424811291782526263E-1Q,
433
 -8.629945436917752003058064731308767664446E-2Q
434
};
435
#define NRDr19 7
436
static const __float128 RDr19[NRDr19 + 1] =
437
{
438
  2.744303447981132701432716278363418643778E3Q,
439
  1.266396359526187065222528050591302171471E3Q,
440
  1.466739461422073351497972255511919814273E3Q,
441
  4.868710570759693955597496520298058147162E2Q,
442
  2.993694301559756046478189634131722579643E2Q,
443
  6.868976819510254139741559102693828237440E1Q,
444
  2.801505816247677193480190483913753613630E1Q,
445
  3.604439909194350263552750347742663954481E0Q,
446
  /* 1.0E0 */
447
};
448
/* erfc(1.0) = C19a + C19b to extra precision.  */
449
static const __float128 C19a = 0.15728759765625Q;
450
static const __float128 C19b = 1.1609394035130658779364917390740703933002E-5Q;
451
 
452
/* erfc(x + 1.125) = erfc(1.125) + x R(x)
453
 
454
   Peak relative error 3.6e-36  */
455
#define NRNr20 8
456
static const __float128 RNr20[NRNr20 + 1] =
457
{
458
 -9.652706916457973956366721379612508047640E2Q,
459
  5.577066396050932776683469951773643880634E2Q,
460
 -4.406335508848496713572223098693575485978E2Q,
461
  5.202893466490242733570232680736966655434E1Q,
462
 -1.931311847665757913322495948705563937159E1Q,
463
 -9.364318268748287664267341457164918090611E-2Q,
464
 -3.306390351286352764891355375882586201069E0Q,
465
  7.573806045289044647727613003096916516475E-1Q,
466
 -9.611744011489092894027478899545635991213E-2Q
467
};
468
#define NRDr20 7
469
static const __float128 RDr20[NRDr20 + 1] =
470
{
471
  3.032829629520142564106649167182428189014E3Q,
472
  1.659648470721967719961167083684972196891E3Q,
473
  1.703545128657284619402511356932569292535E3Q,
474
  6.393465677731598872500200253155257708763E2Q,
475
  3.489131397281030947405287112726059221934E2Q,
476
  8.848641738570783406484348434387611713070E1Q,
477
  3.132269062552392974833215844236160958502E1Q,
478
  4.430131663290563523933419966185230513168E0Q
479
 /* 1.0E0 */
480
};
481
/* erfc(1.125) = C20a + C20b to extra precision.  */
482
static const __float128 C20a = 0.111602783203125Q;
483
static const __float128 C20b = 8.9850951672359304215530728365232161564636E-6Q;
484
 
485
/* erfc(1/x) = 1/x exp (-1/x^2 - 0.5625 + R(1/x^2))
486
   7/8 <= 1/x < 1
487
   Peak relative error 1.4e-35  */
488
#define NRNr8 9
489
static const __float128 RNr8[NRNr8 + 1] =
490
{
491
  3.587451489255356250759834295199296936784E1Q,
492
  5.406249749087340431871378009874875889602E2Q,
493
  2.931301290625250886238822286506381194157E3Q,
494
  7.359254185241795584113047248898753470923E3Q,
495
  9.201031849810636104112101947312492532314E3Q,
496
  5.749697096193191467751650366613289284777E3Q,
497
  1.710415234419860825710780802678697889231E3Q,
498
  2.150753982543378580859546706243022719599E2Q,
499
  8.740953582272147335100537849981160931197E0Q,
500
  4.876422978828717219629814794707963640913E-2Q
501
};
502
#define NRDr8 8
503
static const __float128 RDr8[NRDr8 + 1] =
504
{
505
  6.358593134096908350929496535931630140282E1Q,
506
  9.900253816552450073757174323424051765523E2Q,
507
  5.642928777856801020545245437089490805186E3Q,
508
  1.524195375199570868195152698617273739609E4Q,
509
  2.113829644500006749947332935305800887345E4Q,
510
  1.526438562626465706267943737310282977138E4Q,
511
  5.561370922149241457131421914140039411782E3Q,
512
  9.394035530179705051609070428036834496942E2Q,
513
  6.147019596150394577984175188032707343615E1Q
514
  /* 1.0E0 */
515
};
516
 
517
/* erfc(1/x) = 1/x exp (-1/x^2 - 0.5625 + R(1/x^2))
518
   0.75 <= 1/x <= 0.875
519
   Peak relative error 2.0e-36  */
520
#define NRNr7 9
521
static const __float128 RNr7[NRNr7 + 1] =
522
{
523
 1.686222193385987690785945787708644476545E1Q,
524
 1.178224543567604215602418571310612066594E3Q,
525
 1.764550584290149466653899886088166091093E4Q,
526
 1.073758321890334822002849369898232811561E5Q,
527
 3.132840749205943137619839114451290324371E5Q,
528
 4.607864939974100224615527007793867585915E5Q,
529
 3.389781820105852303125270837910972384510E5Q,
530
 1.174042187110565202875011358512564753399E5Q,
531
 1.660013606011167144046604892622504338313E4Q,
532
 6.700393957480661937695573729183733234400E2Q
533
};
534
#define NRDr7 9
535
static const __float128 RDr7[NRDr7 + 1] =
536
{
537
-1.709305024718358874701575813642933561169E3Q,
538
-3.280033887481333199580464617020514788369E4Q,
539
-2.345284228022521885093072363418750835214E5Q,
540
-8.086758123097763971926711729242327554917E5Q,
541
-1.456900414510108718402423999575992450138E6Q,
542
-1.391654264881255068392389037292702041855E6Q,
543
-6.842360801869939983674527468509852583855E5Q,
544
-1.597430214446573566179675395199807533371E5Q,
545
-1.488876130609876681421645314851760773480E4Q,
546
-3.511762950935060301403599443436465645703E2Q
547
 /* 1.0E0 */
548
};
549
 
550
/* erfc(1/x) = 1/x exp(-1/x^2 - 0.5625 + R(1/x^2))
551
   5/8 <= 1/x < 3/4
552
   Peak relative error 1.9e-35  */
553
#define NRNr6 9
554
static const __float128 RNr6[NRNr6 + 1] =
555
{
556
 1.642076876176834390623842732352935761108E0Q,
557
 1.207150003611117689000664385596211076662E2Q,
558
 2.119260779316389904742873816462800103939E3Q,
559
 1.562942227734663441801452930916044224174E4Q,
560
 5.656779189549710079988084081145693580479E4Q,
561
 1.052166241021481691922831746350942786299E5Q,
562
 9.949798524786000595621602790068349165758E4Q,
563
 4.491790734080265043407035220188849562856E4Q,
564
 8.377074098301530326270432059434791287601E3Q,
565
 4.506934806567986810091824791963991057083E2Q
566
};
567
#define NRDr6 9
568
static const __float128 RDr6[NRDr6 + 1] =
569
{
570
-1.664557643928263091879301304019826629067E2Q,
571
-3.800035902507656624590531122291160668452E3Q,
572
-3.277028191591734928360050685359277076056E4Q,
573
-1.381359471502885446400589109566587443987E5Q,
574
-3.082204287382581873532528989283748656546E5Q,
575
-3.691071488256738343008271448234631037095E5Q,
576
-2.300482443038349815750714219117566715043E5Q,
577
-6.873955300927636236692803579555752171530E4Q,
578
-8.262158817978334142081581542749986845399E3Q,
579
-2.517122254384430859629423488157361983661E2Q
580
 /* 1.00 */
581
};
582
 
583
/* erfc(1/x) = 1/x exp(-1/x^2 - 0.5625 + R(1/x^2))
584
   1/2 <= 1/x < 5/8
585
   Peak relative error 4.6e-36  */
586
#define NRNr5 10
587
static const __float128 RNr5[NRNr5 + 1] =
588
{
589
-3.332258927455285458355550878136506961608E-3Q,
590
-2.697100758900280402659586595884478660721E-1Q,
591
-6.083328551139621521416618424949137195536E0Q,
592
-6.119863528983308012970821226810162441263E1Q,
593
-3.176535282475593173248810678636522589861E2Q,
594
-8.933395175080560925809992467187963260693E2Q,
595
-1.360019508488475978060917477620199499560E3Q,
596
-1.075075579828188621541398761300910213280E3Q,
597
-4.017346561586014822824459436695197089916E2Q,
598
-5.857581368145266249509589726077645791341E1Q,
599
-2.077715925587834606379119585995758954399E0Q
600
};
601
#define NRDr5 9
602
static const __float128 RDr5[NRDr5 + 1] =
603
{
604
 3.377879570417399341550710467744693125385E-1Q,
605
 1.021963322742390735430008860602594456187E1Q,
606
 1.200847646592942095192766255154827011939E2Q,
607
 7.118915528142927104078182863387116942836E2Q,
608
 2.318159380062066469386544552429625026238E3Q,
609
 4.238729853534009221025582008928765281620E3Q,
610
 4.279114907284825886266493994833515580782E3Q,
611
 2.257277186663261531053293222591851737504E3Q,
612
 5.570475501285054293371908382916063822957E2Q,
613
 5.142189243856288981145786492585432443560E1Q
614
 /* 1.0E0 */
615
};
616
 
617
/* erfc(1/x) = 1/x exp(-1/x^2 - 0.5625 + R(1/x^2))
618
   3/8 <= 1/x < 1/2
619
   Peak relative error 2.0e-36  */
620
#define NRNr4 10
621
static const __float128 RNr4[NRNr4 + 1] =
622
{
623
 3.258530712024527835089319075288494524465E-3Q,
624
 2.987056016877277929720231688689431056567E-1Q,
625
 8.738729089340199750734409156830371528862E0Q,
626
 1.207211160148647782396337792426311125923E2Q,
627
 8.997558632489032902250523945248208224445E2Q,
628
 3.798025197699757225978410230530640879762E3Q,
629
 9.113203668683080975637043118209210146846E3Q,
630
 1.203285891339933238608683715194034900149E4Q,
631
 8.100647057919140328536743641735339740855E3Q,
632
 2.383888249907144945837976899822927411769E3Q,
633
 2.127493573166454249221983582495245662319E2Q
634
};
635
#define NRDr4 10
636
static const __float128 RDr4[NRDr4 + 1] =
637
{
638
-3.303141981514540274165450687270180479586E-1Q,
639
-1.353768629363605300707949368917687066724E1Q,
640
-2.206127630303621521950193783894598987033E2Q,
641
-1.861800338758066696514480386180875607204E3Q,
642
-8.889048775872605708249140016201753255599E3Q,
643
-2.465888106627948210478692168261494857089E4Q,
644
-3.934642211710774494879042116768390014289E4Q,
645
-3.455077258242252974937480623730228841003E4Q,
646
-1.524083977439690284820586063729912653196E4Q,
647
-2.810541887397984804237552337349093953857E3Q,
648
-1.343929553541159933824901621702567066156E2Q
649
 /* 1.0E0 */
650
};
651
 
652
/* erfc(1/x) = 1/x exp(-1/x^2 - 0.5625 + R(1/x^2))
653
   1/4 <= 1/x < 3/8
654
   Peak relative error 8.4e-37  */
655
#define NRNr3 11
656
static const __float128 RNr3[NRNr3 + 1] =
657
{
658
-1.952401126551202208698629992497306292987E-6Q,
659
-2.130881743066372952515162564941682716125E-4Q,
660
-8.376493958090190943737529486107282224387E-3Q,
661
-1.650592646560987700661598877522831234791E-1Q,
662
-1.839290818933317338111364667708678163199E0Q,
663
-1.216278715570882422410442318517814388470E1Q,
664
-4.818759344462360427612133632533779091386E1Q,
665
-1.120994661297476876804405329172164436784E2Q,
666
-1.452850765662319264191141091859300126931E2Q,
667
-9.485207851128957108648038238656777241333E1Q,
668
-2.563663855025796641216191848818620020073E1Q,
669
-1.787995944187565676837847610706317833247E0Q
670
};
671
#define NRDr3 10
672
static const __float128 RDr3[NRDr3 + 1] =
673
{
674
 1.979130686770349481460559711878399476903E-4Q,
675
 1.156941716128488266238105813374635099057E-2Q,
676
 2.752657634309886336431266395637285974292E-1Q,
677
 3.482245457248318787349778336603569327521E0Q,
678
 2.569347069372696358578399521203959253162E1Q,
679
 1.142279000180457419740314694631879921561E2Q,
680
 3.056503977190564294341422623108332700840E2Q,
681
 4.780844020923794821656358157128719184422E2Q,
682
 4.105972727212554277496256802312730410518E2Q,
683
 1.724072188063746970865027817017067646246E2Q,
684
 2.815939183464818198705278118326590370435E1Q
685
 /* 1.0E0 */
686
};
687
 
688
/* erfc(1/x) = 1/x exp(-1/x^2 - 0.5625 + R(1/x^2))
689
   1/8 <= 1/x < 1/4
690
   Peak relative error 1.5e-36  */
691
#define NRNr2 11
692
static const __float128 RNr2[NRNr2 + 1] =
693
{
694
-2.638914383420287212401687401284326363787E-8Q,
695
-3.479198370260633977258201271399116766619E-6Q,
696
-1.783985295335697686382487087502222519983E-4Q,
697
-4.777876933122576014266349277217559356276E-3Q,
698
-7.450634738987325004070761301045014986520E-2Q,
699
-7.068318854874733315971973707247467326619E-1Q,
700
-4.113919921935944795764071670806867038732E0Q,
701
-1.440447573226906222417767283691888875082E1Q,
702
-2.883484031530718428417168042141288943905E1Q,
703
-2.990886974328476387277797361464279931446E1Q,
704
-1.325283914915104866248279787536128997331E1Q,
705
-1.572436106228070195510230310658206154374E0Q
706
};
707
#define NRDr2 10
708
static const __float128 RDr2[NRDr2 + 1] =
709
{
710
 2.675042728136731923554119302571867799673E-6Q,
711
 2.170997868451812708585443282998329996268E-4Q,
712
 7.249969752687540289422684951196241427445E-3Q,
713
 1.302040375859768674620410563307838448508E-1Q,
714
 1.380202483082910888897654537144485285549E0Q,
715
 8.926594113174165352623847870299170069350E0Q,
716
 3.521089584782616472372909095331572607185E1Q,
717
 8.233547427533181375185259050330809105570E1Q,
718
 1.072971579885803033079469639073292840135E2Q,
719
 6.943803113337964469736022094105143158033E1Q,
720
 1.775695341031607738233608307835017282662E1Q
721
 /* 1.0E0 */
722
};
723
 
724
/* erfc(1/x) = 1/x exp(-1/x^2 - 0.5625 + R(1/x^2))
725
   1/128 <= 1/x < 1/8
726
   Peak relative error 2.2e-36  */
727
#define NRNr1 9
728
static const __float128 RNr1[NRNr1 + 1] =
729
{
730
-4.250780883202361946697751475473042685782E-8Q,
731
-5.375777053288612282487696975623206383019E-6Q,
732
-2.573645949220896816208565944117382460452E-4Q,
733
-6.199032928113542080263152610799113086319E-3Q,
734
-8.262721198693404060380104048479916247786E-2Q,
735
-6.242615227257324746371284637695778043982E-1Q,
736
-2.609874739199595400225113299437099626386E0Q,
737
-5.581967563336676737146358534602770006970E0Q,
738
-5.124398923356022609707490956634280573882E0Q,
739
-1.290865243944292370661544030414667556649E0Q
740
};
741
#define NRDr1 8
742
static const __float128 RDr1[NRDr1 + 1] =
743
{
744
 4.308976661749509034845251315983612976224E-6Q,
745
 3.265390126432780184125233455960049294580E-4Q,
746
 9.811328839187040701901866531796570418691E-3Q,
747
 1.511222515036021033410078631914783519649E-1Q,
748
 1.289264341917429958858379585970225092274E0Q,
749
 6.147640356182230769548007536914983522270E0Q,
750
 1.573966871337739784518246317003956180750E1Q,
751
 1.955534123435095067199574045529218238263E1Q,
752
 9.472613121363135472247929109615785855865E0Q
753
  /* 1.0E0 */
754
};
755
 
756
 
757
__float128
758
erfq (__float128 x)
759
{
760
  __float128 a, y, z;
761
  int32_t i, ix, sign;
762
  ieee854_float128 u;
763
 
764
  u.value = x;
765
  sign = u.words32.w0;
766
  ix = sign & 0x7fffffff;
767
 
768
  if (ix >= 0x7fff0000)
769
    {                           /* erf(nan)=nan */
770
      i = ((sign & 0xffff0000) >> 31) << 1;
771
      return (__float128) (1 - i) + one / x;    /* erf(+-inf)=+-1 */
772
    }
773
 
774
  if (ix >= 0x3fff0000) /* |x| >= 1.0 */
775
    {
776
      y = erfcq (x);
777
      return (one - y);
778
      /*    return (one - erfcq (x)); */
779
    }
780
  u.words32.w0 = ix;
781
  a = u.value;
782
  z = x * x;
783
  if (ix < 0x3ffec000)  /* a < 0.875 */
784
    {
785
      if (ix < 0x3fc60000) /* |x|<2**-57 */
786
        {
787
          if (ix < 0x00080000)
788
            return 0.125 * (8.0 * x + efx8 * x);        /*avoid underflow */
789
          return x + efx * x;
790
        }
791
      y = a + a * neval (z, TN1, NTN1) / deval (z, TD1, NTD1);
792
    }
793
  else
794
    {
795
      a = a - one;
796
      y = erf_const + neval (a, TN2, NTN2) / deval (a, TD2, NTD2);
797
    }
798
 
799
  if (sign & 0x80000000) /* x < 0 */
800
    y = -y;
801
  return( y );
802
}
803
 
804
 
805
__float128
806
erfcq (__float128 x)
807
{
808
  __float128 y = 0.0Q, z, p, r;
809
  int32_t i, ix, sign;
810
  ieee854_float128 u;
811
 
812
  u.value = x;
813
  sign = u.words32.w0;
814
  ix = sign & 0x7fffffff;
815
  u.words32.w0 = ix;
816
 
817
  if (ix >= 0x7fff0000)
818
    {                           /* erfc(nan)=nan */
819
      /* erfc(+-inf)=0,2 */
820
      return (__float128) (((uint32_t) sign >> 31) << 1) + one / x;
821
    }
822
 
823
  if (ix < 0x3ffd0000) /* |x| <1/4 */
824
    {
825
      if (ix < 0x3f8d0000) /* |x|<2**-114 */
826
        return one - x;
827
      return one - erfq (x);
828
    }
829
  if (ix < 0x3fff4000) /* 1.25 */
830
    {
831
      x = u.value;
832
      i = 8.0 * x;
833
      switch (i)
834
        {
835
        case 2:
836
          z = x - 0.25Q;
837
          y = C13b + z * neval (z, RNr13, NRNr13) / deval (z, RDr13, NRDr13);
838
          y += C13a;
839
          break;
840
        case 3:
841
          z = x - 0.375Q;
842
          y = C14b + z * neval (z, RNr14, NRNr14) / deval (z, RDr14, NRDr14);
843
          y += C14a;
844
          break;
845
        case 4:
846
          z = x - 0.5Q;
847
          y = C15b + z * neval (z, RNr15, NRNr15) / deval (z, RDr15, NRDr15);
848
          y += C15a;
849
          break;
850
        case 5:
851
          z = x - 0.625Q;
852
          y = C16b + z * neval (z, RNr16, NRNr16) / deval (z, RDr16, NRDr16);
853
          y += C16a;
854
          break;
855
        case 6:
856
          z = x - 0.75Q;
857
          y = C17b + z * neval (z, RNr17, NRNr17) / deval (z, RDr17, NRDr17);
858
          y += C17a;
859
          break;
860
        case 7:
861
          z = x - 0.875Q;
862
          y = C18b + z * neval (z, RNr18, NRNr18) / deval (z, RDr18, NRDr18);
863
          y += C18a;
864
          break;
865
        case 8:
866
          z = x - 1.0Q;
867
          y = C19b + z * neval (z, RNr19, NRNr19) / deval (z, RDr19, NRDr19);
868
          y += C19a;
869
          break;
870
        case 9:
871
          z = x - 1.125Q;
872
          y = C20b + z * neval (z, RNr20, NRNr20) / deval (z, RDr20, NRDr20);
873
          y += C20a;
874
          break;
875
        }
876
      if (sign & 0x80000000)
877
        y = 2.0Q - y;
878
      return y;
879
    }
880
  /* 1.25 < |x| < 107 */
881
  if (ix < 0x4005ac00)
882
    {
883
      /* x < -9 */
884
      if ((ix >= 0x40022000) && (sign & 0x80000000))
885
        return two - tiny;
886
 
887
      x = fabsq (x);
888
      z = one / (x * x);
889
      i = 8.0 / x;
890
      switch (i)
891
        {
892
        default:
893
        case 0:
894
          p = neval (z, RNr1, NRNr1) / deval (z, RDr1, NRDr1);
895
          break;
896
        case 1:
897
          p = neval (z, RNr2, NRNr2) / deval (z, RDr2, NRDr2);
898
          break;
899
        case 2:
900
          p = neval (z, RNr3, NRNr3) / deval (z, RDr3, NRDr3);
901
          break;
902
        case 3:
903
          p = neval (z, RNr4, NRNr4) / deval (z, RDr4, NRDr4);
904
          break;
905
        case 4:
906
          p = neval (z, RNr5, NRNr5) / deval (z, RDr5, NRDr5);
907
          break;
908
        case 5:
909
          p = neval (z, RNr6, NRNr6) / deval (z, RDr6, NRDr6);
910
          break;
911
        case 6:
912
          p = neval (z, RNr7, NRNr7) / deval (z, RDr7, NRDr7);
913
          break;
914
        case 7:
915
          p = neval (z, RNr8, NRNr8) / deval (z, RDr8, NRDr8);
916
          break;
917
        }
918
      u.value = x;
919
      u.words32.w3 = 0;
920
      u.words32.w2 &= 0xfe000000;
921
      z = u.value;
922
      r = expq (-z * z - 0.5625) * expq ((z - x) * (z + x) + p);
923
      if ((sign & 0x80000000) == 0)
924
        return r / x;
925
      else
926
        return two - r / x;
927
    }
928
  else
929
    {
930
      if ((sign & 0x80000000) == 0)
931
        return tiny * tiny;
932
      else
933
        return two - tiny;
934
    }
935
}

powered by: WebSVN 2.1.0

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