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-2025 OpenCores.org, equivalent to Oliscience, all rights reserved. OpenCores®, registered trademark.