overby.me
this repo has no description
1/*
2 * Copyright (c) 2002 Apple Computer, Inc. All rights reserved.
3 *
4 * @APPLE_LICENSE_HEADER_START@
5 *
6 * The contents of this file constitute Original Code as defined in and
7 * are subject to the Apple Public Source License Version 1.1 (the
8 * "License"). You may not use this file except in compliance with the
9 * License. Please obtain a copy of the License at
10 * http://www.apple.com/publicsource and read it before using this file.
11 *
12 * This Original Code and all software distributed under the License are
13 * distributed on an "AS IS" basis, WITHOUT WARRANTY OF ANY KIND, EITHER
14 * EXPRESS OR IMPLIED, AND APPLE HEREBY DISCLAIMS ALL SUCH WARRANTIES,
15 * INCLUDING WITHOUT LIMITATION, ANY WARRANTIES OF MERCHANTABILITY,
16 * FITNESS FOR A PARTICULAR PURPOSE OR NON-INFRINGEMENT. Please see the
17 * License for the specific language governing rights and limitations
18 * under the License.
19 *
20 * @APPLE_LICENSE_HEADER_END@
21 */
22// ErfDD.c
23//
24// Double-double Function Library
25// Copyright: � 1995-96 by Apple Computer, Inc., all rights reserved
26//
27// long double erfl( long double x );
28// long double erfcl( long double x );
29//
30
31#include "math.h"
32#include "fp_private.h"
33#if (__WANT_LONG_DOUBLE_FORMAT__ - 0L == 128L)
34#include "DD.h"
35
36// Requires the following functions:
37
38void _d3div(double a, double b, double c, double d, double e, double f,
39 double *high, double *mid, double *low);
40void _d3mul(double a, double b, double c, double d, double e, double f,
41 double *high, double *mid, double *low);
42void _d3add(double a, double b, double c, double d, double e, double f,
43 double *high, double *mid, double *low);
44long double _ExpInnerLD(double alpha, double beta, double gamma,
45 double *extra, int exptype);
46
47// real*8 erfxtr(0:2)
48// FORTRAN: erfxtr(i) C: erfxtr[i]
49static const uint32_t erfxtrData[] __attribute__ ((aligned(8))) = {
50 0x00000000, 0x00000000,
51 0xB9155555, 0x55555555,
52 0x38F99999, 0x9999999A
53};
54
55// real*8 erfpow(1:2, 0:21)
56// FORTRAN: erfpow(i,j) C: erfpow[j][i-1]
57static const uint32_t erfpowData[] __attribute__ ((aligned(8))) = {
58 0x3FF00000, 0x00000000, 0x00000000, 0x00000000,
59 0xBFD55555, 0x55555555, 0xBC755555, 0x55555555,
60 0x3FB99999, 0x9999999A, 0xBC599999, 0x9999999A,
61 0xBF986186, 0x18618618, 0xBC386186, 0x18618613,
62 0x3F72F684, 0xBDA12F68, 0x3C12F684, 0xBDA15772,
63 0xBF48D301, 0x8D3018D3, 0xBB88D301, 0xBA2FF990,
64 0x3F1C01C0, 0x1C01C01C, 0x3B5C01FF, 0xEF4E45A0,
65 0xBEEBBD77, 0x9334EF0B, 0x3B84E5AE, 0x2FF08B10,
66 0x3EB87A00, 0x187A0018, 0x3B5EC7BC, 0xA752F683,
67 0xBE83777C, 0x55568CCD, 0xBB29453D, 0xDF9ADD3E,
68 0x3E4C2E30, 0x54870B52, 0xBADB98FC, 0x71C71C77,
69 0xBE12B673, 0x10AAA27C, 0x3AACE36B, 0x74F03296,
70 0x3DD6F448, 0xE13F19E7, 0xBA7D15ED, 0x097B425F,
71 0xBD9A289E, 0xE7F6D02A, 0xBA21F9AD, 0xD3C0CA47,
72 0x3D5BD577, 0xE7FBCC09, 0xB9E48B0F, 0xCD6E9E08,
73 0xBD1BC625, 0x25C514CA, 0xB9B161F9, 0xADD3C0CA,
74 0x3CDA173A, 0x4FC330FD, 0xB9648B0F, 0xCD6E9E03,
75 0xBC97270A, 0xC344BDA1, 0xB927B425, 0xED097B45,
76 0x3C537610, 0xC735BA78, 0x38D948B0, 0xFCD6E9DE,
77 0xBC0EF15A, 0xDD3C0CA4, 0xB8A61F9A, 0xDD3C0CA5,
78 0x3BC67194, 0x8B0FCD6F, 0xB8687E6B, 0x74F03292,
79 0xBB778E38, 0xE38E38E4, 0x381C71C7, 0x1C71C71C
80};
81
82// real*8 coeff(2,0:28,8)
83// FORTRAN: coeff(i,j,k) C: coeff[k][j][i-1]
84static const uint32_t coeffData[] __attribute__ ((aligned(8))) = {
85 0x3FEFFFFF, 0xF87B4F81, 0xBC74B6DC, 0x82E2B33B,
86 0xBFDFFFF9, 0x651B9CCD, 0xBC77993A, 0x7E74C502,
87 0x3FE7FF4C, 0x63E2B209, 0x3C607D84, 0x888EE934,
88 0xBFFDF3AD, 0x307BCFA0, 0x3C80189F, 0x21239DB8,
89 0x4019F11D, 0x8747B2EA, 0x3CB282E6, 0x775914EA,
90 0xC03BFD72, 0x2186F691, 0x3CD1B8A8, 0xA64D28F3,
91 0x406123E2, 0xF4070A2E, 0xBD0AA69F, 0x832442C8,
92 0xC085E672, 0xCAE6452C, 0xBD264E3F, 0x5B4D5E70,
93 0x40AB48F2, 0xAF9DF505, 0x3D3CFDE8, 0x35619194,
94 0xC0CF99E0, 0xAC6AFA01, 0xBD525B2E, 0x7AEBA4B8,
95 0x40F07ABD, 0x450474EE, 0x3D9F4814, 0x21020E9F,
96 0xC10E4BB4, 0xEAA37C6D, 0x3DADDD80, 0x6CA8C8F2,
97 0x41282796, 0x01A20544, 0x3DC2C51E, 0x7D2AD0DC,
98 0xC1407A5E, 0xA1893145, 0xBDEB8D59, 0x4611C6D5,
99 0x4152F9A4, 0x9A0ED425, 0xBDF53C28, 0x002DF666,
100 0xC16227D6, 0x896B5A9A, 0x3DD50BAE, 0xCA782EDC,
101 0x416C4D5B, 0x1A0268E7, 0xBDE19C5A, 0x4D09E56C,
102 0xC1717E27, 0xC2AB2A36, 0x3E122CE5, 0xA75631BC,
103 0x41707C32, 0x89EC6733, 0xBE1F4B89, 0xEF2106D9,
104 0xC1663DCE, 0x61AF76F4, 0x3E028503, 0xD6CDBB9E,
105 0x41531F97, 0x0C38EEF3, 0xBDE20467, 0x5FB2BCF0,
106 0xC12F7BE6, 0x257D1D65, 0x3D8DFD02, 0xE8C06000,
107 0x00000000, 0x00000000, 0x00000000, 0x00000000,
108 0x00000000, 0x00000000, 0x00000000, 0x00000000,
109 0x00000000, 0x00000000, 0x00000000, 0x00000000,
110 0x00000000, 0x00000000, 0x00000000, 0x00000000,
111 0x00000000, 0x00000000, 0x00000000, 0x00000000,
112 0x00000000, 0x00000000, 0x00000000, 0x00000000,
113 0x00000000, 0x00000000, 0x00000000, 0x00000000,
114
115 0x3FEFFFFF, 0xFFEBC039, 0x3C529C69, 0x6CF540A8,
116 0xBFDFFFFF, 0xE38F77CF, 0x3BFCDCBB, 0x0CD77D00,
117 0x3FE7FFFB, 0x2DCCB103, 0x3C8175D4, 0xB96461CA,
118 0xBFFDFF79, 0x49317263, 0xBC906876, 0xAC321F5A,
119 0x401A3AAB, 0x9BF9D0BD, 0x3CABAD03, 0xD57E6058,
120 0xC03D5E32, 0x130E26B3, 0xBC9E17BB, 0x1F622DC0,
121 0x4063C8DE, 0x39F025C5, 0x3CEE695F, 0x35A91AF0,
122 0xC08E3E6E, 0x953DCEB7, 0x3D24E053, 0xE29EAEBE,
123 0x40B8AA52, 0x3AD56173, 0x3D50EDEF, 0x7FD5BCF1,
124 0xC0E440E6, 0xCEADBFBB, 0x3D71A20A, 0x51AF181A,
125 0x410FEE4B, 0xDCA4E1E0, 0x3D8CCED6, 0xFF021BF0,
126 0xC13752BC, 0xAAA09717, 0xBDDB430A, 0x583D9DFE,
127 0x415EC821, 0x1EEE2921, 0x3DDA8E50, 0xAC97C8D8,
128 0xC182031E, 0xEDA7A8C5, 0x3E23FDA5, 0x816A17A9,
129 0x41A26A35, 0x43351D14, 0x3E4A1C99, 0x4FE13BE4,
130 0xC1C0397D, 0x3FF8253E, 0xBE55B81D, 0xBAD43C68,
131 0x41D84B8E, 0x7B8414AB, 0x3E797BE5, 0x69A73302,
132 0xC1EE6B6B, 0x349E71AE, 0xBE7EB248, 0x45D29614,
133 0x41FF34DA, 0xF64919BF, 0xBE87460E, 0x4926C028,
134 0xC20984E7, 0x864FAEBA, 0x3EA1DD11, 0x07234030,
135 0x420FF82D, 0xCF64D788, 0x3E9B32FC, 0x90159470,
136 0xC20CCAA6, 0x72BB40E4, 0x3E619A97, 0x18E25640,
137 0x420096A8, 0x68135449, 0xBEAFFD28, 0xED6BBBFA,
138 0xC1E25D70, 0xCAB7C34C, 0xBE5B42DB, 0xEC71D440,
139 0x00000000, 0x00000000, 0x00000000, 0x00000000,
140 0x00000000, 0x00000000, 0x00000000, 0x00000000,
141 0x00000000, 0x00000000, 0x00000000, 0x00000000,
142 0x00000000, 0x00000000, 0x00000000, 0x00000000,
143 0x00000000, 0x00000000, 0x00000000, 0x00000000,
144
145 0x3FEFFFFF, 0xFFFFE639, 0xBC76B4F8, 0xF22A1D6A,
146 0xBFDFFFFF, 0xFFC4926E, 0x3C74AA2B, 0x05602AA4,
147 0x3FE7FFFF, 0xEF8F13A8, 0xBC802572, 0x96640478,
148 0xBFFDFFFD, 0x175A90A3, 0x3C8448AC, 0x5447AC1C,
149 0x401A3FD0, 0xA3F30955, 0x3CBDC702, 0x0AF0E74C,
150 0xC03D85B0, 0xDE15683B, 0x3CD775BC, 0x8D3FFDA1,
151 0x406441FD, 0xA9731CF2, 0xBCEC6B02, 0xE82A19B0,
152 0xC0904FCE, 0xB4270577, 0x3CF73CA2, 0xBC7037A0,
153 0x40BDA90A, 0xA9B95F59, 0xBD5B5C9C, 0x5ABA6422,
154 0xC0ED1C94, 0x9666E63C, 0xBD5C789C, 0x5A7928DC,
155 0x411D56F4, 0x20DB5218, 0x3DB320CC, 0xDE074D1E,
156 0xC14CEC6A, 0x9FD76B11, 0x3DE9DE29, 0x03D42F2C,
157 0x417AC1F5, 0xF5446DF4, 0xBE031819, 0x3468BE58,
158 0xC1A67958, 0xAC7171B2, 0xBE3581A4, 0xB51D7A74,
159 0x41D0B077, 0x6E96B0B2, 0x3E574081, 0x998F5AA4,
160 0xC1F56C00, 0xFD9F8D53, 0xBE94D297, 0xFD7D3E08,
161 0x42173F0C, 0x0EB8D9D7, 0x3EBE2714, 0x71045A4A,
162 0xC234D11D, 0x4DA55FFA, 0x3EC4EA18, 0x03DADE62,
163 0x424DDAFE, 0x5B2DFE46, 0x3EE6D6D1, 0x58677012,
164 0xC2607209, 0x5B35C7C9, 0x3F04E6A4, 0x0B2429BA,
165 0x426A1766, 0x45D71E26, 0xBF0B66BC, 0xFE7F12C0,
166 0xC26A80F6, 0x5FAAB021, 0x3F042804, 0xC0D4DFB5,
167 0x4259DF2B, 0x9341D4CC, 0xBEE2816A, 0x58CBD636,
168 0x00000000, 0x00000000, 0x00000000, 0x00000000,
169 0x00000000, 0x00000000, 0x00000000, 0x00000000,
170 0x00000000, 0x00000000, 0x00000000, 0x00000000,
171 0x00000000, 0x00000000, 0x00000000, 0x00000000,
172 0x00000000, 0x00000000, 0x00000000, 0x00000000,
173 0x00000000, 0x00000000, 0x00000000, 0x00000000,
174
175 0x3FF00000, 0x00000000, 0xBB75CA29, 0xC48EA000,
176 0xBFE00000, 0x00000000, 0x3C1A999F, 0xD44FCE00,
177 0x3FE7FFFF, 0xFFFFFFFE, 0x3C642A4A, 0xB64D8068,
178 0xBFFDFFFF, 0xFFFFFE99, 0xBC9ABDD0, 0x98308E2C,
179 0x401A3FFF, 0xFFFFA343, 0x3CB85B18, 0x87B93172,
180 0xC03D87FF, 0xFFEDDBFB, 0x3CD309B5, 0xEEE35E41,
181 0x40644D7F, 0xFE9A5199, 0x3CFC5B15, 0x3519E834,
182 0xC0907EF7, 0xE9B322F4, 0xBD35FCB1, 0x58754A9E,
183 0x40BEEE0E, 0xB41A926D, 0xBD57FB38, 0x347B51CF,
184 0xC0F06E6C, 0x572C5E36, 0x3D720A1F, 0x8B0D1254,
185 0x412382B8, 0xB55E07E1, 0xBDCABBAE, 0x93440AE8,
186 0xC1599872, 0x5E40B785, 0x3DF51615, 0xACBF5F19,
187 0x41925B4B, 0x43808517, 0xBE3C17EF, 0x5CEF40F8,
188 0xC1CC74A4, 0x3D911FB7, 0x3E172705, 0xD877BCC0,
189 0x42077587, 0x54FEECDC, 0xBEA98313, 0x07804B47,
190 0xC2441B83, 0xD7BACF31, 0xBEEA3073, 0x0D0DBA05,
191 0x428165FB, 0x78789F6C, 0xBED45B8E, 0xD3B8907C,
192 0xC2BD68C3, 0xB0580084, 0x3F1BFE68, 0xE7B64DE0,
193 0x42F78115, 0xFD9904B1, 0xBF8FED1E, 0x29D74809,
194 0xC3313FA3, 0xDBF81756, 0x3FDC570B, 0xC88AFF19,
195 0x4366A5A6, 0xB262F5E4, 0x3FFA9257, 0x0A3AD3F8,
196 0xC399F8EF, 0xB5A4E015, 0x401B7ECE, 0x1B301370,
197 0x43C9687F, 0xB88D7E83, 0xC06CE8BA, 0x85792A6C,
198 0xC3F4A8E4, 0x451D5574, 0x409B5C8E, 0x70046666,
199 0x441B0DCA, 0xB2FAB836, 0x40996671, 0xF8BB2152,
200 0xC43B52DB, 0xC1E7B2EE, 0x40DB5AF6, 0x615EB7EE,
201 0x4453ECD3, 0xF252D951, 0x40FC199B, 0x4FAD14EA,
202 0xC462A264, 0xDCCB32FF, 0xC0E17D1D, 0x52615E34,
203 0x4460C115, 0xAC203924, 0x40E780D0, 0xA1F70E00,
204
205 0x3FF00000, 0x00000000, 0x3A9A41C2, 0x18000000,
206 0xBFE00000, 0x00000000, 0xBB5DCD07, 0x8C200000,
207 0x3FE80000, 0x00000000, 0x3C0EBA46, 0xBCD71000,
208 0xBFFE0000, 0x00000001, 0xBC88EE2A, 0x6486B1E4,
209 0x401A4000, 0x00000101, 0x3CB25362, 0xC332A010,
210 0xC03D8800, 0x00009B65, 0x3CD119CB, 0x15915886,
211 0x40644D80, 0x0022EC9B, 0xBCE9E9C4, 0x0897A1AC,
212 0xC0907EF8, 0x05F94723, 0x3D2AB506, 0xB989F798,
213 0x40BEEE12, 0x939C8382, 0x3D5AD230, 0x743AEFD8,
214 0xC0F06E8D, 0xB86AB92B, 0x3D98E5AE, 0xBE0C9174,
215 0x412384D6, 0x161AD32F, 0x3DBB5BF4, 0xEF6A7DBF,
216 0xC159B642, 0x7281154F, 0x3DFECD4C, 0xAA7F47EC,
217 0x419305F9, 0xC7BF95B4, 0xBE349B0E, 0x9D9581DE,
218 0xC1D12C24, 0x5D3084BD, 0xBE5E4995, 0xB96B303C,
219 0x4215451C, 0xDF186F15, 0x3EA8548E, 0x669919A8,
220 0xC25F3153, 0x01C60944, 0x3ED326BD, 0x842F8568,
221 0x42A25C17, 0xD2287EF6, 0x3F493E23, 0x133CA4B2,
222 0xC2D613F5, 0x225EEB4A, 0x3F4A86AC, 0x4B33D440,
223 0x00000000, 0x00000000, 0x00000000, 0x00000000,
224 0x00000000, 0x00000000, 0x00000000, 0x00000000,
225 0x00000000, 0x00000000, 0x00000000, 0x00000000,
226 0x00000000, 0x00000000, 0x00000000, 0x00000000,
227 0x00000000, 0x00000000, 0x00000000, 0x00000000,
228 0x00000000, 0x00000000, 0x00000000, 0x00000000,
229 0x00000000, 0x00000000, 0x00000000, 0x00000000,
230 0x00000000, 0x00000000, 0x00000000, 0x00000000,
231 0x00000000, 0x00000000, 0x00000000, 0x00000000,
232 0x00000000, 0x00000000, 0x00000000, 0x00000000,
233 0x00000000, 0x00000000, 0x00000000, 0x00000000,
234
235 0xBFA4AFAB, 0x9B6E1105, 0x3C463A99, 0xD0AF98EE,
236 0xBFD4F526, 0x85059209, 0xBC73ECA6, 0xE2948B70,
237 0xBFEA28C7, 0xA68495F3, 0xBC81EAF6, 0x808C4753,
238 0xBFA980E7, 0x79855BE9, 0x3C30F43B, 0xBDDFAE98,
239 0x3F8895D7, 0xF2B826D4, 0x3C250FDF, 0xAF15221E,
240 0xBF6222D1, 0xDE8EFFF2, 0x3BFC3D53, 0x8B04AE74,
241 0x3F26522D, 0x3A826B8E, 0xBBBFB440, 0x1723702C,
242 0x3F15289F, 0x67803BEA, 0x3BBA9594, 0x12B030A4,
243 0xBF08627F, 0xFCF37770, 0x3BAD7A54, 0xCCB615F2,
244 0x3EEBEB7F, 0xBC45B14A, 0xBB715091, 0x8DB6D7B4,
245 0xBEBE0FDE, 0x30CC0705, 0xBB3ED352, 0x5BE0D640,
246 0xBE9B7052, 0x9A2F9B59, 0xBB29B1CC, 0xF66FE0B8,
247 0x3E97B4BD, 0x41755DF3, 0xBB10EE4A, 0xCBA6B5AC,
248 0xBE807536, 0x0AA20FDE, 0xBAE9E0B1, 0xF22797C0,
249 0x3E57D20D, 0xE4C07B40, 0x3AFADF83, 0x4F8A74E3,
250 0x3E1E7ADC, 0x1A5BC732, 0xBAACBD1B, 0xDB944806,
251 0xBE29CB52, 0xFF49FDEF, 0xBACE7432, 0x1EF280D9,
252 0x3E14D367, 0x091F84DD, 0x3ABD9BF9, 0x8C93A4B0,
253 0xBDF32CB0, 0xA2D996E6, 0x3A947D3D, 0xF3852A64,
254 0x3DC1332B, 0x63A393F5, 0xBA666775, 0x70308F4A,
255 0x00000000, 0x00000000, 0x00000000, 0x00000000,
256 0x00000000, 0x00000000, 0x00000000, 0x00000000,
257 0x00000000, 0x00000000, 0x00000000, 0x00000000,
258 0x00000000, 0x00000000, 0x00000000, 0x00000000,
259 0x00000000, 0x00000000, 0x00000000, 0x00000000,
260 0x00000000, 0x00000000, 0x00000000, 0x00000000,
261 0x00000000, 0x00000000, 0x00000000, 0x00000000,
262 0x00000000, 0x00000000, 0x00000000, 0x00000000,
263 0x00000000, 0x00000000, 0x00000000, 0x00000000,
264
265 0xBFD09B84, 0xB58C74B1, 0x3C649EAC, 0xBE7125A0,
266 0xBFE04818, 0x9E7A4AF5, 0xBC7247BB, 0x9C14B6B2,
267 0xBFEC9082, 0x7462A177, 0xBC65A4D5, 0x2B14D05A,
268 0xBF9A71A4, 0xBDE00D7D, 0xBC3D297D, 0x6BADB471,
269 0x3F789A87, 0xB48FE138, 0x3BF9AE10, 0x1027D54D,
270 0xBF549926, 0x2376DBFA, 0x3BDC0F4D, 0x2A750CAB,
271 0x3F2C4E86, 0xAA32FCDB, 0xBBCBD7A9, 0xE43684BF,
272 0xBEF790BC, 0x075477C9, 0x3B89898A, 0xADCD6E9C,
273 0xBEC63409, 0x8491E3A5, 0xBB572A20, 0x12F684BF,
274 0x3EC3E354, 0x0DCFCA75, 0x3B485BE3, 0x87E6B74C,
275 0xBEAB790E, 0xE236BFDB, 0x3B40DEF2, 0x9161F9AE,
276 0x3E89DA01, 0x98D7CE6C, 0xBB2327A2, 0xC3F35BA6,
277 0xBE5DEF8D, 0x200AB94C, 0x3AF95BA7, 0x81948B10,
278 0xBE129E35, 0x62080668, 0x3AB6F684, 0xBDA12F68,
279 0x3E24ACF9, 0x5B57BC54, 0x3AC55555, 0x55555556,
280 0xBE10CB5E, 0x9DE40846, 0x3AAC71C7, 0x1C71C71F,
281 0x3DF1C294, 0xD79A2DA8, 0xBA9F9ADD, 0x3C0CA45A,
282 0xBDC806EC, 0x95F17FF3, 0xBA66E9E0, 0x6522C3F3,
283 0x3D63D790, 0x6BE97B42, 0x3A07B425, 0xED097B40,
284 0x3D8A0F87, 0x44E09E06, 0x3A248B0F, 0xCD6E9E04,
285 0xBD780C79, 0x96291620, 0x3A1948B0, 0xFCD6E9E0,
286 0x3D5B9C07, 0xD9F03291, 0x39F87E6B, 0x74F03291,
287 0xBD34D565, 0xE06522C4, 0x39A948B0, 0xFCD6E9DF,
288 0x3CF1D800, 0xD6E9E065, 0x398161F9, 0xADD3C0C9,
289 0x3CF13B9A, 0xDD3C0CA4, 0x39961F9A, 0xDD3C0CA4,
290 0xBCE41FE1, 0xF9ADD3C1, 0x397ADD3C, 0x0CA45880,
291 0x3CC7D948, 0xB0FCD6EA, 0xB94F9ADD, 0x3C0CA458,
292 0x00000000, 0x00000000, 0x00000000, 0x00000000,
293 0x00000000, 0x00000000, 0x00000000, 0x00000000,
294
295 0xBF7E4EAA, 0x825F8588, 0x3C0E9883, 0x7D311FB0,
296 0xBF9B20DB, 0x38454E52, 0x3C282460, 0xA5A22510,
297 0xBFE92D49, 0x5DE062C7, 0x3C842242, 0x550B45F8,
298 0xBFAE856A, 0x106722F3, 0xBC3E5242, 0x6EC1567E,
299 0x3F8CF900, 0xD71B5137, 0xBC29B9D7, 0x4D921F9C,
300 0xBF630F14, 0x2DC98724, 0x3BEFD901, 0xC6097B43,
301 0x3EE6A533, 0x044E582C, 0x3B653307, 0x1C71C6A0,
302 0x3F261616, 0xBB9808F2, 0x3BC00838, 0xA4587E6A,
303 0xBF12A42E, 0x193904D9, 0x3B7FAA5E, 0xD097B432,
304 0x3EEFFC81, 0x1FC57602, 0xBB657CD6, 0xE9E0651E,
305 0x3E8108C1, 0x35437559, 0xBAFF9ADD, 0x3C0CA484,
306 0xBEBA9B0B, 0xE9B78F13, 0x3B555555, 0x55555554,
307 0x3EA7C3AD, 0xEB179410, 0xBB3948B0, 0xFCD6E9DF,
308 0xBE853AAC, 0x329694F0, 0xBB1948B0, 0xFCD6E9E2,
309 0xBE27591E, 0x10C71C72, 0x3ABC71C7, 0x1C71C700,
310 0x3E53D507, 0xBCC71C72, 0xBAEC71C7, 0x1C71C720,
311 0xBE4220D2, 0xDA12F685, 0x3AE097B4, 0x25ED097C,
312 0x3E2224EC, 0xBDA12F68, 0x3AC2F684, 0xBDA12F68,
313 0xBDF19BC0, 0xCA4587E7, 0x3A922C3F, 0x35BA781A
314 // Last block shorter to save space
315};
316
317typedef double DD[2];
318typedef DD DDarray[29];
319
320static const Double ca = {{0x3FF20DD7, 0x50429B6D}};
321static const Double cb = {{0x3C71AE3A, 0x914FED80}};
322static const Double cc = {{0xB903CBBE, 0xBF65F145}};
323
324static long double ErfCommon( double dhead, double dtail, int ierf )
325{
326 double temp, temp2, temp3, argsq, argsq2, sum, suml, suminf, reslow, resmid, restiny;
327 double q, p, prodlow, addend, res, prod1,prod2,px, center, adjust, times;
328 double oarg, oarg2, a1, a2, a3, ab, ac, asq1,asq2,asq, arg2,argx;
329 double rx, arg, as, asx, resx, qx, exp;
330 DblDbl resdd;
331 int i, ir, is, ie;
332
333 double *erfxtr = (double *)erfxtrData;
334 DD *erfpow = (DD *)erfpowData;
335 DDarray *coeff = (DDarray *)(coeffData - 116);
336
337 if (fabs(dhead) <= 0.75) { // power series expansion
338
339 temp = (2.0*dhead)*dtail; // series is odd series-
340 argsq = __FMADD( dhead, dhead, temp ); // dhead*dhead + temp; // most of expansion is
341 argsq2 = __FMSUB( dhead, dhead, argsq + temp ) + // dhead*dhead - argsq+temp +
342 __FMSUB( (2.0*dhead), dtail, temp ); // ((2.0*dhead)*dtail-temp); // is in terms of arg^2
343 sum = erfpow[21][0]; // sum of series init
344 suml = erfpow[21][1];
345 for (i = 20; i >= 3; --i) { // all but the 3 h.o.
346 temp = __FMADD( sum, argsq, erfpow[i][0] ); // erfpow[i][0] + sum*argsq; // terms executed in
347 // q&d quad precision
348 /* suml = erfpow[i][0] - temp + sum*argsq +
349 (erfpow[i][1] + sum*argsq2 + suml*argsq); */
350 suml = __FMADD( sum, argsq, erfpow[i][0] - temp ) +
351 __FMADD( suml, argsq, __FMADD( sum, argsq2, erfpow[i][1] ) );
352 sum = temp;
353 }
354 suminf = 0.0e0;
355 for (i = 2; i >= 0; --i) {
356 prodlow = __FMADD( suml, argsq, (sum*argsq2) ); // (suml*argsq) + (sum*argsq2); // mult. by arg^2
357 /* addend = (suml*argsq - (suml*argsq)) + (sum*argsq2 - (sum*argsq2)) +
358 ((suml*argsq) - prodlow + (sum*argsq2)) + suminf*argsq; */
359 addend = __FMSUB( suml, argsq, (suml*argsq) ) + __FMSUB( sum, argsq2, (sum*argsq2) ) +
360 __FMADD( suminf, argsq, __FMADD( sum, argsq2, __FMSUB( suml, argsq, prodlow ) ) );
361 temp = __FMADD( sum, argsq, prodlow ); // sum*argsq + prodlow; // last 3 terms computed
362 temp2 = __FADD( __FMSUB( sum, argsq, temp), prodlow );// sum*argsq - temp + prodlow; // with greater care
363 temp3 = __FMSUB( sum, argsq, temp) - temp2 + prodlow + addend;//(sum*argsq - temp) - temp2 + prodlow + addend;
364 res = __FADD( erfpow[i][0], temp );
365 reslow = (erfpow[i][0] - res + temp);
366 resmid = __FADD( erfpow[i][1], temp2 );
367 restiny = erfpow[i][1] - resmid + temp2 + erfxtr[i] + temp3;
368 p = __FADD( reslow, resmid ); // sum of middle terms
369 q = reslow - p + resmid; // possible residual
370 sum = __FADD( res, p );
371 suml = __FADD( (res - sum + p), (q + restiny) );
372 suminf = (res - sum + p) - suml + (q + restiny);
373 }
374
375 //********************************************************************
376 // To complete the job, the series (so far summed in even powers), *
377 // must be multiplied by the argument and by 2/PI^.5. *
378 // These products are executed in sextuple precision to avoid *
379 // excessive rounding errors. *
380 //********************************************************************
381
382 _d3mul( sum, suml, suminf, dhead, dtail, 0.0, &prod1, &prod2, &px );
383 _d3mul( prod1, prod2, px, ca.f, cb.f, cc.f, &(resdd.d[0]), &(resdd.d[1]), &qx );
384
385 if (ierf == 0) {
386 // erfc(x) = 1 - erf(x)
387 _d3add( 1.0, 0.0, 0.0, -resdd.d[0], -resdd.d[1], -qx,
388 &(resdd.d[0]), &(resdd.d[1]), &qx );
389 }
390 return resdd.f;
391 }
392 else if (fabs(dhead) <= 2.0) {
393
394 //********************************************************************
395 // For arguments 1 < |x| <= 2.0, erfc is computed as: *
396 // *
397 // erfc(x)= e^(-x^2+series) *
398 // *
399 // WHERE: "series" is a polynomial approximation devised to *
400 // make the above formula correct to within 105 bits of *
401 // precision. The actual series is in powers of *
402 // |x|-center. *
403 //********************************************************************
404
405 if (fabs(dhead) < 1.00) { // 3 subintervals used:
406 center = 0.625; // 0.75 < x < 1.0
407 adjust = 0.28125; // 1.0 <= x < 1.25
408 times = 2.0; // 1.25 <= x <= 2.0
409 ir = 8;
410 is = 18;
411 ie = 15;
412 }
413 else if (fabs(dhead) < 1.25) {
414 center = 0.8125;
415 adjust = 0.28125;
416 times = 2.0;
417 ir = 6;
418 is = 19;
419 ie = 16;
420 }
421 else {
422 center = 1.5;
423 adjust = 1.375;
424 times = 3.0;
425 ir = 7;
426 is = 26;
427 ie = 14;
428 }
429
430 if (dhead < 0.0) {
431 oarg = dhead;
432 oarg2 = dtail;
433 temp = -center - dhead; // Compute |x| - center
434 arg2 = -dtail;
435 }
436 else {
437 temp = dhead - center;
438 arg2 = dtail;
439 oarg = -dhead;
440 oarg2 = -dtail;
441 }
442 arg = __FADD( temp, arg2 );
443 arg2 = temp - arg + arg2;
444 sum = coeff[ir][is][0];
445 suml = coeff[ir][is][1];
446 for (i = is-1; i >= ie; --i) { // computed using q&d
447 temp = __FMADD( sum, arg, coeff[ir][i][0] ); // coeff[ir][i][0] + sum*arg; // quad precision
448 /* suml = coeff[ir][i][0] - temp + sum*arg +
449 (coeff[ir][i][1] + sum*arg2 + suml*arg); */
450 suml = __FMADD( sum, arg, coeff[ir][i][0] - temp ) +
451 __FMADD( suml, arg, __FMADD( sum, arg2, coeff[ir][i][1] ) );
452 sum = temp;
453 }
454 suminf = 0.0e0; // remaining terms require full
455 for (i = ie-1; i >= 0; --i) { // sextuple precision arithmetic
456 _d3mul( sum, suml, suminf, arg, arg2, 0.0, &prod1, &prod2, &px );
457 _d3add( prod1, prod2, px, coeff[ir][i][0], coeff[ir][i][1], 0.0,
458 &sum, &suml, &suminf );
459 // NOTE: The original FORTRAN (apparently in error) was:
460 // _sum=_q_a6(%val(_prod),%val(px),%val(coeff(1,i,ir)),
461 // + %val(coeff(2,i,ir)),suminf) <-- wrong number of parameters!
462 }
463 a2 = times*oarg2;
464 a3 = __FMSUB( times, oarg2, a2 ); // times*oarg2 - a2; // series is adjusted by
465 a1 = oarg*times; // adding times*x-adjust
466 ab = __FADD( __FMSUB( times, oarg2, a1 ), a2 ); // oarg*times - a1 + a2;
467 ac = __FMSUB( times, oarg2, a1 ) - ab + a2 + a3; // (oarg*times - a1) - ab + a2 + a3;
468 _d3add( sum, suml, suminf, adjust, 0.0, 0.0, &sum, &suml, &as );
469 _d3add( sum, suml, as, a1, ab, ac, &sum, &suml, &asx );
470 resdd.f = _ExpInnerLD(sum, suml, asx, &resx, 4);
471
472 // Erfc(|x|) done!
473 if (ierf == 0) {
474 if (dhead <= 0.0) // For pos args, done!
475 // Erfc(-x) = 2 - Erfc(x)
476 _d3add( 2.0, 0.0, 0.0, -resdd.d[0], -resdd.d[1], -resx,
477 &(resdd.d[0]), &(resdd.d[1]), &qx );
478 }
479 else { // Caller wanted Erf.
480 // Use Erf(x) = 1 - Erfc(x)
481 if (dhead > 0.0)
482 _d3add( 1.0, 0.0, 0.0, -resdd.d[0], -resdd.d[1], -resx,
483 &(resdd.d[0]), &(resdd.d[1]), &qx );
484 else
485 _d3add(-1.0, 0.0, 0.0, resdd.d[0], resdd.d[1], resx,
486 &(resdd.d[0]), &(resdd.d[1]), &qx );
487 }
488 return resdd.f;
489 }
490 else {
491
492 //********************************************************************
493 // For all remaining arguments, ERFC(x) is calculated as follows: *
494 // *
495 // erfc(x)=e^(-x^2) / x*cf, *
496 // *
497 // WHERE: cf is a continued fraction that converges in the left *
498 // half plane. *
499 // *
500 // VARIABLE NAME MEANING *
501 // cf over various intervals, is *
502 // approximated by polynomials *
503 // *
504 // ir identifies the polynomial *
505 // *
506 // is the degree of the polynomial *
507 // *
508 // ie the lower bound after which *
509 // sextuple precision must be used *
510 // *
511 // The approximating polynomial is in powers of (1/x)^2 *
512 //********************************************************************
513
514 _d3mul( dhead, dtail, 0.0, dhead, dtail, 0.0,
515 &asq1, &asq2, &asq ); // arg^2
516 _d3div( 1.0, 0.0, 0.0, asq1, asq2, asq,
517 &arg, &arg2, &argx ); // 1/arg^2
518 if (arg >= 0.18e0) { // .18 <= _arg < .25
519 ir = 1; // which polynomial
520 is = 21; // How many terms
521 ie = 10; // Last term done q&d
522 }
523 else if (arg >= 0.11e0) {
524 ir = 2; // .11 <= _arg < .18
525 is = 23;
526 ie = 12;
527 }
528 else if (arg >= 0.06e0) {
529 ir = 3; // .06 <= _arg < .11
530 is = 22;
531 ie = 6;
532 if (ierf != 1) ie = 11;
533 }
534 else if (arg >= 0.01e0) {
535 ir = 4; // .01 <= _arg < .06
536 is = 28;
537 ie = 10;
538 }
539 else if (arg >= 0.00138e0) {
540 ir = 5; // .00138 <= _arg < .01
541 is = 17;
542 ie = 4;
543 if (ierf != 1) ie = 8;
544 }
545 else {
546
547 //***********************************************************
548 // For larger arguments, erfc is so small that it cannot *
549 // be represented, so the value of zero is used for erfc *
550 //***********************************************************
551
552 resdd.d[0] = 0.0;
553 resdd.d[1] = 0.0;
554 rx = 0.0;
555 if (dhead == dhead) goto commontail;
556
557 // For NaN inputs ->NaNQ
558 resdd.d[0] = dhead;
559 return resdd.f;
560 }
561 if ((ierf == 1) && (arg < .01)) { // For erf(x), x > 10,
562 resdd.d[0] = 0.0;
563 resdd.d[1] = 0.0; // the result is 1.
564 rx = 0.0; // (to save time)
565 goto commontail;
566 }
567 sum = coeff[ir][is][0];
568 suml = coeff[ir][is][1]; // q&d quad prec.
569 for (i = is-1; i >= ie; --i) { // part of polynomial
570 temp = __FMADD( sum, arg, coeff[ir][i][0] ); // coeff[ir][i][0] + sum*arg; // evaluation
571 /* suml = coeff[ir][i][0] - temp + sum*arg +
572 (coeff[ir][i][1] + sum*arg2 + suml*arg); */
573 suml = __FMADD( sum, arg, coeff[ir][i][0] - temp ) +
574 __FMADD( suml, arg, __FMADD( sum, arg2, coeff[ir][i][1] ) );
575 sum = temp;
576 }
577 suminf = 0.0e0; // sextuple prec eval of
578 // remaining terms of the
579 for (i = ie-1; i >= 0; --i) { // approximating polynomial
580 _d3mul( sum, suml, suminf, arg, arg2, argx, &prod1, &prod2, &px);
581 _d3add( prod1, prod2, px, coeff[ir][i][0], coeff[ir][i][1], 0.0,
582 &sum, &suml, &suminf );
583 // NOTE: The original FORTRAN (apparently in error) was:
584 // _sum=_q_a6(%val(_prod),%val(px),%val(coeff(1,i,ir)),
585 // + %val(coeff(2,i,ir)),suminf) <-- wrong number of parameters!
586 }
587 _d3mul( sum, suml, suminf, ca.f, cb.f, cc.f,
588 &prod1, &prod2, &px ); // multiply by 2/sqrt(pi)
589 resdd.f = _ExpInnerLD(-asq1, -asq2, -asq, &exp, 4);
590 _d3mul( prod1, prod2, px, resdd.d[0], resdd.d[1], exp,
591 &prod1, &prod2, &px );
592 _d3div( prod1, prod2, px, dhead, dtail, 0.0,
593 &(resdd.d[0]), &(resdd.d[1]), &rx );
594 resdd.f = 0.5 * resdd.f;
595 rx = 0.5*rx;
596
597 commontail:
598
599 if (dhead > 0.0) { // to return ERF or ERFC
600 if (ierf == 1) // for pos and neg args
601 _d3add( 1.0, 0.0, 0.0, -resdd.d[0], -resdd.d[1], -rx,
602 &(resdd.d[0]), &(resdd.d[1]), &qx );
603 }
604 else {
605 if (ierf == 1)
606 _d3add(-1.0, 0.0, 0.0, -resdd.d[0], -resdd.d[1], -rx,
607 &(resdd.d[0]), &(resdd.d[1]), &qx );
608 else
609 _d3add( 2.0, 0.0, 0.0, resdd.d[0], resdd.d[1], rx,
610 &(resdd.d[0]), &(resdd.d[1]), &qx );
611 }
612 return resdd.f;
613 }
614}
615
616
617long double erfl( long double x )
618{
619 double head, fenv;
620 DblDbl t;
621
622 FEGETENVD(fenv);
623 FESETENVD(0.0);
624
625 t.f = x;
626
627 t.f = ErfCommon(t.d[0], t.d[1], 1);
628
629 head = __FADD( t.d[0], t.d[1] );
630 t.d[1] = (t.d[0] - head) + t.d[1]; // cannonize tail
631 t.d[0] = head; // cannonize head
632
633 FESETENVD(fenv);
634 return t.f;
635}
636
637long double erfcl( long double x )
638{
639 double head, fenv;
640 DblDbl t;
641
642 FEGETENVD(fenv);
643 FESETENVD(0.0);
644
645 t.f = x;
646
647 t.f = ErfCommon(t.d[0], t.d[1], 0);
648
649 head = __FADD( t.d[0], t.d[1] );
650 t.d[1] = (t.d[0] - head) + t.d[1]; // cannonize tail
651 t.d[0] = head; // cannonize head
652
653 FESETENVD(fenv);
654 return t.f;
655}
656#endif