Graph Framework
Loading...
Searching...
No Matches
special_functions.hpp
Go to the documentation of this file.
1//------------------------------------------------------------------------------
6//------------------------------------------------------------------------------
7
8#ifndef special_functions_h
9#define special_functions_h
10
11#ifdef CUDA_DEVICE_CODE
12#include <cuda/std/complex>
13
14using namespace cuda::std;
15
17template<typename T>
18using complex_type = complex<T>;
19
20#else
21#include <iostream>
22#include <complex>
23#include <cfloat>
24#include <cmath>
25#include <cstdint>
26#include <limits>
27
28using namespace std;
29
31#if __cplusplus >= 202002L
32template<std::floating_point T>
33#else
34template<typename T>
35#endif
36using complex_type = std::complex<T>;
37
38#endif
39
41namespace special {
42
44#if __cplusplus >= 202002L
45 template<std::floating_point T>
46#else
47 template<typename T>
48#endif
49 constexpr complex_type<T> i(static_cast<T> (0),
50 static_cast<T> (1));
51
52//------------------------------------------------------------------------------
59//------------------------------------------------------------------------------
60#if __cplusplus >= 202002L
61 template<std::floating_point T>
62#else
63 template<typename T>
64#endif
65 complex_type<T> erfcx(complex_type<T> z) {
66 return w(i<T>*z);
67 }
68
69//------------------------------------------------------------------------------
76//------------------------------------------------------------------------------
77#if __cplusplus >= 202002L
78 template<std::floating_point T>
79#else
80 template<typename T>
81#endif
82 T sq(T x) {
83 return x*x;
84 }
85
86//------------------------------------------------------------------------------
100//------------------------------------------------------------------------------
101#if __cplusplus >= 202002L
102 template<std::floating_point T>
103#else
104 template<typename T>
105#endif
106 T w_im_y100(T x) {
107 const T y100 = static_cast<T> (100)/(static_cast<T> (1) + x);
108 switch (static_cast<uint8_t> (y100)) {
109 case 0: {
110 const T t = static_cast<T> (2)*y100 - static_cast<T> (1);
111 return static_cast<T> (0.28351593328822191546E-2) + (static_cast<T> (0.28494783221378400759E-2) + (static_cast<T> (0.14427470563276734183E-4) + (static_cast<T> (0.10939723080231588129E-6) + (static_cast<T> (0.92474307943275042045E-9) + (static_cast<T> (0.89128907666450075245E-11) + static_cast<T> (0.92974121935111111110E-13)*t)*t)*t)*t)*t)*t;
112 }
113 case 1: {
114 const T t = static_cast<T> (2)*y100 - static_cast<T> (3);
115 return static_cast<T> (0.85927161243940350562E-2) + (static_cast<T> (0.29085312941641339862E-2) + (static_cast<T> (0.15106783707725582090E-4) + (static_cast<T> (0.11716709978531327367E-6) + (static_cast<T> (0.10197387816021040024E-8) + (static_cast<T> (0.10122678863073360769E-10) + static_cast<T> (0.10917479678400000000E-12)*t)*t)*t)*t)*t)*t;
116 }
117 case 2: {
118 const T t = static_cast<T> (2)*y100 - static_cast<T> (5);
119 return static_cast<T> (0.14471159831187703054E-1) + (static_cast<T> (0.29703978970263836210E-2) + (static_cast<T> (0.15835096760173030976E-4) + (static_cast<T> (0.12574803383199211596E-6) + (static_cast<T> (0.11278672159518415848E-8) + (static_cast<T> (0.11547462300333495797E-10) + static_cast<T> (0.12894535335111111111E-12)*t)*t)*t)*t)*t)*t;
120 }
121 case 3: {
122 const T t = static_cast<T> (2)*y100 - static_cast<T> (7);
123 return static_cast<T> (0.20476320420324610618E-1) + (0.30352843012898665856E-2 + (0.16617609387003727409E-4 + (0.13525429711163116103E-6 + (0.12515095552507169013E-8 + (0.13235687543603382345E-10 + 0.15326595042666666667E-12*t)*t)*t)*t)*t)*t;
124 }
125 case 4: {
126 const T t = static_cast<T> (2)*y100 - static_cast<T> (9);
127 return static_cast<T> (0.26614461952489004566E-1) + (0.31034189276234947088E-2 + (0.17460268109986214274E-4 + (0.14582130824485709573E-6 + (0.13935959083809746345E-8 + (0.15249438072998932900E-10 + 0.18344741882133333333E-12*t)*t)*t)*t)*t)*t;
128 }
129 case 5: {
130 const T t = static_cast<T> (2)*y100 - static_cast<T> (11);
131 return static_cast<T> (0.32892330248093586215E-1) + (0.31750557067975068584E-2 + (0.18369907582308672632E-4 + (0.15761063702089457882E-6 + (0.15577638230480894382E-8 + (0.17663868462699097951E-10 + (0.22126732680711111111E-12 + 0.30273474177737853668E-14*t)*t)*t)*t)*t)*t)*t;
132 }
133 case 6: {
134 const T t = static_cast<T> (2)*y100 - static_cast<T> (13);
135 return static_cast<T> (0.39317207681134336024E-1) + (0.32504779701937539333E-2 + (0.19354426046513400534E-4 + (0.17081646971321290539E-6 + (0.17485733959327106250E-8 + (0.20593687304921961410E-10 + (0.26917401949155555556E-12 + 0.38562123837725712270E-14*t)*t)*t)*t)*t)*t)*t;
136 }
137 case 7: {
138 const T t = static_cast<T> (2)*y100 - static_cast<T> (15);
139 return static_cast<T> (0.45896976511367738235E-1) + (0.33300031273110976165E-2 + (0.20423005398039037313E-4 + (0.18567412470376467303E-6 + (0.19718038363586588213E-8 + (0.24175006536781219807E-10 + (0.33059982791466666666E-12 + 0.49756574284439426165E-14*t)*t)*t)*t)*t)*t)*t;
140 }
141 case 8: {
142 const T t = static_cast<T> (2)*y100 - static_cast<T> (17);
143 return static_cast<T> (0.52640192524848962855E-1) + (0.34139883358846720806E-2 + (0.21586390240603337337E-4 + (0.20247136501568904646E-6 + (0.22348696948197102935E-8 + (0.28597516301950162548E-10 + (0.41045502119111111110E-12 + 0.65151614515238361946E-14*t)*t)*t)*t)*t)*t)*t;
144 }
145 case 9: {
146 const T t = static_cast<T> (2)*y100 - static_cast<T> (19);
147 return static_cast<T> (0.59556171228656770456E-1) + (0.35028374386648914444E-2 + (0.22857246150998562824E-4 + (0.22156372146525190679E-6 + (0.25474171590893813583E-8 + (0.34122390890697400584E-10 + (0.51593189879111111110E-12 + 0.86775076853908006938E-14*t)*t)*t)*t)*t)*t)*t;
148 }
149 case 10: {
150 const T t = static_cast<T> (2)*y100 - static_cast<T> (21);
151 return static_cast<T> (0.66655089485108212551E-1) + (0.35970095381271285568E-2 + (0.24250626164318672928E-4 + (0.24339561521785040536E-6 + (0.29221990406518411415E-8 + (0.41117013527967776467E-10 + (0.65786450716444444445E-12 + 0.11791885745450623331E-13*t)*t)*t)*t)*t)*t)*t;
152 }
153 case 11: {
154 const T t = static_cast<T> (2)*y100 - static_cast<T> (23);
155 return static_cast<T> (0.73948106345519174661E-1) + (0.36970297216569341748E-2 + (0.25784588137312868792E-4 + (0.26853012002366752770E-6 + (0.33763958861206729592E-8 + (0.50111549981376976397E-10 + (0.85313857496888888890E-12 + 0.16417079927706899860E-13*t)*t)*t)*t)*t)*t)*t;
156 }
157 case 12: {
158 const T t = static_cast<T> (2)*y100 - static_cast<T> (25);
159 return static_cast<T> (0.81447508065002963203E-1) + (0.38035026606492705117E-2 + (0.27481027572231851896E-4 + (0.29769200731832331364E-6 + (0.39336816287457655076E-8 + (0.61895471132038157624E-10 + (0.11292303213511111111E-11 + 0.23558532213703884304E-13*t)*t)*t)*t)*t)*t)*t;
160 }
161 case 13: {
162 const T t = static_cast<T> (2)*y100 - static_cast<T> (27);
163 return static_cast<T> (0.89166884027582716628E-1) + (0.39171301322438946014E-2 + (0.29366827260422311668E-4 + (0.33183204390350724895E-6 + (0.46276006281647330524E-8 + (0.77692631378169813324E-10 + (0.15335153258844444444E-11 + 0.35183103415916026911E-13*t)*t)*t)*t)*t)*t)*t;
164 }
165 case 14: {
166 const T t = static_cast<T> (2)*y100 - static_cast<T> (29);
167 return static_cast<T> (0.97121342888032322019E-1) + (0.40387340353207909514E-2 + (0.31475490395950776930E-4 + (0.37222714227125135042E-6 + (0.55074373178613809996E-8 + (0.99509175283990337944E-10 + (0.21552645758222222222E-11 + 0.55728651431872687605E-13*t)*t)*t)*t)*t)*t)*t;
168 }
169 case 15: {
170 const T t = static_cast<T> (2)*y100 - static_cast<T> (31);
171 return static_cast<T> (0.10532778218603311137) + (0.41692873614065380607E-2 + (0.33849549774889456984E-4 + (0.42064596193692630143E-6 + (0.66494579697622432987E-8 + (0.13094103581931802337E-9 + (0.31896187409777777778E-11 + 0.97271974184476560742E-13*t)*t)*t)*t)*t)*t)*t;
172 }
173 case 16: {
174 const T t = static_cast<T> (2)*y100 - static_cast<T> (33);
175 return static_cast<T> (0.11380523107427108222) + (0.43099572287871821013E-2 + (0.36544324341565929930E-4 + (0.47965044028581857764E-6 + (0.81819034238463698796E-8 + (0.17934133239549647357E-9 + (0.50956666166186293627E-11 + (0.18850487318190638010E-12 + 0.79697813173519853340E-14*t)*t)*t)*t)*t)*t)*t)*t;
176 }
177 case 17: {
178 const T t = static_cast<T> (2)*y100 - static_cast<T> (35);
179 return static_cast<T> (0.12257529703447467345) + (0.44621675710026986366E-2 + (0.39634304721292440285E-4 + (0.55321553769873381819E-6 + (0.10343619428848520870E-7 + (0.26033830170470368088E-9 + (0.87743837749108025357E-11 + (0.34427092430230063401E-12 + 0.10205506615709843189E-13*t)*t)*t)*t)*t)*t)*t)*t;
180 }
181 case 18: {
182 const T t = static_cast<T> (2)*y100 - static_cast<T> (37);
183 return 0.13166276955656699478 + (0.46276970481783001803E-2 + (0.43225026380496399310E-4 + (0.64799164020016902656E-6 + (0.13580082794704641782E-7 + (0.39839800853954313927E-9 + (0.14431142411840000000E-10 + 0.42193457308830027541E-12*t)*t)*t)*t)*t)*t)*t;
184 }
185 case 19: {
186 const T t = static_cast<T> (2)*y100 - static_cast<T> (39);
187 return 0.14109647869803356475 + (0.48088424418545347758E-2 + (0.47474504753352150205E-4 + (0.77509866468724360352E-6 + (0.18536851570794291724E-7 + (0.60146623257887570439E-9 + (0.18533978397305276318E-10 + (0.41033845938901048380E-13 - 0.46160680279304825485E-13*t)*t)*t)*t)*t)*t)*t)*t;
188 }
189 case 20: {
190 const T t = static_cast<T> (2)*y100 - static_cast<T> (41);
191 return 0.15091057940548936603 + (0.50086864672004685703E-2 + (0.52622482832192230762E-4 + (0.95034664722040355212E-6 + (0.25614261331144718769E-7 + (0.80183196716888606252E-9 + (0.12282524750534352272E-10 + (-0.10531774117332273617E-11 - 0.86157181395039646412E-13*t)*t)*t)*t)*t)*t)*t)*t;
192 }
193 case 21: {
194 const T t = static_cast<T> (2)*y100 - static_cast<T> (43);
195 return 0.16114648116017010770 + (0.52314661581655369795E-2 + (0.59005534545908331315E-4 + (0.11885518333915387760E-5 + (0.33975801443239949256E-7 + (0.82111547144080388610E-9 + (-0.12357674017312854138E-10 + (-0.24355112256914479176E-11 - 0.75155506863572930844E-13*t)*t)*t)*t)*t)*t)*t)*t;
196 }
197 case 22: {
198 const T t = static_cast<T> (2)*y100 - static_cast<T> (45);
199 return 0.17185551279680451144 + (0.54829002967599420860E-2 + (0.67013226658738082118E-4 + (0.14897400671425088807E-5 + (0.40690283917126153701E-7 + (0.44060872913473778318E-9 + (-0.52641873433280000000E-10 - 0.30940587864543343124E-11*t)*t)*t)*t)*t)*t)*t;
200 }
201 case 23: {
202 const T t = static_cast<T> (2)*y100 - static_cast<T> (47);
203 return 0.18310194559815257381 + (0.57701559375966953174E-2 + (0.76948789401735193483E-4 + (0.18227569842290822512E-5 + (0.41092208344387212276E-7 + (-0.44009499965694442143E-9 + (-0.92195414685628803451E-10 + (-0.22657389705721753299E-11 + 0.10004784908106839254E-12*t)*t)*t)*t)*t)*t)*t)*t;
204 }
205 case 24: {
206 const T t = static_cast<T> (2)*y100 - static_cast<T> (49);
207 return 0.19496527191546630345 + (0.61010853144364724856E-2 + (0.88812881056342004864E-4 + (0.21180686746360261031E-5 + (0.30652145555130049203E-7 + (-0.16841328574105890409E-8 + (-0.11008129460612823934E-9 + (-0.12180794204544515779E-12 + 0.15703325634590334097E-12*t)*t)*t)*t)*t)*t)*t)*t;
208 }
209 case 25: {
210 const T t = static_cast<T> (2)*y100 - static_cast<T> (51);
211 return 0.20754006813966575720 + (0.64825787724922073908E-2 + (0.10209599627522311893E-3 + (0.22785233392557600468E-5 + (0.73495224449907568402E-8 + (-0.29442705974150112783E-8 + (-0.94082603434315016546E-10 + (0.23609990400179321267E-11 + 0.14141908654269023788E-12*t)*t)*t)*t)*t)*t)*t)*t;
212 }
213 case 26: {
214 const T t = static_cast<T> (2)*y100 - static_cast<T> (53);
215 return 0.22093185554845172146 + (0.69182878150187964499E-2 + (0.11568723331156335712E-3 + (0.22060577946323627739E-5 + (-0.26929730679360840096E-7 + (-0.38176506152362058013E-8 + (-0.47399503861054459243E-10 + (0.40953700187172127264E-11 + 0.69157730376118511127E-13*t)*t)*t)*t)*t)*t)*t)*t;
216 }
217 case 27: {
218 const T t = static_cast<T> (2)*y100 - static_cast<T> (55);
219 return 0.23524827304057813918 + (0.74063350762008734520E-2 + (0.12796333874615790348E-3 + (0.18327267316171054273E-5 + (-0.66742910737957100098E-7 + (-0.40204740975496797870E-8 + (0.14515984139495745330E-10 + (0.44921608954536047975E-11 - 0.18583341338983776219E-13*t)*t)*t)*t)*t)*t)*t)*t;
220 }
221 case 28: {
222 const T t = static_cast<T> (2)*y100 - static_cast<T> (57);
223 return 0.25058626331812744775 + (0.79377285151602061328E-2 + (0.13704268650417478346E-3 + (0.11427511739544695861E-5 + (-0.10485442447768377485E-6 + (-0.34850364756499369763E-8 + (0.72656453829502179208E-10 + (0.36195460197779299406E-11 - 0.84882136022200714710E-13*t)*t)*t)*t)*t)*t)*t)*t;
224 }
225 case 29: {
226 const T t = static_cast<T> (2)*y100 - static_cast<T> (59);
227 return 0.26701724900280689785 + (0.84959936119625864274E-2 + (0.14112359443938883232E-3 + (0.17800427288596909634E-6 + (-0.13443492107643109071E-6 + (-0.23512456315677680293E-8 + (0.11245846264695936769E-9 + (0.19850501334649565404E-11 - 0.11284666134635050832E-12*t)*t)*t)*t)*t)*t)*t)*t;
228 }
229 case 30: {
230 const T t = static_cast<T> (2)*y100 - static_cast<T> (61);
231 return 0.28457293586253654144 + (0.90581563892650431899E-2 + (0.13880520331140646738E-3 + (-0.97262302362522896157E-6 + (-0.15077100040254187366E-6 + (-0.88574317464577116689E-9 + (0.12760311125637474581E-9 + (0.20155151018282695055E-12 - 0.10514169375181734921E-12*t)*t)*t)*t)*t)*t)*t)*t;
232 }
233 case 31: {
234 const T t = static_cast<T> (2)*y100 - static_cast<T> (63);
235 return 0.30323425595617385705 + (0.95968346790597422934E-2 + (0.12931067776725883939E-3 + (-0.21938741702795543986E-5 + (-0.15202888584907373963E-6 + (0.61788350541116331411E-9 + (0.11957835742791248256E-9 + (-0.12598179834007710908E-11 - 0.75151817129574614194E-13*t)*t)*t)*t)*t)*t)*t)*t;
236 }
237 case 32: {
238 const T t = static_cast<T> (2)*y100 - static_cast<T> (65);
239 return 0.32292521181517384379 + (0.10082957727001199408E-1 + (0.11257589426154962226E-3 + (-0.33670890319327881129E-5 + (-0.13910529040004008158E-6 + (0.19170714373047512945E-8 + (0.94840222377720494290E-10 + (-0.21650018351795353201E-11 - 0.37875211678024922689E-13*t)*t)*t)*t)*t)*t)*t)*t;
240 }
241 case 33: {
242 const T t = static_cast<T> (2)*y100 - static_cast<T> (67);
243 return 0.34351233557911753862 + (0.10488575435572745309E-1 + (0.89209444197248726614E-4 + (-0.43893459576483345364E-5 + (-0.11488595830450424419E-6 + (0.28599494117122464806E-8 + (0.61537542799857777779E-10 - 0.24935749227658002212E-11*t)*t)*t)*t)*t)*t)*t;
244 }
245 case 34: {
246 const T t = static_cast<T> (2)*y100 - static_cast<T> (69);
247 return 0.36480946642143669093 + (0.10789304203431861366E-1 + (0.60357993745283076834E-4 + (-0.51855862174130669389E-5 + (-0.83291664087289801313E-7 + (0.33898011178582671546E-8 + (0.27082948188277716482E-10 + (-0.23603379397408694974E-11 + 0.19328087692252869842E-13*t)*t)*t)*t)*t)*t)*t)*t;
248 }
249 case 35: {
250 const T t = static_cast<T> (2)*y100 - static_cast<T> (71);
251 return 0.38658679935694939199 + (0.10966119158288804999E-1 + (0.27521612041849561426E-4 + (-0.57132774537670953638E-5 + (-0.48404772799207914899E-7 + (0.35268354132474570493E-8 + (-0.32383477652514618094E-11 + (-0.19334202915190442501E-11 + 0.32333189861286460270E-13*t)*t)*t)*t)*t)*t)*t)*t;
252 }
253 case 36: {
254 const T t = static_cast<T> (2)*y100 - static_cast<T> (73);
255 return 0.40858275583808707870 + (0.11006378016848466550E-1 + (-0.76396376685213286033E-5 + (-0.59609835484245791439E-5 + (-0.13834610033859313213E-7 + (0.33406952974861448790E-8 + (-0.26474915974296612559E-10 + (-0.13750229270354351983E-11 + 0.36169366979417390637E-13*t)*t)*t)*t)*t)*t)*t)*t;
256 }
257 case 37: {
258 const T t = static_cast<T> (2)*y100 - static_cast<T> (75);
259 return 0.43051714914006682977 + (0.10904106549500816155E-1 + (-0.43477527256787216909E-4 + (-0.59429739547798343948E-5 + (0.17639200194091885949E-7 + (0.29235991689639918688E-8 + (-0.41718791216277812879E-10 + (-0.81023337739508049606E-12 + 0.33618915934461994428E-13*t)*t)*t)*t)*t)*t)*t)*t;
260 }
261 case 38: {
262 const T t = static_cast<T> (2)*y100 - static_cast<T> (77);
263 return 0.45210428135559607406 + (0.10659670756384400554E-1 + (-0.78488639913256978087E-4 + (-0.56919860886214735936E-5 + (0.44181850467477733407E-7 + (0.23694306174312688151E-8 + (-0.49492621596685443247E-10 + (-0.31827275712126287222E-12 + 0.27494438742721623654E-13*t)*t)*t)*t)*t)*t)*t)*t;
264 }
265 case 39: {
266 const T t = static_cast<T> (2)*y100 - static_cast<T> (79);
267 return 0.47306491195005224077 + (0.10279006119745977570E-1 + (-0.11140268171830478306E-3 + (-0.52518035247451432069E-5 + (0.64846898158889479518E-7 + (0.17603624837787337662E-8 + (-0.51129481592926104316E-10 + (0.62674584974141049511E-13 + 0.20055478560829935356E-13*t)*t)*t)*t)*t)*t)*t)*t;
268 }
269 case 40: {
270 const T t = static_cast<T> (2)*y100 - static_cast<T> (81);
271 return 0.49313638965719857647 + (0.97725799114772017662E-2 + (-0.14122854267291533334E-3 + (-0.46707252568834951907E-5 + (0.79421347979319449524E-7 + (0.11603027184324708643E-8 + (-0.48269605844397175946E-10 + (0.32477251431748571219E-12 + 0.12831052634143527985E-13*t)*t)*t)*t)*t)*t)*t)*t;
272 }
273 case 41: {
274 const T t = static_cast<T> (2)*y100 - static_cast<T> (83);
275 return 0.51208057433416004042 + (0.91542422354009224951E-2 + (-0.16726530230228647275E-3 + (-0.39964621752527649409E-5 + (0.88232252903213171454E-7 + (0.61343113364949928501E-9 + (-0.42516755603130443051E-10 + (0.47910437172240209262E-12 + 0.66784341874437478953E-14*t)*t)*t)*t)*t)*t)*t)*t;
276 }
277 case 42: {
278 const T t = static_cast<T> (2)*y100 - static_cast<T> (85);
279 return 0.52968945458607484524 + (0.84400880445116786088E-2 + (-0.18908729783854258774E-3 + (-0.32725905467782951931E-5 + (0.91956190588652090659E-7 + (0.14593989152420122909E-9 + (-0.35239490687644444445E-10 + 0.54613829888448694898E-12*t)*t)*t)*t)*t)*t)*t;
280 }
281 case 43: {
282 const T t = static_cast<T> (2)*y100 - static_cast<T> (87);
283 return 0.54578857454330070965 + (0.76474155195880295311E-2 + (-0.20651230590808213884E-3 + (-0.25364339140543131706E-5 + (0.91455367999510681979E-7 + (-0.23061359005297528898E-9 + (-0.27512928625244444444E-10 + 0.54895806008493285579E-12*t)*t)*t)*t)*t)*t)*t;
284 }
285 case 44: {
286 const T t = static_cast<T> (2)*y100 - static_cast<T> (89);
287 return 0.56023851910298493910 + (0.67938321739997196804E-2 + (-0.21956066613331411760E-3 + (-0.18181127670443266395E-5 + (0.87650335075416845987E-7 + (-0.51548062050366615977E-9 + (-0.20068462174044444444E-10 + 0.50912654909758187264E-12*t)*t)*t)*t)*t)*t)*t;
288 }
289 case 45: {
290 const T t = static_cast<T> (2)*y100 - static_cast<T> (91);
291 return 0.57293478057455721150 + (0.58965321010394044087E-2 + (-0.22841145229276575597E-3 + (-0.11404605562013443659E-5 + (0.81430290992322326296E-7 + (-0.71512447242755357629E-9 + (-0.13372664928000000000E-10 + 0.44461498336689298148E-12*t)*t)*t)*t)*t)*t)*t;
292 }
293 case 46: {
294 const T t = static_cast<T> (2)*y100 - static_cast<T> (93);
295 return 0.58380635448407827360 + (0.49717469530842831182E-2 + (-0.23336001540009645365E-3 + (-0.51952064448608850822E-6 + (0.73596577815411080511E-7 + (-0.84020916763091566035E-9 + (-0.76700972702222222221E-11 + 0.36914462807972467044E-12*t)*t)*t)*t)*t)*t)*t;
296 }
297 case 47: {
298 const T t = static_cast<T> (2)*y100 - static_cast<T> (95);
299 return 0.59281340237769489597 + (0.40343592069379730568E-2 + (-0.23477963738658326185E-3 + (0.34615944987790224234E-7 + (0.64832803248395814574E-7 + (-0.90329163587627007971E-9 + (-0.30421940400000000000E-11 + 0.29237386653743536669E-12*t)*t)*t)*t)*t)*t)*t;
300 }
301 case 48: {
302 const T t = static_cast<T> (2)*y100 - static_cast<T> (97);
303 return 0.59994428743114271918 + (0.30976579788271744329E-2 + (-0.23308875765700082835E-3 + (0.51681681023846925160E-6 + (0.55694594264948268169E-7 + (-0.91719117313243464652E-9 + (0.53982743680000000000E-12 + 0.22050829296187771142E-12*t)*t)*t)*t)*t)*t)*t;
304 }
305 case 49: {
306 const T t = static_cast<T> (2)*y100 - static_cast<T> (99);
307 return 0.60521224471819875444 + (0.21732138012345456060E-2 + (-0.22872428969625997456E-3 + (0.92588959922653404233E-6 + (0.46612665806531930684E-7 + (-0.89393722514414153351E-9 + (0.31718550353777777778E-11 + 0.15705458816080549117E-12*t)*t)*t)*t)*t)*t)*t;
308 }
309 case 50: {
310 const T t = static_cast<T> (2)*y100 - static_cast<T> (101);
311 return 0.60865189969791123620 + (0.12708480848877451719E-2 + (-0.22212090111534847166E-3 + (0.12636236031532793467E-5 + (0.37904037100232937574E-7 + (-0.84417089968101223519E-9 + (0.49843180828444444445E-11 + 0.10355439441049048273E-12*t)*t)*t)*t)*t)*t)*t;
312 }
313 case 51: {
314 const T t = static_cast<T> (2)*y100 - static_cast<T> (103);
315 return 0.61031580103499200191 + (0.39867436055861038223E-3 + (-0.21369573439579869291E-3 + (0.15339402129026183670E-5 + (0.29787479206646594442E-7 + (-0.77687792914228632974E-9 + (0.61192452741333333334E-11 + 0.60216691829459295780E-13*t)*t)*t)*t)*t)*t)*t;
316 }
317 case 52: {
318 const T t = static_cast<T> (2)*y100 - static_cast<T> (105);
319 return 0.61027109047879835868 + (-0.43680904508059878254E-3 + (-0.20383783788303894442E-3 + (0.17421743090883439959E-5 + (0.22400425572175715576E-7 + (-0.69934719320045128997E-9 + (0.67152759655111111110E-11 + 0.26419960042578359995E-13*t)*t)*t)*t)*t)*t)*t;
320 }
321 case 53: {
322 const T t = static_cast<T> (2)*y100 - static_cast<T> (107);
323 return 0.60859639489217430521 + (-0.12305921390962936873E-2 + (-0.19290150253894682629E-3 + (0.18944904654478310128E-5 + (0.15815530398618149110E-7 + (-0.61726850580964876070E-9 + 0.68987888999111111110E-11*t)*t)*t)*t)*t)*t;
324 }
325 case 54: {
326 const T t = static_cast<T> (2)*y100 - static_cast<T> (109);
327 return 0.60537899426486075181 + (-0.19790062241395705751E-2 + (-0.18120271393047062253E-3 + (0.19974264162313241405E-5 + (0.10055795094298172492E-7 + (-0.53491997919318263593E-9 + (0.67794550295111111110E-11 - 0.17059208095741511603E-13*t)*t)*t)*t)*t)*t)*t;
328 }
329 case 55: {
330 const T t = static_cast<T> (2)*y100 - static_cast<T> (111);
331 return 0.60071229457904110537 + (-0.26795676776166354354E-2 + (-0.16901799553627508781E-3 + (0.20575498324332621581E-5 + (0.51077165074461745053E-8 + (-0.45536079828057221858E-9 + (0.64488005516444444445E-11 - 0.29311677573152766338E-13*t)*t)*t)*t)*t)*t)*t;
332 }
333 case 56: {
334 const T t = static_cast<T> (2)*y100 - static_cast<T> (113);
335 return 0.59469361520112714738 + (-0.33308208190600993470E-2 + (-0.15658501295912405679E-3 + (0.20812116912895417272E-5 + (0.93227468760614182021E-9 + (-0.38066673740116080415E-9 + (0.59806790359111111110E-11 - 0.36887077278950440597E-13*t)*t)*t)*t)*t)*t)*t;
336 }
337 case 57: {
338 const T t = static_cast<T> (2)*y100 - static_cast<T> (115);
339 return 0.58742228631775388268 + (-0.39321858196059227251E-2 + (-0.14410441141450122535E-3 + (0.20743790018404020716E-5 + (-0.25261903811221913762E-8 + (-0.31212416519526924318E-9 + (0.54328422462222222221E-11 - 0.40864152484979815972E-13*t)*t)*t)*t)*t)*t)*t;
340 }
341 case 58: {
342 const T t = static_cast<T> (2)*y100 - static_cast<T> (117);
343 return 0.57899804200033018447 + (-0.44838157005618913447E-2 + (-0.13174245966501437965E-3 + (0.20425306888294362674E-5 + (-0.53330296023875447782E-8 + (-0.25041289435539821014E-9 + (0.48490437205333333334E-11 - 0.42162206939169045177E-13*t)*t)*t)*t)*t)*t)*t;
344 }
345 case 59: {
346 const T t = static_cast<T> (2)*y100 - static_cast<T> (119);
347 return 0.56951968796931245974 + (-0.49864649488074868952E-2 + (-0.11963416583477567125E-3 + (0.19906021780991036425E-5 + (-0.75580140299436494248E-8 + (-0.19576060961919820491E-9 + (0.42613011928888888890E-11 - 0.41539443304115604377E-13*t)*t)*t)*t)*t)*t)*t;
348 }
349 case 60: {
350 const T t = static_cast<T> (2)*y100 - static_cast<T> (121);
351 return 0.55908401930063918964 + (-0.54413711036826877753E-2 + (-0.10788661102511914628E-3 + (0.19229663322982839331E-5 + (-0.92714731195118129616E-8 + (-0.14807038677197394186E-9 + (0.36920870298666666666E-11 - 0.39603726688419162617E-13*t)*t)*t)*t)*t)*t)*t;
352 }
353 case 61: {
354 const T t = static_cast<T> (2)*y100 - static_cast<T> (123);
355 return 0.54778496152925675315 + (-0.58501497933213396670E-2 + (-0.96582314317855227421E-4 + (0.18434405235069270228E-5 + (-0.10541580254317078711E-7 + (-0.10702303407788943498E-9 + (0.31563175582222222222E-11 - 0.36829748079110481422E-13*t)*t)*t)*t)*t)*t)*t;
356 }
357 case 62: {
358 const T t = static_cast<T> (2)*y100 - static_cast<T> (125);
359 return 0.53571290831682823999 + (-0.62147030670760791791E-2 + (-0.85782497917111760790E-4 + (0.17553116363443470478E-5 + (-0.11432547349815541084E-7 + (-0.72157091369041330520E-10 + (0.26630811607111111111E-11 - 0.33578660425893164084E-13*t)*t)*t)*t)*t)*t)*t;
360 }
361 case 63: {
362 const T t = static_cast<T> (2)*y100 - static_cast<T> (127);
363 return 0.52295422962048434978 + (-0.65371404367776320720E-2 + (-0.75530164941473343780E-4 + (0.16613725797181276790E-5 + (-0.12003521296598910761E-7 + (-0.42929753689181106171E-10 + (0.22170894940444444444E-11 - 0.30117697501065110505E-13*t)*t)*t)*t)*t)*t)*t;
364 }
365 case 64: {
366 const T t = static_cast<T> (2)*y100 - static_cast<T> (129);
367 return 0.50959092577577886140 + (-0.68197117603118591766E-2 + (-0.65852936198953623307E-4 + (0.15639654113906716939E-5 + (-0.12308007991056524902E-7 + (-0.18761997536910939570E-10 + (0.18198628922666666667E-11 - 0.26638355362285200932E-13*t)*t)*t)*t)*t)*t)*t;
368 }
369 case 65: {
370 const T t = static_cast<T> (2)*y100 - static_cast<T> (131);
371 return 0.49570040481823167970 + (-0.70647509397614398066E-2 + (-0.56765617728962588218E-4 + (0.14650274449141448497E-5 + (-0.12393681471984051132E-7 + (0.92904351801168955424E-12 + (0.14706755960177777778E-11 - 0.23272455351266325318E-13*t)*t)*t)*t)*t)*t)*t;
372 }
373 case 66: {
374 const T t = static_cast<T> (2)*y100 - static_cast<T> (133);
375 return 0.48135536250935238066 + (-0.72746293327402359783E-2 + (-0.48272489495730030780E-4 + (0.13661377309113939689E-5 + (-0.12302464447599382189E-7 + (0.16707760028737074907E-10 + (0.11672928324444444444E-11 - 0.20105801424709924499E-13*t)*t)*t)*t)*t)*t)*t;
376 }
377 case 67: {
378 const T t = static_cast<T> (2)*y100 - static_cast<T> (135);
379 return 0.46662374675511439448 + (-0.74517177649528487002E-2 + (-0.40369318744279128718E-4 + (0.12685621118898535407E-5 + (-0.12070791463315156250E-7 + (0.29105507892605823871E-10 + (0.90653314645333333334E-12 - 0.17189503312102982646E-13*t)*t)*t)*t)*t)*t)*t;
380 }
381 case 68: {
382 const T t = static_cast<T> (2)*y100 - static_cast<T> (137);
383 return 0.45156879030168268778 + (-0.75983560650033817497E-2 + (-0.33045110380705139759E-4 + (0.11732956732035040896E-5 + (-0.11729986947158201869E-7 + (0.38611905704166441308E-10 + (0.68468768305777777779E-12 - 0.14549134330396754575E-13*t)*t)*t)*t)*t)*t)*t;
384 }
385 case 69: {
386 const T t = static_cast<T> (2)*y100 - static_cast<T> (139);
387 return 0.43624909769330896904 + (-0.77168291040309554679E-2 + (-0.26283612321339907756E-4 + (0.10811018836893550820E-5 + (-0.11306707563739851552E-7 + (0.45670446788529607380E-10 + (0.49782492549333333334E-12 - 0.12191983967561779442E-13*t)*t)*t)*t)*t)*t)*t;
388 }
389 case 70: {
390 const T t = static_cast<T> (2)*y100 - static_cast<T> (141);
391 return 0.42071877443548481181 + (-0.78093484015052730097E-2 + (-0.20064596897224934705E-4 + (0.99254806680671890766E-6 + (-0.10823412088884741451E-7 + (0.50677203326904716247E-10 + (0.34200547594666666666E-12 - 0.10112698698356194618E-13*t)*t)*t)*t)*t)*t)*t;
392 }
393 case 71: {
394 const T t = static_cast<T> (2)*y100 - static_cast<T> (143);
395 return 0.40502758809710844280 + (-0.78780384460872937555E-2 + (-0.14364940764532853112E-4 + (0.90803709228265217384E-6 + (-0.10298832847014466907E-7 + (0.53981671221969478551E-10 + (0.21342751381333333333E-12 - 0.82975901848387729274E-14*t)*t)*t)*t)*t)*t)*t;
396 }
397 case 72: {
398 const T t = static_cast<T> (2)*y100 - static_cast<T> (145);
399 return 0.38922115269731446690 + (-0.79249269708242064120E-2 + (-0.91595258799106970453E-5 + (0.82783535102217576495E-6 + (-0.97484311059617744437E-8 + (0.55889029041660225629E-10 + (0.10851981336888888889E-12 - 0.67278553237853459757E-14*t)*t)*t)*t)*t)*t)*t;
400 }
401 case 73: {
402 const T t = static_cast<T> (2)*y100 - static_cast<T> (147);
403 return 0.37334112915460307335 + (-0.79519385109223148791E-2 + (-0.44219833548840469752E-5 + (0.75209719038240314732E-6 + (-0.91848251458553190451E-8 + (0.56663266668051433844E-10 + (0.23995894257777777778E-13 - 0.53819475285389344313E-14*t)*t)*t)*t)*t)*t)*t;
404 }
405 case 74: {
406 const T t = static_cast<T> (2)*y100 - static_cast<T> (149);
407 return 0.35742543583374223085 + (-0.79608906571527956177E-2 + (-0.12530071050975781198E-6 + (0.68088605744900552505E-6 + (-0.86181844090844164075E-8 + (0.56530784203816176153E-10 + (-0.43120012248888888890E-13 - 0.42372603392496813810E-14*t)*t)*t)*t)*t)*t)*t;
408 }
409 case 75: {
410 const T t = static_cast<T> (2)*y100 - static_cast<T> (151);
411 return 0.34150846431979618536 + (-0.79534924968773806029E-2 + (0.37576885610891515813E-5 + (0.61419263633090524326E-6 + (-0.80565865409945960125E-8 + (0.55684175248749269411E-10 + (-0.95486860764444444445E-13 - 0.32712946432984510595E-14*t)*t)*t)*t)*t)*t)*t;
412 }
413 case 76: {
414 const T t = static_cast<T> (2)*y100 - static_cast<T> (153);
415 return 0.32562129649136346824 + (-0.79313448067948884309E-2 + (0.72539159933545300034E-5 + (0.55195028297415503083E-6 + (-0.75063365335570475258E-8 + (0.54281686749699595941E-10 - 0.13545424295111111111E-12*t)*t)*t)*t)*t)*t;
416 }
417 case 77: {
418 const T t = static_cast<T> (2)*y100 - static_cast<T> (155);
419 return 0.30979191977078391864 + (-0.78959416264207333695E-2 + (0.10389774377677210794E-4 + (0.49404804463196316464E-6 + (-0.69722488229411164685E-8 + (0.52469254655951393842E-10 - 0.16507860650666666667E-12*t)*t)*t)*t)*t)*t;
420 }
421 case 78: {
422 const T t = static_cast<T> (2)*y100 - static_cast<T> (157);
423 return 0.29404543811214459904 + (-0.78486728990364155356E-2 + (0.13190885683106990459E-4 + (0.44034158861387909694E-6 + (-0.64578942561562616481E-8 + (0.50354306498006928984E-10 - 0.18614473550222222222E-12*t)*t)*t)*t)*t)*t;
424 }
425 case 79: {
426 const T t = static_cast<T> (2)*y100 - static_cast<T> (159);
427 return 0.27840427686253660515 + (-0.77908279176252742013E-2 + (0.15681928798708548349E-4 + (0.39066226205099807573E-6 + (-0.59658144820660420814E-8 + (0.48030086420373141763E-10 - 0.20018995173333333333E-12*t)*t)*t)*t)*t)*t;
428 }
429 case 80: {
430 const T t = static_cast<T> (2)*y100 - static_cast<T> (161);
431 return 0.26288838011163800908 + (-0.77235993576119469018E-2 + (0.17886516796198660969E-4 + (0.34482457073472497720E-6 + (-0.54977066551955420066E-8 + (0.45572749379147269213E-10 - 0.20852924954666666667E-12*t)*t)*t)*t)*t)*t;
432 }
433 case 81: {
434 const T t = static_cast<T> (2)*y100 - static_cast<T> (163);
435 return 0.24751539954181029717 + (-0.76480877165290370975E-2 + (0.19827114835033977049E-4 + (0.30263228619976332110E-6 + (-0.50545814570120129947E-8 + (0.43043879374212005966E-10 - 0.21228012028444444444E-12*t)*t)*t)*t)*t)*t;
436 }
437 case 82: {
438 const T t = static_cast<T> (2)*y100 - static_cast<T> (165);
439 return 0.23230087411688914593 + (-0.75653060136384041587E-2 + (0.21524991113020016415E-4 + (0.26388338542539382413E-6 + (-0.46368974069671446622E-8 + (0.40492715758206515307E-10 - 0.21238627815111111111E-12*t)*t)*t)*t)*t)*t;
440 }
441 case 83: {
442 const T t = static_cast<T> (2)*y100 - static_cast<T> (167);
443 return 0.21725840021297341931 + (-0.74761846305979730439E-2 + (0.23000194404129495243E-4 + (0.22837400135642906796E-6 + (-0.42446743058417541277E-8 + (0.37958104071765923728E-10 - 0.20963978568888888889E-12*t)*t)*t)*t)*t)*t;
444 }
445 case 84: {
446 const T t = static_cast<T> (2)*y100 - static_cast<T> (169);
447 return 0.20239979200788191491 + (-0.73815761980493466516E-2 + (0.24271552727631854013E-4 + (0.19590154043390012843E-6 + (-0.38775884642456551753E-8 + (0.35470192372162901168E-10 - 0.20470131678222222222E-12*t)*t)*t)*t)*t)*t;
448 }
449 case 85: {
450 const T t = static_cast<T> (2)*y100 - static_cast<T> (171);
451 return 0.18773523211558098962 + (-0.72822604530339834448E-2 + (0.25356688567841293697E-4 + (0.16626710297744290016E-6 + (-0.35350521468015310830E-8 + (0.33051896213898864306E-10 - 0.19811844544000000000E-12*t)*t)*t)*t)*t)*t;
452 }
453 case 86: {
454 const T t = static_cast<T> (2)*y100 - static_cast<T> (173);
455 return 0.17327341258479649442 + (-0.71789490089142761950E-2 + (0.26272046822383820476E-4 + (0.13927732375657362345E-6 + (-0.32162794266956859603E-8 + (0.30720156036105652035E-10 - 0.19034196304000000000E-12*t)*t)*t)*t)*t)*t;
456 }
457 case 87: {
458 const T t = static_cast<T> (2)*y100 - static_cast<T> (175);
459 return 0.15902166648328672043 + (-0.70722899934245504034E-2 + (0.27032932310132226025E-4 + (0.11474573347816568279E-6 + (-0.29203404091754665063E-8 + (0.28487010262547971859E-10 - 0.18174029063111111111E-12*t)*t)*t)*t)*t)*t;
460 }
461 case 88: {
462 const T t = static_cast<T> (2)*y100 - static_cast<T> (177);
463 return 0.14498609036610283865 + (-0.69628725220045029273E-2 + (0.27653554229160596221E-4 + (0.92493727167393036470E-7 + (-0.26462055548683583849E-8 + (0.26360506250989943739E-10 - 0.17261211260444444444E-12*t)*t)*t)*t)*t)*t;
464 }
465 case 89: {
466 const T t = static_cast<T> (2)*y100 - static_cast<T> (179);
467 return 0.13117165798208050667 + (-0.68512309830281084723E-2 + (0.28147075431133863774E-4 + (0.72351212437979583441E-7 + (-0.23927816200314358570E-8 + (0.24345469651209833155E-10 - 0.16319736960000000000E-12*t)*t)*t)*t)*t)*t;
468 }
469 case 90: {
470 const T t = static_cast<T> (2)*y100 - static_cast<T> (181);
471 return static_cast<T> (0.11758232561160626306) + (-0.67378491192463392927E-2 + (0.28525664781722907847E-4 + (0.54156999310046790024E-7 + (-0.21589405340123827823E-8 + (0.22444150951727334619E-10 - 0.15368675584000000000E-12*t)*t)*t)*t)*t)*t;
472 }
473 case 91: {
474 const T t = static_cast<T> (2)*y100 - static_cast<T> (183);
475 return static_cast<T> (0.10422112945361673560) + (-0.66231638959845581564E-2 + (0.28800551216363918088E-4 + (0.37758983397952149613E-7 + (-0.19435423557038933431E-8 + (0.20656766125421362458E-10 - 0.14422990012444444444E-12*t)*t)*t)*t)*t)*t;
476 }
477 case 92: {
478 const T t = static_cast<T> (2)*y100 - static_cast<T> (185);
479 return static_cast<T> (0.91090275493541084785E-1) + (-0.65075691516115160062E-2 + (0.28982078385527224867E-4 + (0.23014165807643012781E-7 + (-0.17454532910249875958E-8 + (0.18981946442680092373E-10 - 0.13494234691555555556E-12*t)*t)*t)*t)*t)*t;
480 }
481 case 93: {
482 const T t = static_cast<T> (2)*y100 - static_cast<T> (187);
483 return static_cast<T> (0.78191222288771379358E-1) + (-0.63914190297303976434E-2 + (0.29079759021299682675E-4 + (0.97885458059415717014E-8 + (-0.15635596116134296819E-8 + (0.17417110744051331974E-10 - 0.12591151763555555556E-12*t)*t)*t)*t)*t)*t;
484 }
485 case 94: {
486 const T t = static_cast<T> (2)*y100 - static_cast<T> (189);
487 return static_cast<T> (0.65524757106147402224E-1) + (-0.62750311956082444159E-2 + (0.29102328354323449795E-4 + (-0.20430838882727954582E-8 + (-0.13967781903855367270E-8 + (0.15958771833747057569E-10 - 0.11720175765333333333E-12*t)*t)*t)*t)*t)*t;
488 }
489 case 95: {
490 const T t = static_cast<T> (2)*y100 - static_cast<T> (191);
491 return static_cast<T> (0.53091065838453612773E-1) + (-0.61586898417077043662E-2 + (0.29057796072960100710E-4 + (-0.12597414620517987536E-7 + (-0.12440642607426861943E-8 + (0.14602787128447932137E-10 - 0.10885859114666666667E-12*t)*t)*t)*t)*t)*t;
492 }
493 case 96: {
494 const T t = static_cast<T> (2)*y100 - static_cast<T> (193);
495 return static_cast<T> (0.40889797115352738582E-1) + (-0.60426484889413678200E-2 + (0.28953496450191694606E-4 + (-0.21982952021823718400E-7 + (-0.11044169117553026211E-8 + (0.13344562332430552171E-10 - 0.10091231402844444444E-12*t)*t)*t)*t)*t)*t;
496 }
497 case 97:
498 case 98:
499 case 99:
500 case 100: {
501// Use Taylor expansion for small x (|x| <= 0.0309...)
502// (2/sqrt(pi)) * (x - 2/3 x^3 + 4/15 x^5 - 8/105 x^7 + 16/945 x^9)
503 const T x2 = sq(x);
504 return x*(static_cast<T> (1.1283791670955125739) -
505 x2*(static_cast<T> (0.75225277806367504925) -
506 x2*(static_cast<T> (0.30090111122547001970) -
507 x2*(static_cast<T> (0.085971746064420005629) -
508 x2*static_cast<T> (0.016931216931216931217)))));
509 }
510 }
511
512// Since 0 <= y100 < 101, this is only reached if x is Nan, in which case we
513// should return NaN.
514 return numeric_limits<T>::signaling_NaN();
515 }
516
517//------------------------------------------------------------------------------
521//------------------------------------------------------------------------------
522 template<typename T>
523 T isqpi() {
524 return static_cast<T> (1)/sqrt(static_cast<T> (M_PI));
525 }
526
527//------------------------------------------------------------------------------
536//------------------------------------------------------------------------------
537#if __cplusplus >= 202002L
538 template<std::floating_point T>
539#else
540 template<typename T>
541#endif
542 T w_im(T x) {
543// Continued-fraction expansion is faster.
544 if (abs(x) > static_cast<T> (45)) {
545// 1-term expansion, important to avoid overflow.
546 if (abs(x) > static_cast<T> (5e7)) {
547 return isqpi<T> ()/x;
548 }
549
550// 5-term expansion (rely on compiler for CSE), simplified from:
551// ispi / (x-0.5/(x-1/(x-1.5/(x-2/x))))
552 return isqpi<T> ()*(sq(x)*(sq(x) - static_cast<T> (4.5)) + static_cast<T> (2)) /
553 (x*(sq(x)*(sq(x) - static_cast<T> (5)) + static_cast<T> (3.75)));
554 }
555
556 const T prefix = x >= 0 ? static_cast<T> (1) :
557 static_cast<T> (-1);
558 return prefix*w_im_y100(prefix*x);
559 }
560
561//------------------------------------------------------------------------------
605//------------------------------------------------------------------------------
606#if __cplusplus >= 202002L
607 template<std::floating_point T>
608#else
609 template<typename T>
610#endif
611 T erfcx_y100(const T y100) {
612 switch (static_cast<uint8_t> (y100)) {
613 case 0: {
614 const T t = static_cast<T> (2)*y100 - static_cast<T> (1);
615 return static_cast<T> (0.70878032454106438663E-3) + (0.71234091047026302958E-3 + (0.35779077297597742384E-5 + (0.17403143962587937815E-7 + (0.81710660047307788845E-10 + (0.36885022360434957634E-12 + 0.15917038551111111111E-14*t)*t)*t)*t)*t)*t;
616 }
617 case 1: {
618 const T t = static_cast<T> (2)*y100 - static_cast<T> (3);
619 return static_cast<T> (0.21479143208285144230E-2) + (0.72686402367379996033E-3 + (0.36843175430938995552E-5 + (0.18071841272149201685E-7 + (0.85496449296040325555E-10 + (0.38852037518534291510E-12 + 0.16868473576888888889E-14*t)*t)*t)*t)*t)*t;
620 }
621 case 2: {
622 const T t = static_cast<T> (2)*y100 - static_cast<T> (5);
623 return static_cast<T> (0.36165255935630175090E-2) + (0.74182092323555510862E-3 + (0.37948319957528242260E-5 + (0.18771627021793087350E-7 + (0.89484715122415089123E-10 + (0.40935858517772440862E-12 + 0.17872061464888888889E-14*t)*t)*t)*t)*t)*t;
624 }
625 case 3: {
626 const T t = static_cast<T> (2)*y100 - static_cast<T> (7);
627 return static_cast<T> (0.51154983860031979264E-2) + (0.75722840734791660540E-3 + (0.39096425726735703941E-5 + (0.19504168704300468210E-7 + (0.93687503063178993915E-10 + (0.43143925959079664747E-12 + 0.18939926435555555556E-14*t)*t)*t)*t)*t)*t;
628 }
629 case 4: {
630 const T t = static_cast<T> (2)*y100 - static_cast<T> (9);
631 return static_cast<T> (0.66457513172673049824E-2) + (0.77310406054447454920E-3 + (0.40289510589399439385E-5 + (0.20271233238288381092E-7 + (0.98117631321709100264E-10 + (0.45484207406017752971E-12 + 0.20076352213333333333E-14*t)*t)*t)*t)*t)*t;
632 }
633 case 5: {
634 const T t = static_cast<T> (2)*y100 - static_cast<T> (11);
635 return static_cast<T> (0.82082389970241207883E-2) + (0.78946629611881710721E-3 + (0.41529701552622656574E-5 + (0.21074693344544655714E-7 + (0.10278874108587317989E-9 + (0.47965201390613339638E-12 + 0.21285907413333333333E-14*t)*t)*t)*t)*t)*t;
636 }
637 case 6: {
638 const T t = static_cast<T> (2)*y100 - static_cast<T> (13);
639 return static_cast<T> (0.98039537275352193165E-2) + (0.80633440108342840956E-3 + (0.42819241329736982942E-5 + (0.21916534346907168612E-7 + (0.10771535136565470914E-9 + (0.50595972623692822410E-12 + 0.22573462684444444444E-14*t)*t)*t)*t)*t)*t;
640 }
641 case 7: {
642 const T t = static_cast<T> (2)*y100 - static_cast<T> (15);
643 return static_cast<T> (0.11433927298290302370E-1) + (0.82372858383196561209E-3 + (0.44160495311765438816E-5 + (0.22798861426211986056E-7 + (0.11291291745879239736E-9 + (0.53386189365816880454E-12 + 0.23944209546666666667E-14*t)*t)*t)*t)*t)*t;
644 }
645 case 8: {
646 const T t = static_cast<T> (2)*y100 - static_cast<T> (17);
647 return static_cast<T> (0.13099232878814653979E-1) + (0.84167002467906968214E-3 + (0.45555958988457506002E-5 + (0.23723907357214175198E-7 + (0.11839789326602695603E-9 + (0.56346163067550237877E-12 + 0.25403679644444444444E-14*t)*t)*t)*t)*t)*t;
648 }
649 case 9: {
650 const T t = static_cast<T> (2)*y100 - static_cast<T> (19);
651 return 0.14800987015587535621E-1 + (0.86018092946345943214E-3 + (0.47008265848816866105E-5 + (0.24694040760197315333E-7 + (0.12418779768752299093E-9 + (0.59486890370320261949E-12 + 0.26957764568888888889E-14*t)*t)*t)*t)*t)*t;
652 }
653 case 10: {
654 const T t = static_cast<T> (2)*y100 - static_cast<T> (21);
655 return 0.16540351739394069380E-1 + (0.87928458641241463952E-3 + (0.48520195793001753903E-5 + (0.25711774900881709176E-7 + (0.13030128534230822419E-9 + (0.62820097586874779402E-12 + 0.28612737351111111111E-14*t)*t)*t)*t)*t)*t;
656 }
657 case 11: {
658 const T t = static_cast<T> (2)*y100 - static_cast<T> (23);
659 return 0.18318536789842392647E-1 + (0.89900542647891721692E-3 + (0.50094684089553365810E-5 + (0.26779777074218070482E-7 + (0.13675822186304615566E-9 + (0.66358287745352705725E-12 + 0.30375273884444444444E-14*t)*t)*t)*t)*t)*t;
660 }
661 case 12: {
662 const T t = static_cast<T> (2)*y100 - static_cast<T> (25);
663 return 0.20136801964214276775E-1 + (0.91936908737673676012E-3 + (0.51734830914104276820E-5 + (0.27900878609710432673E-7 + (0.14357976402809042257E-9 + (0.70114790311043728387E-12 + 0.32252476000000000000E-14*t)*t)*t)*t)*t)*t;
664 }
665 case 13: {
666 const T t = static_cast<T> (2)*y100 - static_cast<T> (27);
667 return 0.21996459598282740954E-1 + (0.94040248155366777784E-3 + (0.53443911508041164739E-5 + (0.29078085538049374673E-7 + (0.15078844500329731137E-9 + (0.74103813647499204269E-12 + 0.34251892320000000000E-14*t)*t)*t)*t)*t)*t;
668 }
669 case 14: {
670 const T t = static_cast<T> (2)*y100 - static_cast<T> (29);
671 return 0.23898877187226319502E-1 + (0.96213386835900177540E-3 + (0.55225386998049012752E-5 + (0.30314589961047687059E-7 + (0.15840826497296335264E-9 + (0.78340500472414454395E-12 + 0.36381553564444444445E-14*t)*t)*t)*t)*t)*t;
672 }
673 case 15: {
674 const T t = static_cast<T> (2)*y100 - static_cast<T> (31);
675 return 0.25845480155298518485E-1 + (0.98459293067820123389E-3 + (0.57082915920051843672E-5 + (0.31613782169164830118E-7 + (0.16646478745529630813E-9 + (0.82840985928785407942E-12 + 0.38649975768888888890E-14*t)*t)*t)*t)*t)*t;
676 }
677 case 16: {
678 const T t = static_cast<T> (2)*y100 - static_cast<T> (33);
679 return 0.27837754783474696598E-1 + (0.10078108563256892757E-2 + (0.59020366493792212221E-5 + (0.32979263553246520417E-7 + (0.17498524159268458073E-9 + (0.87622459124842525110E-12 + 0.41066206488888888890E-14*t)*t)*t)*t)*t)*t;
680 }
681 case 17: {
682 const T t = static_cast<T> (2)*y100 - static_cast<T> (35);
683 return 0.29877251304899307550E-1 + (0.10318204245057349310E-2 + (0.61041829697162055093E-5 + (0.34414860359542720579E-7 + (0.18399863072934089607E-9 + (0.92703227366365046533E-12 + 0.43639844053333333334E-14*t)*t)*t)*t)*t)*t;
684 }
685 case 18: {
686 const T t = static_cast<T> (2)*y100 - static_cast<T> (37);
687 return 0.31965587178596443475E-1 + (0.10566560976716574401E-2 + (0.63151633192414586770E-5 + (0.35924638339521924242E-7 + (0.19353584758781174038E-9 + (0.98102783859889264382E-12 + 0.46381060817777777779E-14*t)*t)*t)*t)*t)*t;
688 }
689 case 19: {
690 const T t = static_cast<T> (2)*y100 - static_cast<T> (39);
691 return 0.34104450552588334840E-1 + (0.10823541191350532574E-2 + (0.65354356159553934436E-5 + (0.37512918348533521149E-7 + (0.20362979635817883229E-9 + (0.10384187833037282363E-11 + 0.49300625262222222221E-14*t)*t)*t)*t)*t)*t;
692 }
693 case 20: {
694 const T t = static_cast<T> (2)*y100 - static_cast<T> (41);
695 return 0.36295603928292425716E-1 + (0.11089526167995268200E-2 + (0.67654845095518363577E-5 + (0.39184292949913591646E-7 + (0.21431552202133775150E-9 + (0.10994259106646731797E-11 + 0.52409949102222222221E-14*t)*t)*t)*t)*t)*t;
696 }
697 case 21: {
698 const T t = static_cast<T> (2)*y100 - static_cast<T> (43);
699 return 0.38540888038840509795E-1 + (0.11364917134175420009E-2 + (0.70058230641246312003E-5 + (0.40943644083718586939E-7 + (0.22563034723692881631E-9 + (0.11642841011361992885E-11 + 0.55721092871111111110E-14*t)*t)*t)*t)*t)*t;
700 }
701 case 22: {
702 const T t = static_cast<T> (2)*y100 - static_cast<T> (45);
703 return 0.40842225954785960651E-1 + (0.11650136437945673891E-2 + (0.72569945502343006619E-5 + (0.42796161861855042273E-7 + (0.23761401711005024162E-9 + (0.12332431172381557035E-11 + 0.59246802364444444445E-14*t)*t)*t)*t)*t)*t;
704 }
705 case 23: {
706 const T t = static_cast<T> (2)*y100 - static_cast<T> (47);
707 return 0.43201627431540222422E-1 + (0.11945628793917272199E-2 + (0.75195743532849206263E-5 + (0.44747364553960993492E-7 + (0.25030885216472953674E-9 + (0.13065684400300476484E-11 + 0.63000532853333333334E-14*t)*t)*t)*t)*t)*t;
708 }
709 case 24: {
710 const T t = static_cast<T> (2)*y100 - static_cast<T> (49);
711 return 0.45621193513810471438E-1 + (0.12251862608067529503E-2 + (0.77941720055551920319E-5 + (0.46803119830954460212E-7 + (0.26375990983978426273E-9 + (0.13845421370977119765E-11 + 0.66996477404444444445E-14*t)*t)*t)*t)*t)*t;
712 }
713 case 25: {
714 const T t = static_cast<T> (2)*y100 - static_cast<T> (51);
715 return 0.48103121413299865517E-1 + (0.12569331386432195113E-2 + (0.80814333496367673980E-5 + (0.48969667335682018324E-7 + (0.27801515481905748484E-9 + (0.14674637611609884208E-11 + 0.71249589351111111110E-14*t)*t)*t)*t)*t)*t;
716 }
717 case 26: {
718 const T t = static_cast<T> (2)*y100 - static_cast<T> (53);
719 return 0.50649709676983338501E-1 + (0.12898555233099055810E-2 + (0.83820428414568799654E-5 + (0.51253642652551838659E-7 + (0.29312563849675507232E-9 + (0.15556512782814827846E-11 + 0.75775607822222222221E-14*t)*t)*t)*t)*t)*t;
720 }
721 case 27: {
722 const T t = static_cast<T> (2)*y100 - static_cast<T> (55);
723 return 0.53263363664388864181E-1 + (0.13240082443256975769E-2 + (0.86967260015007658418E-5 + (0.53662102750396795566E-7 + (0.30914568786634796807E-9 + (0.16494420240828493176E-11 + 0.80591079644444444445E-14*t)*t)*t)*t)*t)*t;
724 }
725 case 28: {
726 const T t = static_cast<T> (2)*y100 - static_cast<T> (57);
727 return 0.55946601353500013794E-1 + (0.13594491197408190706E-2 + (0.90262520233016380987E-5 + (0.56202552975056695376E-7 + (0.32613310410503135996E-9 + (0.17491936862246367398E-11 + 0.85713381688888888890E-14*t)*t)*t)*t)*t)*t;
728 }
729 case 29: {
730 const T t = static_cast<T> (2)*y100 - static_cast<T> (59);
731 return 0.58702059496154081813E-1 + (0.13962391363223647892E-2 + (0.93714365487312784270E-5 + (0.58882975670265286526E-7 + (0.34414937110591753387E-9 + (0.18552853109751857859E-11 + 0.91160736711111111110E-14*t)*t)*t)*t)*t)*t;
732 }
733 case 30: {
734 const T t = static_cast<T> (2)*y100 - static_cast<T> (61);
735 return 0.61532500145144778048E-1 + (0.14344426411912015247E-2 + (0.97331446201016809696E-5 + (0.61711860507347175097E-7 + (0.36325987418295300221E-9 + (0.19681183310134518232E-11 + 0.96952238400000000000E-14*t)*t)*t)*t)*t)*t;
736 }
737 case 31: {
738 const T t = static_cast<T> (2)*y100 - static_cast<T> (63);
739 return 0.64440817576653297993E-1 + (0.14741275456383131151E-2 + (0.10112293819576437838E-4 + (0.64698236605933246196E-7 + (0.38353412915303665586E-9 + (0.20881176114385120186E-11 + 0.10310784480000000000E-13*t)*t)*t)*t)*t)*t;
740 }
741 case 32: {
742 const T t = static_cast<T> (2)*y100 - static_cast<T> (65);
743 return 0.67430045633130393282E-1 + (0.15153655418916540370E-2 + (0.10509857606888328667E-4 + (0.67851706529363332855E-7 + (0.40504602194811140006E-9 + (0.22157325110542534469E-11 + 0.10964842115555555556E-13*t)*t)*t)*t)*t)*t;
744 }
745 case 33: {
746 const T t = static_cast<T> (2)*y100 - static_cast<T> (67);
747 return 0.70503365513338850709E-1 + (0.15582323336495709827E-2 + (0.10926868866865231089E-4 + (0.71182482239613507542E-7 + (0.42787405890153386710E-9 + (0.23514379522274416437E-11 + 0.11659571751111111111E-13*t)*t)*t)*t)*t)*t;
748 }
749 case 34: {
750 const T t = static_cast<T> (2)*y100 - static_cast<T> (69);
751 return 0.73664114037944596353E-1 + (0.16028078812438820413E-2 + (0.11364423678778207991E-4 + (0.74701423097423182009E-7 + (0.45210162777476488324E-9 + (0.24957355004088569134E-11 + 0.12397238257777777778E-13*t)*t)*t)*t)*t)*t;
752 }
753 case 35: {
754 const T t = static_cast<T> (2)*y100 - static_cast<T> (71);
755 return 0.76915792420819562379E-1 + (0.16491766623447889354E-2 + (0.11823685320041302169E-4 + (0.78420075993781544386E-7 + (0.47781726956916478925E-9 + (0.26491544403815724749E-11 + 0.13180196462222222222E-13*t)*t)*t)*t)*t)*t;
756 }
757 case 36: {
758 const T t = static_cast<T> (2)*y100 - static_cast<T> (73);
759 return 0.80262075578094612819E-1 + (0.16974279491709504117E-2 + (0.12305888517309891674E-4 + (0.82350717698979042290E-7 + (0.50511496109857113929E-9 + (0.28122528497626897696E-11 + 0.14010889635555555556E-13*t)*t)*t)*t)*t)*t;
760 }
761 case 37: {
762 const T t = static_cast<T> (2)*y100 - static_cast<T> (75);
763 return 0.83706822008980357446E-1 + (0.17476561032212656962E-2 + (0.12812343958540763368E-4 + (0.86506399515036435592E-7 + (0.53409440823869467453E-9 + (0.29856186620887555043E-11 + 0.14891851591111111111E-13*t)*t)*t)*t)*t)*t;
764 }
765 case 38: {
766 const T t = static_cast<T> (2)*y100 - static_cast<T> (77);
767 return 0.87254084284461718231E-1 + (0.17999608886001962327E-2 + (0.13344443080089492218E-4 + (0.90900994316429008631E-7 + (0.56486134972616465316E-9 + (0.31698707080033956934E-11 + 0.15825697795555555556E-13*t)*t)*t)*t)*t)*t;
768 }
769 case 39: {
770 const T t = static_cast<T> (2)*y100 - static_cast<T> (79);
771 return 0.90908120182172748487E-1 + (0.18544478050657699758E-2 + (0.13903663143426120077E-4 + (0.95549246062549906177E-7 + (0.59752787125242054315E-9 + (0.33656597366099099413E-11 + 0.16815130613333333333E-13*t)*t)*t)*t)*t)*t;
772 }
773 case 40: {
774 const T t = static_cast<T> (2)*y100 - static_cast<T> (81);
775 return 0.94673404508075481121E-1 + (0.19112284419887303347E-2 + (0.14491572616545004930E-4 + (0.10046682186333613697E-6 + (0.63221272959791000515E-9 + (0.35736693975589130818E-11 + 0.17862931591111111111E-13*t)*t)*t)*t)*t)*t;
776 }
777 case 41: {
778 const T t = static_cast<T> (2)*y100 - static_cast<T> (83);
779 return 0.98554641648004456555E-1 + (0.19704208544725622126E-2 + (0.15109836875625443935E-4 + (0.10567036667675984067E-6 + (0.66904168640019354565E-9 + (0.37946171850824333014E-11 + 0.18971959040000000000E-13*t)*t)*t)*t)*t)*t;
780 }
781 case 42: {
782 const T t = static_cast<T> (2)*y100 - static_cast<T> (85);
783 return 0.10255677889470089531 + (0.20321499629472857418E-2 + (0.15760224242962179564E-4 + (0.11117756071353507391E-6 + (0.70814785110097658502E-9 + (0.40292553276632563925E-11 + 0.20145143075555555556E-13*t)*t)*t)*t)*t)*t;
784 }
785 case 43: {
786 const T t = static_cast<T> (2)*y100 - static_cast<T> (87);
787 return 0.10668502059865093318 + (0.20965479776148731610E-2 + (0.16444612377624983565E-4 + (0.11700717962026152749E-6 + (0.74967203250938418991E-9 + (0.42783716186085922176E-11 + 0.21385479360000000000E-13*t)*t)*t)*t)*t)*t;
788 }
789 case 44: {
790 const T t = static_cast<T> (2)*y100 - static_cast<T> (89);
791 return 0.11094484319386444474 + (0.21637548491908170841E-2 + (0.17164995035719657111E-4 + (0.12317915750735938089E-6 + (0.79376309831499633734E-9 + (0.45427901763106353914E-11 + 0.22696025653333333333E-13*t)*t)*t)*t)*t)*t;
792 }
793 case 45: {
794 const T t = static_cast<T> (2)*y100 - static_cast<T> (91);
795 return 0.11534201115268804714 + (0.22339187474546420375E-2 + (0.17923489217504226813E-4 + (0.12971465288245997681E-6 + (0.84057834180389073587E-9 + (0.48233721206418027227E-11 + 0.24079890062222222222E-13*t)*t)*t)*t)*t)*t;
796 }
797 case 46: {
798 const T t = static_cast<T> (2)*y100 - static_cast<T> (93);
799 return 0.11988259392684094740 + (0.23071965691918689601E-2 + (0.18722342718958935446E-4 + (0.13663611754337957520E-6 + (0.89028385488493287005E-9 + (0.51210161569225846701E-11 + 0.25540227111111111111E-13*t)*t)*t)*t)*t)*t;
800 }
801 case 47: {
802 const T t = static_cast<T> (2)*y100 - static_cast<T> (95);
803 return 0.12457298393509812907 + (0.23837544771809575380E-2 + (0.19563942105711612475E-4 + (0.14396736847739470782E-6 + (0.94305490646459247016E-9 + (0.54366590583134218096E-11 + 0.27080225920000000000E-13*t)*t)*t)*t)*t)*t;
804 }
805 case 48: {
806 const T t = static_cast<T> (2)*y100 - static_cast<T> (97);
807 return 0.12941991566142438816 + (0.24637684719508859484E-2 + (0.20450821127475879816E-4 + (0.15173366280523906622E-6 + (0.99907632506389027739E-9 + (0.57712760311351625221E-11 + 0.28703099555555555556E-13*t)*t)*t)*t)*t)*t;
808 }
809 case 49: {
810 const T t = static_cast<T> (2)*y100 - static_cast<T> (99);
811 return 0.13443048593088696613 + (0.25474249981080823877E-2 + (0.21385669591362915223E-4 + (0.15996177579900443030E-6 + (0.10585428844575134013E-8 + (0.61258809536787882989E-11 + 0.30412080142222222222E-13*t)*t)*t)*t)*t)*t;
812 }
813 case 50: {
814 const T t = static_cast<T> (2)*y100 - static_cast<T> (101);
815 return 0.13961217543434561353 + (0.26349215871051761416E-2 + (0.22371342712572567744E-4 + (0.16868008199296822247E-6 + (0.11216596910444996246E-8 + (0.65015264753090890662E-11 + 0.32210394506666666666E-13*t)*t)*t)*t)*t)*t;
816 }
817 case 51: {
818 const T t = static_cast<T> (2)*y100 - static_cast<T> (103);
819 return 0.14497287157673800690 + (0.27264675383982439814E-2 + (0.23410870961050950197E-4 + (0.17791863939526376477E-6 + (0.11886425714330958106E-8 + (0.68993039665054288034E-11 + 0.34101266222222222221E-13*t)*t)*t)*t)*t)*t;
820 }
821 case 52: {
822 const T t = static_cast<T> (2)*y100 - static_cast<T> (105);
823 return 0.15052089272774618151 + (0.28222846410136238008E-2 + (0.24507470422713397006E-4 + (0.18770927679626136909E-6 + (0.12597184587583370712E-8 + (0.73203433049229821618E-11 + 0.36087889048888888890E-13*t)*t)*t)*t)*t)*t;
824 }
825 case 53: {
826 const T t = static_cast<T> (2)*y100 - static_cast<T> (107);
827 return 0.15626501395774612325 + (0.29226079376196624949E-2 + (0.25664553693768450545E-4 + (0.19808568415654461964E-6 + (0.13351257759815557897E-8 + (0.77658124891046760667E-11 + 0.38173420035555555555E-13*t)*t)*t)*t)*t)*t;
828 }
829 case 54: {
830 const T t = static_cast<T> (2)*y100 - static_cast<T> (109);
831 return 0.16221449434620737567 + (0.30276865332726475672E-2 + (0.26885741326534564336E-4 + (0.20908350604346384143E-6 + (0.14151148144240728728E-8 + (0.82369170665974313027E-11 + 0.40360957457777777779E-13*t)*t)*t)*t)*t)*t;
832 }
833 case 55: {
834 const T t = static_cast<T> (2)*y100 - static_cast<T> (111);
835 return 0.16837910595412130659 + (0.31377844510793082301E-2 + (0.28174873844911175026E-4 + (0.22074043807045782387E-6 + (0.14999481055996090039E-8 + (0.87348993661930809254E-11 + 0.42653528977777777779E-13*t)*t)*t)*t)*t)*t;
836 }
837 case 56: {
838 const T t = static_cast<T> (2)*y100 - static_cast<T> (113);
839 return 0.17476916455659369953 + (0.32531815370903068316E-2 + (0.29536024347344364074E-4 + (0.23309632627767074202E-6 + (0.15899007843582444846E-8 + (0.92610375235427359475E-11 + 0.45054073102222222221E-13*t)*t)*t)*t)*t)*t;
840 }
841 case 57: {
842 const T t = static_cast<T> (2)*y100 - static_cast<T> (115);
843 return 0.18139556223643701364 + (0.33741744168096996041E-2 + (0.30973511714709500836E-4 + (0.24619326937592290996E-6 + (0.16852609412267750744E-8 + (0.98166442942854895573E-11 + 0.47565418097777777779E-13*t)*t)*t)*t)*t)*t;
844 }
845 case 58: {
846 const T t = static_cast<T> (2)*y100 - static_cast<T> (117);
847 return 0.18826980194443664549 + (0.35010775057740317997E-2 + (0.32491914440014267480E-4 + (0.26007572375886319028E-6 + (0.17863299617388376116E-8 + (0.10403065638343878679E-10 + 0.50190265831111111110E-13*t)*t)*t)*t)*t)*t;
848 }
849 case 59: {
850 const T t = static_cast<T> (2)*y100 - static_cast<T> (119);
851 return 0.19540403413693967350 + (0.36342240767211326315E-2 + (0.34096085096200907289E-4 + (0.27479061117017637474E-6 + (0.18934228504790032826E-8 + (0.11021679075323598664E-10 + 0.52931171733333333334E-13*t)*t)*t)*t)*t)*t;
852 }
853 case 60: {
854 const T t = static_cast<T> (2)*y100 - static_cast<T> (121);
855 return 0.20281109560651886959 + (0.37739673859323597060E-2 + (0.35791165457592409054E-4 + (0.29038742889416172404E-6 + (0.20068685374849001770E-8 + (0.11673891799578381999E-10 + 0.55790523093333333334E-13*t)*t)*t)*t)*t)*t;
856 }
857 case 61: {
858 const T t = static_cast<T> (2)*y100 - static_cast<T> (123);
859 return 0.21050455062669334978 + (0.39206818613925652425E-2 + (0.37582602289680101704E-4 + (0.30691836231886877385E-6 + (0.21270101645763677824E-8 + (0.12361138551062899455E-10 + 0.58770520160000000000E-13*t)*t)*t)*t)*t)*t;
860 }
861 case 62: {
862 const T t = static_cast<T> (2)*y100 - static_cast<T> (125);
863 return 0.21849873453703332479 + (0.40747643554689586041E-2 + (0.39476163820986711501E-4 + (0.32443839970139918836E-6 + (0.22542053491518680200E-8 + (0.13084879235290858490E-10 + 0.61873153262222222221E-13*t)*t)*t)*t)*t)*t;
864 }
865 case 63: {
866 const T t = static_cast<T> (2)*y100 - static_cast<T> (127);
867 return 0.22680879990043229327 + (0.42366354648628516935E-2 + (0.41477956909656896779E-4 + (0.34300544894502810002E-6 + (0.23888264229264067658E-8 + (0.13846596292818514601E-10 + 0.65100183751111111110E-13*t)*t)*t)*t)*t)*t;
868 }
869 case 64: {
870 const T t = static_cast<T> (2)*y100 - static_cast<T> (129);
871 return 0.23545076536988703937 + (0.44067409206365170888E-2 + (0.43594444916224700881E-4 + (0.36268045617760415178E-6 + (0.25312606430853202748E-8 + (0.14647791812837903061E-10 + 0.68453122631111111110E-13*t)*t)*t)*t)*t)*t;
872 }
873 case 65: {
874 const T t = static_cast<T> (2)*y100 - static_cast<T> (131);
875 return 0.24444156740777432838 + (0.45855530511605787178E-2 + (0.45832466292683085475E-4 + (0.38352752590033030472E-6 + (0.26819103733055603460E-8 + (0.15489984390884756993E-10 + 0.71933206364444444445E-13*t)*t)*t)*t)*t)*t;
876 }
877 case 66: {
878 const T t = static_cast<T> (2)*y100 - static_cast<T> (133);
879 return 0.25379911500634264643 + (0.47735723208650032167E-2 + (0.48199253896534185372E-4 + (0.40561404245564732314E-6 + (0.28411932320871165585E-8 + (0.16374705736458320149E-10 + 0.75541379822222222221E-13*t)*t)*t)*t)*t)*t;
880 }
881 case 67: {
882 const T t = static_cast<T> (2)*y100 - static_cast<T> (135);
883 return 0.26354234756393613032 + (0.49713289477083781266E-2 + (0.50702455036930367504E-4 + (0.42901079254268185722E-6 + (0.30095422058900481753E-8 + (0.17303497025347342498E-10 + 0.79278273368888888890E-13*t)*t)*t)*t)*t)*t;
884 }
885 case 68: {
886 const T t = static_cast<T> (2)*y100 - static_cast<T> (137);
887 return 0.27369129607732343398 + (0.51793846023052643767E-2 + (0.53350152258326602629E-4 + (0.45379208848865015485E-6 + (0.31874057245814381257E-8 + (0.18277905010245111046E-10 + 0.83144182364444444445E-13*t)*t)*t)*t)*t)*t;
888 }
889 case 69: {
890 const T t = static_cast<T> (2)*y100 - static_cast<T> (139);
891 return 0.28426714781640316172 + (0.53983341916695141966E-2 + (0.56150884865255810638E-4 + (0.48003589196494734238E-6 + (0.33752476967570796349E-8 + (0.19299477888083469086E-10 + 0.87139049137777777779E-13*t)*t)*t)*t)*t)*t;
892 }
893 case 70: {
894 const T t = static_cast<T> (2)*y100 - static_cast<T> (141);
895 return 0.29529231465348519920 + (0.56288077305420795663E-2 + (0.59113671189913307427E-4 + (0.50782393781744840482E-6 + (0.35735475025851713168E-8 + (0.20369760937017070382E-10 + 0.91262442613333333334E-13*t)*t)*t)*t)*t)*t;
896 }
897 case 71: {
898 const T t = static_cast<T> (2)*y100 - static_cast<T> (143);
899 return 0.30679050522528838613 + (0.58714723032745403331E-2 + (0.62248031602197686791E-4 + (0.53724185766200945789E-6 + (0.37827999418960232678E-8 + (0.21490291930444538307E-10 + 0.95513539182222222221E-13*t)*t)*t)*t)*t)*t;
900 }
901 case 72: {
902 const T t = static_cast<T> (2)*y100 - static_cast<T> (145);
903 return static_cast<T> (0.31878680111173319425) + (0.61270341192339103514E-2 + (0.65564012259707640976E-4 + (0.56837930287837738996E-6 + (0.40035151353392378882E-8 + (0.22662596341239294792E-10 + 0.99891109760000000000E-13*t)*t)*t)*t)*t)*t;
904 }
905 case 73: {
906 const T t = static_cast<T> (2)*y100 - static_cast<T> (147);
907 return static_cast<T> (0.33130773722152622027) + (0.63962406646798080903E-2 + (0.69072209592942396666E-4 + (0.60133006661885941812E-6 + (0.42362183765883466691E-8 + (0.23888182347073698382E-10 + 0.10439349811555555556E-12*t)*t)*t)*t)*t)*t;
908 }
909 case 74: {
910 const T t = static_cast<T> (2)*y100 - static_cast<T> (149);
911 return static_cast<T> (0.34438138658041336523) + (0.66798829540414007258E-2 + (0.72783795518603561144E-4 + (0.63619220443228800680E-6 + (0.44814499336514453364E-8 + (0.25168535651285475274E-10 + 0.10901861383111111111E-12*t)*t)*t)*t)*t)*t;
912 }
913 case 75: {
914 const T t = static_cast<T> (2)*y100 - static_cast<T> (151);
915 return static_cast<T> (0.35803744972380175583) + (0.69787978834882685031E-2 + (0.76710543371454822497E-4 + (0.67306815308917386747E-6 + (0.47397647975845228205E-8 + (0.26505114141143050509E-10 + 0.11376390933333333333E-12*t)*t)*t)*t)*t)*t;
916 }
917 case 76: {
918 const T t = static_cast<T> (2)*y100 - static_cast<T> (153);
919 return static_cast<T> (0.37230734890119724188) + (0.72938706896461381003E-2 + (0.80864854542670714092E-4 + (0.71206484718062688779E-6 + (0.50117323769745883805E-8 + (0.27899342394100074165E-10 + 0.11862637614222222222E-12*t)*t)*t)*t)*t)*t;
920 }
921 case 77: {
922 const T t = static_cast<T> (2)*y100 - static_cast<T> (155);
923 return static_cast<T> (0.38722432730555448223) + (0.76260375162549802745E-2 + (0.85259785810004603848E-4 + (0.75329383305171327677E-6 + (0.52979361368388119355E-8 + (0.29352606054164086709E-10 + 0.12360253370666666667E-12*t)*t)*t)*t)*t)*t;
924 }
925 case 78: {
926 const T t = static_cast<T> (2)*y100 - static_cast<T> (157);
927 return static_cast<T> (0.40282355354616940667) + (0.79762880915029728079E-2 + (0.89909077342438246452E-4 + (0.79687137961956194579E-6 + (0.55989731807360403195E-8 + (0.30866246101464869050E-10 + 0.12868841946666666667E-12*t)*t)*t)*t)*t)*t;
928 }
929 case 79: {
930 const T t = static_cast<T> (2)*y100 - static_cast<T> (159);
931 return static_cast<T> (0.41914223158913787649) + (0.83456685186950463538E-2 + (0.94827181359250161335E-4 + (0.84291858561783141014E-6 + (0.59154537751083485684E-8 + (0.32441553034347469291E-10 + 0.13387957943111111111E-12*t)*t)*t)*t)*t)*t;
932 }
933 case 80: {
934 const T t = static_cast<T> (2)*y100 - static_cast<T> (161);
935 return static_cast<T> (0.43621971639463786896) + (0.87352841828289495773E-2 + (0.10002929142066799966E-3 + (0.89156148280219880024E-6 + (0.62480008150788597147E-8 + (0.34079760983458878910E-10 + 0.13917107176888888889E-12*t)*t)*t)*t)*t)*t;
936 }
937 case 81: {
938 const T t = static_cast<T> (2)*y100 - static_cast<T> (163);
939 return static_cast<T> (0.45409763548534330981) + (0.91463027755548240654E-2 + (0.10553137232446167258E-3 + (0.94293113464638623798E-6 + (0.65972492312219959885E-8 + (0.35782041795476563662E-10 + 0.14455745872000000000E-12*t)*t)*t)*t)*t)*t;
940 }
941 case 82: {
942 const T t = static_cast<T> (2)*y100 - static_cast<T> (165);
943 return static_cast<T> (0.47282001668512331468) + (0.95799574408860463394E-2 + (0.11135019058000067469E-3 + (0.99716373005509038080E-6 + (0.69638453369956970347E-8 + (0.37549499088161345850E-10 + 0.15003280712888888889E-12*t)*t)*t)*t)*t)*t;
944 }
945 case 83: {
946 const T t = static_cast<T> (2)*y100 - static_cast<T> (167);
947 return static_cast<T> (0.49243342227179841649) + (0.10037550043909497071E-1 + (0.11750334542845234952E-3 + (0.10544006716188967172E-5 + (0.73484461168242224872E-8 + (0.39383162326435752965E-10 + 0.15559069118222222222E-12*t)*t)*t)*t)*t)*t;
948 }
949 case 84: {
950 const T t = static_cast<T> (2)*y100 - static_cast<T> (169);
951 return static_cast<T> (0.51298708979209258326) + (0.10520454564612427224E-1 + (0.12400930037494996655E-3 + (0.11147886579371265246E-5 + (0.77517184550568711454E-8 + (0.41283980931872622611E-10 + 0.16122419680000000000E-12*t)*t)*t)*t)*t)*t;
952 }
953 case 85: {
954 const T t = static_cast<T> (2)*y100 - static_cast<T> (171);
955 return static_cast<T> (0.53453307979101369843) + (0.11030120618800726938E-1 + (0.13088741519572269581E-3 + (0.11784797595374515432E-5 + (0.81743383063044825400E-8 + (0.43252818449517081051E-10 + 0.16692592640000000000E-12*t)*t)*t)*t)*t)*t;
956 }
957 case 86: {
958 const T t = static_cast<T> (2)*y100 - static_cast<T> (173);
959 return static_cast<T> (0.55712643071169299478) + (0.11568077107929735233E-1 + (0.13815797838036651289E-3 + (0.12456314879260904558E-5 + (0.86169898078969313597E-8 + (0.45290446811539652525E-10 + 0.17268801084444444444E-12*t)*t)*t)*t)*t)*t;
960 }
961 case 87: {
962 const T t = static_cast<T> (2)*y100 - static_cast<T> (175);
963 return static_cast<T> (0.58082532122519320968) + (0.12135935999503877077E-1 + (0.14584223996665838559E-3 + (0.13164068573095710742E-5 + (0.90803643355106020163E-8 + (0.47397540713124619155E-10 + 0.17850211608888888889E-12*t)*t)*t)*t)*t)*t;
964 }
965 case 88: {
966 const T t = static_cast<T> (2)*y100 - static_cast<T> (177);
967 return static_cast<T> (0.60569124025293375554) + (0.12735396239525550361E-1 + (0.15396244472258863344E-3 + (0.13909744385382818253E-5 + (0.95651595032306228245E-8 + (0.49574672127669041550E-10 + 0.18435945564444444444E-12*t)*t)*t)*t)*t)*t;
968 }
969 case 89: {
970 const T t = static_cast<T> (2)*y100 - static_cast<T> (179);
971 return static_cast<T> (0.63178916494715716894) + (0.13368247798287030927E-1 + (0.16254186562762076141E-3 + (0.14695084048334056083E-5 + (0.10072078109604152350E-7 + (0.51822304995680707483E-10 + 0.19025081422222222222E-12*t)*t)*t)*t)*t)*t;
972 }
973 case 90: {
974 const T t = static_cast<T> (2)*y100 - static_cast<T> (181);
975 return static_cast<T> (0.65918774689725319200) + (0.14036375850601992063E-1 + (0.17160483760259706354E-3 + (0.15521885688723188371E-5 + (0.10601827031535280590E-7 + (0.54140790105837520499E-10 + 0.19616655146666666667E-12*t)*t)*t)*t)*t)*t;
976 }
977 case 91: {
978 const T t = static_cast<T> (2)*y100 - static_cast<T> (183);
979 return static_cast<T> (0.68795950683174433822) + (0.14741765091365869084E-1 + (0.18117679143520433835E-3 + (0.16392004108230585213E-5 + (0.11155116068018043001E-7 + (0.56530360194925690374E-10 + 0.20209663662222222222E-12*t)*t)*t)*t)*t)*t;
980 }
981 case 92: {
982 const T t = static_cast<T> (2)*y100 - static_cast<T> (185);
983 return static_cast<T> (0.71818103808729967036) + (0.15486504187117112279E-1 + (0.19128428784550923217E-3 + (0.17307350969359975848E-5 + (0.11732656736113607751E-7 + (0.58991125287563833603E-10 + 0.20803065333333333333E-12*t)*t)*t)*t)*t)*t;
984 }
985 case 93: {
986 const T t = static_cast<T> (2)*y100 - static_cast<T> (187);
987 return static_cast<T> (0.74993321911726254661) + (0.16272790364044783382E-1 + (0.20195505163377912645E-3 + (0.18269894883203346953E-5 + (0.12335161021630225535E-7 + (0.61523068312169087227E-10 + 0.21395783431111111111E-12*t)*t)*t)*t)*t)*t;
988 }
989 case 94: {
990 const T t = static_cast<T> (2)*y100 - static_cast<T> (189);
991 return static_cast<T> (0.78330143531283492729) + (0.17102934132652429240E-1 + (0.21321800585063327041E-3 + (0.19281661395543913713E-5 + (0.12963340087354341574E-7 + (0.64126040998066348872E-10 + 0.21986708942222222222E-12*t)*t)*t)*t)*t)*t;
992 }
993 case 95: {
994 const T t = static_cast<T> (2)*y100 - static_cast<T> (191);
995 return static_cast<T> (0.81837581041023811832) + (0.17979364149044223802E-1 + (0.22510330592753129006E-3 + (0.20344732868018175389E-5 + (0.13617902941839949718E-7 + (0.66799760083972474642E-10 + 0.22574701262222222222E-12*t)*t)*t)*t)*t)*t;
996 }
997 case 96: {
998 const T t = static_cast<T> (2)*y100 - static_cast<T> (193);
999 return static_cast<T> (0.85525144775685126237) + (0.18904632212547561026E-1 + (0.23764237370371255638E-3 + (0.21461248251306387979E-5 + (0.14299555071870523786E-7 + (0.69543803864694171934E-10 + 0.23158593688888888889E-12*t)*t)*t)*t)*t)*t;
1000 }
1001 case 97: {
1002 const T t = static_cast<T> (2)*y100 - static_cast<T> (195);
1003 return static_cast<T> (0.89402868170849933734) + (0.19881418399127202569E-1 + (0.25086793128395995798E-3 + (0.22633402747585233180E-5 + (0.15008997042116532283E-7 + (0.72357609075043941261E-10 + 0.23737194737777777778E-12*t)*t)*t)*t)*t)*t;
1004 }
1005 case 98: {
1006 const T t = static_cast<T> (2)*y100 - static_cast<T> (197);
1007 return static_cast<T> (0.93481333942870796363) + (0.20912536329780368893E-1 + (0.26481403465998477969E-3 + (0.23863447359754921676E-5 + (0.15746923065472184451E-7 + (0.75240468141720143653E-10 + 0.24309291271111111111E-12*t)*t)*t)*t)*t)*t;
1008 }
1009 case 99: {
1010 const T t = static_cast<T> (2)*y100 - static_cast<T> (199);
1011 return static_cast<T> (0.97771701335885035464) + (static_cast<T> (0.22000938572830479551E-1) + (static_cast<T> (0.27951610702682383001E-3) + (static_cast<T> (0.25153688325245314530E-5) + (static_cast<T> (0.16514019547822821453E-7) + (static_cast<T> (0.78191526829368231251E-10) + static_cast<T> (0.24873652355555555556E-12)*t)*t)*t)*t)*t)*t;
1012 }
1013 }
1014
1015// We only get here if y = 1, i.e. |x| < 4*eps, in which case erfcx is within
1016// 1E-15 of 1..
1017 return static_cast<T> (1);
1018 }
1019
1020//------------------------------------------------------------------------------
1027//------------------------------------------------------------------------------
1028#if __cplusplus >= 202002L
1029 template<std::floating_point T>
1030#else
1031 template<typename T>
1032#endif
1033 T erfcx(T x) {
1034 if (x >= static_cast<T> (0)) {
1035 if (x > static_cast<T> (50)) {
1036// Continued-fraction expansion is faster.
1037 if (x > static_cast<T> (5.0e7)) {
1038// 1-term expansion, important to avoid overflow
1039 return isqpi<T> ()/x;
1040 }
1041
1042// 5-term expansion (rely on compiler for CSE), simplified from:
1043// ispi / (x+0.5/(x+1/(x+1.5/(x+2/x))))
1044 return isqpi<T> ()*(sq(x)*(sq(x) + static_cast<T> (4.5)) + static_cast<T> (2))/(x*(sq(x) * (sq(x) + static_cast<T> (5)) + static_cast<T> (3.75)));
1045 }
1046 return erfcx_y100(static_cast<T> (400)/(static_cast<T> (4) + x));
1047 } else {
1048 return x < static_cast<T> (-26.7) ? numeric_limits<T>::max() :
1049 (x < static_cast<T> (-6.1) ? static_cast<T> (2)*exp(sq(x)) :
1050 static_cast<T> (2)*exp(sq(x)) - erfcx_y100(static_cast<T> (400)/(static_cast<T> (4) - x)));
1051 }
1052 }
1053
1055#if __cplusplus >= 202002L
1056 template<std::floating_point T>
1057#else
1058 template<typename T>
1059#endif
1060 constexpr T expa2n2[] = {
1061 7.64405281671221563E-01,
1062 3.41424527166548425E-01,
1063 8.91072646929412548E-02,
1064 1.35887299055460086E-02,
1065 1.21085455253437481E-03,
1066 6.30452613933449404E-05,
1067 1.91805156577114683E-06,
1068 3.40969447714832381E-08,
1069 3.54175089099469393E-10,
1070 2.14965079583260682E-12,
1071 7.62368911833724354E-15,
1072 1.57982797110681093E-17,
1073 1.91294189103582677E-20,
1074 1.35344656764205340E-23,
1075 5.59535712428588720E-27,
1076 1.35164257972401769E-30,
1077 1.90784582843501167E-34,
1078 1.57351920291442930E-38,
1079 7.58312432328032845E-43,
1080 2.13536275438697082E-47,
1081 3.51352063787195769E-52,
1082 3.37800830266396920E-57,
1083 1.89769439468301000E-62,
1084 6.22929926072668851E-68,
1085 1.19481172006938722E-73,
1086 1.33908181133005953E-79,
1087 8.76924303483223939E-86,
1088 3.35555576166254986E-92,
1089 7.50264110688173024E-99,
1090 9.80192200745410268E-106,
1091 7.48265412822268959E-113,
1092 3.33770122566809425E-120,
1093 8.69934598159861140E-128,
1094 1.32486951484088852E-135,
1095 1.17898144201315253E-143,
1096 6.13039120236180012E-152,
1097 1.86258785950822098E-160,
1098 3.30668408201432783E-169,
1099 3.43017280887946235E-178,
1100 2.07915397775808219E-187,
1101 7.36384545323984966E-197,
1102 1.52394760394085741E-206,
1103 1.84281935046532100E-216,
1104 1.30209553802992923E-226,
1105 5.37588903521080531E-237,
1106 1.29689584599763145E-247,
1107 1.82813078022866562E-258,
1108 1.50576355348684241E-269,
1109 7.24692320799294194E-281,
1110 2.03797051314726829E-292,
1111 3.34880215927873807E-304,
1112 0.0
1113// Underflow (also prevents reads past array end, below)
1114 };
1115
1116//------------------------------------------------------------------------------
1125//------------------------------------------------------------------------------
1126#if __cplusplus >= 202002L
1127 template<std::floating_point T>
1128#else
1129 template<typename T>
1130#endif
1131 T sinh_taylor(T x) {
1132 const T x2 = sq<T> (x);
1133 return x*(static_cast<T> (1) + x2*(static_cast<T> (0.1666666666666666666667) +
1134 static_cast<T> (0.00833333333333333333333)*x2));
1135 }
1136
1137//------------------------------------------------------------------------------
1149//------------------------------------------------------------------------------
1150#if __cplusplus >= 202002L
1151 template<std::floating_point T>
1152#else
1153 template<typename T>
1154#endif
1155 T sinc(T x, T sinx) {
1156 return abs(x) < static_cast<T> (1.0E-4) ? static_cast<T> (1) - static_cast<T> (0.1666666666666666666667)*sq(x) : sinx/x;
1157 }
1158
1159//------------------------------------------------------------------------------
1166//------------------------------------------------------------------------------
1167#if __cplusplus >= 202002L
1168 template<std::floating_point T>
1169#else
1170 template<typename T>
1171#endif
1172 complex_type<T> w(complex_type<T> z) {
1173 if (real(z) == static_cast<T> (0)) {
1174// Give correct sign of 0 in cimag(w)
1175 return complex_type<T> (erfcx(imag(z)), real(z));
1176 } else if (imag(z) == static_cast<T> (0)) {
1177 return complex_type<T> (abs(real(z)) > static_cast<T> (27) ? static_cast<T> (0) : exp(-sq<T> (real(z))),
1178 w_im(real(z)));
1179 }
1180
1181 const T a = static_cast<T> (M_PI)/sqrt(-log(numeric_limits<T>::epsilon()*static_cast<T> (0.5))); // pi / sqrt(-log(eps*0.5))
1182 const T c = static_cast<T> (2)/static_cast<T> (M_PI)*a; // (2/pi) * a;
1183 const T a2 = sq(a);
1184
1185 const T x = abs(real(z));
1186 const T y = imag(z);
1187 const T ya = abs(y);
1188
1189// Return value.
1190 complex_type<T> ret(0.0, 0.0);
1191
1192 T sum1 = static_cast<T> (0);
1193 T sum2 = static_cast<T> (0);
1194 T sum3 = static_cast<T> (0);
1195 T sum4 = static_cast<T> (0);
1196 T sum5 = static_cast<T> (0);
1197
1198// Continued fraction is faster
1199// As pointed out by M. Zaghloul, the continued fraction seems to give a large
1200// relative error in Re w(z) for |x| ~ 6 and small |y|, so use algorithm 816 in
1201// this region:
1202 if (ya > static_cast<T> (7) || (x > static_cast<T> (6) &&
1203 (ya > static_cast<T> (0.1) ||
1204 (x > static_cast<T> (8) && ya > static_cast<T> (1.0E-10)) ||
1205 x > static_cast<T> (28)))) {
1206
1207// Poppe & Wijers suggest using a number of terms
1208//
1209// nu = 3 + 1442 / (26*rho + 77)
1210//
1211// where rho = sqrt((x/x0)^2 + (y/y0)^2) where x0=6.3, y0=4.4.
1212// (They only use this expansion for rho >= 1, but rho a little less than 1
1213// seems okay too.)
1214// Instead, I did my own fit to a slightly different function that avoids the
1215// hypotenuse calculation, using NLopt to minimize the sum of the squares of
1216// the errors in nu with the constraint that the estimated nu be >= minimum nu
1217// to attain machine precision. I also separate the regions where nu == 2 and
1218// nu == 1.
1219 const T xs = y < static_cast<T> (0) ? -real(z) :
1220 real(z); // compute for -z if y < 0
1221 if (x + ya > static_cast<T> (4000)) {
1222// nu <= 2
1223 if (x + ya > static_cast<T> (1e7)) {
1224// nu == 1, w(z) = i/sqrt(pi) / z
1225// scale to avoid overflow
1226 if (x > ya) {
1227 T yax = ya / xs;
1228 T denom = isqpi<T> ()/(xs + yax*ya);
1229 ret = complex_type<T> (denom*yax, denom);
1230 } else if (isinf(ya)) {
1231 return (isnan(x) || y < static_cast<T> (0)) ? complex_type<T> (numeric_limits<T>::signaling_NaN(),
1232 numeric_limits<T>::signaling_NaN()) :
1233 complex_type<T> (0.0, 0.0);
1234 } else {
1235 const T xya = xs/ya;
1236 const T denom = isqpi<T> ()/(xya*xs + ya);
1237 ret = complex_type<T> (denom, denom*xya);
1238 }
1239 } else {
1240// nu == 2, w(z) = i/sqrt(pi) * z / (z*z - 0.5)
1241 const T dr = xs*xs - ya*ya - static_cast<T> (0.5);
1242 const T di = static_cast<T> (2)*xs*ya;
1243 const T denom = isqpi<T> ()/(dr*dr + di*di);
1244 ret = complex_type<T> (denom*(xs*di - ya*dr), denom*(xs*dr + ya*di));
1245 }
1246 } else {
1247// compute nu(z) estimate and do general continued fraction
1248 const T c0 = static_cast<T> (3.9);
1249 const T c1 = static_cast<T> (11.398);
1250 const T c2 = static_cast<T> (0.08254);
1251 const T c3 = static_cast<T> (0.1421);
1252 const T c4 = static_cast<T> (0.2023);
1253 T nu = floor(c0 + c1/(c2*x + c3*ya + c4));
1254 T wr = xs;
1255 T wi = ya;
1256 for (nu = static_cast<T> (0.5)*(nu - static_cast<T> (1)); nu > static_cast<T> (0.4); nu -= static_cast<T> (0.5)) {
1257// w <- z - nu/w:
1258 const T denom = nu/(wr*wr + wi*wi);
1259 wr = xs - wr*denom;
1260 wi = ya + wi*denom;
1261 }
1262// w(z) = i/sqrt(pi) / w:
1263 const T denom2 = isqpi<T> ()/(wr*wr + wi*wi);
1264 ret = complex_type<T> (denom2*wi, denom2*wr);
1265 }
1266 if (y < static_cast<T> (0)) {
1267// use w(z) = 2.0*exp(-z*z) - w(-z),
1268// but be careful of overflow in exp(-z*z)
1269// = exp(-(xs*xs-ya*ya) -2*i*xs*ya)
1270 return static_cast<T> (2)*exp(complex_type<T> ((ya - xs)*(xs + ya),
1271 static_cast<T> (2)*xs*y)) - ret;
1272 } else {
1273 return ret;
1274 }
1275
1276// Note: The test that seems to be suggested in the paper is
1277// x < sqrt(-log(DBL_MIN)), about 26.6, since otherwise exp(-x^2) underflows to
1278// zero and sum1,sum2,sum4 are zero. However, long before this occurs, the
1279// sum1,sum2,sum4 contributions are negligible in double precision; I find that
1280// this happens for x > about 6, for all y. On the other hand, I find that the
1281// case where we compute all of the sums is faster (at least with the
1282// precomputed expa2n2 table) until about x=10. Furthermore, if we try to
1283// compute all of the sums for x > 20, I find that we sometimes run into
1284// numerical problems because underflow/overflow problems start to appear in
1285// the various coefficients of the sums, below. Therefore, we use x < 10 here.
1286 } else if (x < static_cast<T> (10)) {
1287 T prod2ax = static_cast<T> (1);
1288 T prodm2ax = static_cast<T> (1);
1289 T expx2;
1290
1291 if (isnan(y)) {
1292 return complex_type<T> (y, y);
1293 }
1294
1295// Somewhat ugly copy-and-paste duplication here, but I see significant
1296// speedups from using the special-case code with the precomputed exponential,
1297// and the x < 5E-4 special case is needed for accuracy.
1298
1299// Use precomputed exp(-a2*(n*n)) table
1300 if (x < static_cast<T> (5.0E-4)) { // compute sum4 and sum5 together as sum5-sum4
1301 const T x2 = sq(x);
1302// Exp(-x*x) via Taylor
1303 expx2 = static_cast<T> (1) - x2*(static_cast<T> (1) - static_cast<T> (0.5)*x2);
1304// Compute exp(2*a*x) and exp(-2*a*x) via Taylor, to double precision
1305 const T ax2 = static_cast<T> (1.036642960860171859744)*x; // 2*a*x
1306 const T exp2ax = static_cast<T> (1) + ax2*(static_cast<T> (1) + ax2*(static_cast<T> (0.5) + static_cast<T> (0.166666666666666666667)*ax2));
1307 const T expm2ax = static_cast<T> (1) - ax2*(static_cast<T> (1) - ax2*(static_cast<T> (0.5) - static_cast<T> (0.166666666666666666667)*ax2));
1308 for (int n = 1; n < 53; ++n) {
1309 const T coef = expa2n2<T> [n - 1]*expx2/(a2*static_cast<T> (n*n) + sq(y));
1310 prod2ax *= exp2ax;
1311 prodm2ax *= expm2ax;
1312 sum1 += coef;
1313 sum2 += coef*prodm2ax;
1314 sum3 += coef*prod2ax;
1315
1316// really = sum5 - sum4
1317 sum5 += coef*static_cast<T> (2)*a*n*sinh_taylor(static_cast<T> (2)*a*n*x);
1318
1319// test convergence via sum3
1320 if (coef*prod2ax < numeric_limits<T>::epsilon()*sum3) {
1321 break;
1322 }
1323 }
1324 } else {
1325// x > 5E-4, compute sum4 and sum5 separately
1326 expx2 = exp(-sq(x));
1327 const T exp2ax = exp(static_cast<T> (2)*a*x);
1328 const T expm2ax = static_cast<T> (1)/exp2ax;
1329 for (int n = 1; n < 53; ++n) {
1330 const T coef = expa2n2<T> [n - 1]*expx2/(a2*static_cast<T> (n*n) + sq(y));
1331 prod2ax *= exp2ax;
1332 prodm2ax *= expm2ax;
1333 sum1 += coef;
1334 sum2 += coef*prodm2ax;
1335 sum4 += coef*prodm2ax*a*static_cast<T> (n);
1336 sum3 += coef*prod2ax;
1337 sum5 += coef*prod2ax*a*static_cast<T> (n);
1338
1339// test convergence via sum5, since this sum has the slowest decay
1340 if (coef*prod2ax*a*static_cast<T> (n) < numeric_limits<T>::epsilon()*sum5) {
1341 break;
1342 }
1343 }
1344 }
1345
1346// avoid spurious overflow for large negative y
1347// for y < -6, erfcx(y) = 2*exp(y*y) to double precision
1348 const T expx2erfcxy = y > static_cast<T> (-6) ? expx2*erfcx(y) :
1349 static_cast<T> (2)*exp(sq(y) - sq(x));
1350 if (y > static_cast<T> (5)) { // imaginary terms cancel
1351 const T sinxy = sin(x*y);
1352 ret = (expx2erfcxy - c*y*sum1)*cos(static_cast<T> (2)*x*y)
1353 + c*x*expx2*sinxy*sinc(x*y, sinxy);
1354 } else {
1355 const T xs = real(z);
1356 const T sinxy = sin(xs*y);
1357 const T sin2xy = sin(static_cast<T> (2)*xs*y);
1358 const T cos2xy = cos(static_cast<T> (2)*xs*y);
1359 const T coef1 = expx2erfcxy - c*y*sum1;
1360 const T coef2 = c*xs*expx2;
1361 ret = complex_type<T> (coef1*cos2xy + coef2*sinxy*sinc(xs*y, sinxy),
1362 coef2*sinc(static_cast<T> (2)*xs*y, sin2xy) - coef1*sin2xy);
1363 }
1364 } else {
1365// x large: only sum3 & sum5 contribute (see above note)
1366 if (isnan(x)) {
1367 return complex_type<T> (x, x);
1368 }
1369 if (isnan(y)) {
1370 return complex_type<T> (y, y);
1371 }
1372
1373// |y| < 1E-10, so we only need exp(-x*x) term
1374 ret = exp(-sq(x));
1375
1376// (round instead of ceil as in original paper; note that x/a > 1 here)
1377 const T n0 = floor(x/a + static_cast<T> (0.5)); // sum in both directions, starting at n0
1378 const T dx = a*n0 - x;
1379 sum3 = exp(-dx*dx) / (a2*(n0*n0) + sq(y));
1380 sum5 = a*n0 * sum3;
1381 const T exp1 = exp(static_cast<T> (4)*a*dx);
1382 T exp1dn = static_cast<T> (1);
1383 int dn;
1384 for (dn = 1; n0 - static_cast<T> (dn) > static_cast<T> (0); ++dn) {
1385// loop over n0-dn and n0+dn terms
1386 const T np = n0 + static_cast<T> (dn);
1387 const T nm = n0 - static_cast<T> (dn);
1388 T tp = exp(-sq<T> (a*static_cast<T> (dn) + dx));
1389 T tm = tp*(exp1dn *= exp1); // trick to get tm from tp
1390 tp /= (a2*np*np + sq(y));
1391 tm /= (a2*nm*nm + sq(y));
1392 sum3 += tp + tm;
1393 sum5 += a*(np*tp + nm*tm);
1394 if (a*(np*tp + nm*tm) < numeric_limits<T>::epsilon()*sum5) {
1395 return ret + complex_type<T> (static_cast<T> (0.5)*c*y*(sum2 + sum3),
1396 static_cast<T> (0.5)*c*copysign(sum5 - sum4,
1397 real(z)));
1398 }
1399 }
1400 while (true) {
1401// Loop over n0+dn terms only (since n0-dn <= 0)
1402 const T np = n0 + static_cast<T> (dn++);
1403 const T tp = exp(-sq<T> (a*static_cast<T> (dn) + dx))/(a2*(np*np) + sq(y));
1404 sum3 += tp;
1405 sum5 += a*np*tp;
1406 if (a*np*tp < numeric_limits<T>::epsilon()*sum5) {
1407 return ret + complex_type<T> (static_cast<T> (0.5)*c*y*(sum2 + sum3),
1408 static_cast<T> (0.5)*c*copysign(sum5 - sum4,
1409 real(z)));
1410 }
1411 }
1412 }
1413
1414 return ret + complex_type<T> (static_cast<T> (0.5)*c*y*(sum2 + sum3),
1415 static_cast<T> (0.5)*c*copysign(sum5 - sum4,
1416 real(z)));
1417 }
1418
1419//------------------------------------------------------------------------------
1431//------------------------------------------------------------------------------
1432#if __cplusplus >= 202002L
1433 template<std::floating_point T>
1434#else
1435 template<typename T>
1436#endif
1437 complex_type<T> taylor(const complex_type<T> z, const T mRe_z2, const T mIm_z2) {
1438 const complex_type<T> mz2(mRe_z2, mIm_z2); // -z^2
1439 return z*(static_cast<T> (1.1283791670955125739) +
1440 mz2*(static_cast<T> (0.37612638903183752464) +
1441 mz2*(static_cast<T> (0.11283791670955125739) +
1442 mz2*(static_cast<T> (0.026866170645131251760) +
1443 mz2*static_cast<T> (0.0052239776254421878422)))));
1444 }
1445
1446//------------------------------------------------------------------------------
1463//------------------------------------------------------------------------------
1464 template<typename T>
1465 complex_type<T> taylor_erfi(const T x, const T y) {
1466 const T x2 = sq(x);
1467 const T y2 = sq(y);
1468 const T expy2 = exp(y2);
1469 return expy2*complex_type<T> (x*(static_cast<T> (1.1283791670955125739) -
1470 x2*(static_cast<T> (0.37612638903183752464) +
1471 static_cast<T> (0.75225277806367504925)*y2) +
1472 x2*x2*(static_cast<T> (0.11283791670955125739) +
1473 y2*(static_cast<T> (0.45135166683820502956) +
1474 static_cast<T> (0.15045055561273500986)*y2))),
1475 (w_im(y) - x2*y*(static_cast<T> (1.1283791670955125739) -
1476 x2*(static_cast<T> (0.56418958354775628695) +
1477 static_cast<T> (0.37612638903183752464)*y2))));
1478 }
1479
1480//------------------------------------------------------------------------------
1489//------------------------------------------------------------------------------
1490#if __cplusplus >= 202002L
1491 template<std::floating_point T>
1492#else
1493 template<typename T>
1494#endif
1495 complex_type<T> erf_complex(const complex_type<T> z) {
1496 T x = real(z);
1497 T y = imag(z);
1498
1499 if (y == static_cast<T> (0)) {
1500#ifdef CUDA_DEVICE_CODE
1501 return complex_type<T> (erf(x), y);
1502#else
1503 return complex_type<T> (std::erf(x), y);
1504#endif
1505 } else if (x == static_cast<T> (0)) {
1506 const T y2 = sq(y);
1507 return complex_type<T> (x, y2 > static_cast<T> (720) ? (y > static_cast<T> (0) ? numeric_limits<T>::max() :
1508 -numeric_limits<T>::max()) :
1509 exp(y2)*w_im(y));
1510 }
1511
1512// In this case, the denominator of the zeta function is zero so Exp(-zeta^2)
1513// is zero.
1514 if (isinf(x) && isinf(y)) {
1515 return complex_type<T> (static_cast<T> (0.0),
1516 static_cast<T> (0.0));
1517 }
1518
1519 T mRe_z2 = (y - x)*(x + y);
1520 T mIm_z2 = -static_cast<T> (2)*x*y;
1521
1522// Underflow
1523 if (mRe_z2 < -static_cast<T> (750)) {
1524 return (x >= static_cast<T> (0) ? static_cast<T> (1) :
1525 static_cast<T> (-1));
1526 }
1527
1528// Handle positive and negative x via different formulas, using the mirror
1529// symmetries of w, to avoid overflow/underflow problems from multiplying
1530// exponentially large and small quantities.
1531 if (x >= static_cast<T> (0)) {
1532 if (x < static_cast<T> (8.0E-2)) {
1533 if (abs(y) < static_cast<T> (1.0E-2)) {
1534 return taylor(z, mRe_z2, mIm_z2);
1535 } else if (abs(mIm_z2) < static_cast<T> (5.0E-3) &&
1536 x < static_cast<T> (5.0E-3)) {
1537 return taylor_erfi(x, y);
1538 }
1539 }
1540// Don't use complex exp function, since that will produce spurious NaN
1541// values when multiplying w in an overflow situation.
1542 return static_cast<T> (1) -
1543 exp(mRe_z2)*(complex_type<T> (cos(mIm_z2),
1544 sin(mIm_z2))*w(complex_type<T> (-y, x)));
1545 } else {
1546// x < 0
1547 if (x > static_cast<T> (-8.0E-2)) {
1548// Duplicate from above to avoid abs(x) call
1549 if (abs(y) < static_cast<T> (1.0E-2)) {
1550 return taylor(z, mRe_z2, mIm_z2);
1551 } else if (abs(mIm_z2) < 5e-3 && x > -5e-3) {
1552 return taylor_erfi(x, y);
1553 }
1554 } else if (isnan(x)) {
1555 return complex_type<T> (numeric_limits<T>::signaling_NaN(),
1556 y == static_cast<T> (0) ? static_cast<T> (0) :
1557 numeric_limits<T>::signaling_NaN());
1558 }
1559// Don't use complex exp function, since that will produce spurious NaN
1560// values when multiplying w in an overflow situation.
1561 return exp(mRe_z2)*complex_type<T> (cos(mIm_z2),
1562 sin(mIm_z2))*w(complex_type<T> (y, -x)) -
1563 static_cast<T> (1);
1564 }
1565 }
1566
1567//------------------------------------------------------------------------------
1574//------------------------------------------------------------------------------
1575#if __cplusplus >= 202002L
1576 template<std::floating_point T>
1577#else
1578 template<typename T>
1579#endif
1580 complex_type<T> erfi(const complex_type<T> z) {
1581// Avoids NaN instead of doing i<T>*z and -i*temp;
1582 const complex_type<T> temp = erf_complex<T> (complex_type<T> (-imag<T> (z), real<T> (z)));
1583 return complex_type<T> (imag(temp), -real(temp));
1584 }
1585}
1586
1587#endif /* special_functions_h */
special
Name space for special functions.
Definition special_functions.hpp:41
special::w_im_y100
T w_im_y100(T x)
Compute a scaled Dawson integral.
Definition special_functions.hpp:106
special::sinc
T sinc(T x, T sinx)
sinc(x) = sin(x)/x
Definition special_functions.hpp:1155
special::erfcx
complex_type< T > erfcx(complex_type< T > z)
Compute erfcx(z) = exp(z^2)*erfz(z)
Definition special_functions.hpp:65
special::expa2n2
constexpr T expa2n2[]
Precomputed table of expa2n2[n-1] = exp(-a2*n*n) for double-precision a2 = 0.26865....
Definition special_functions.hpp:1060
special::isqpi
T isqpi()
Define 1.0/sqrt(pi)
Definition special_functions.hpp:523
special::w
complex_type< T > w(complex_type< T > z)
Compute w(z) = exp(z^2)*erfz(z).
Definition special_functions.hpp:1172
special::sq
T sq(T x)
x^2
Definition special_functions.hpp:82
special::w_im
T w_im(T x)
Compute specal case of w(z) = exp(z^2)*erfz(z).
Definition special_functions.hpp:542
special::erf_complex
complex_type< T > erf_complex(const complex_type< T > z)
Compute the error function erf(z)
Definition special_functions.hpp:1495
special::erfi
complex_type< T > erfi(const complex_type< T > z)
erfi(z) = -i erf(iz)
Definition special_functions.hpp:1580
special::erfcx_y100
T erfcx_y100(const T y100)
erfcx(x) = exp(x^2)*erfc(x) function, for real x
Definition special_functions.hpp:611
special::taylor_erfi
complex_type< T > taylor_erfi(const T x, const T y)
Taylor erfi series.
Definition special_functions.hpp:1465
special::taylor
complex_type< T > taylor(const complex_type< T > z, const T mRe_z2, const T mIm_z2)
Taylor series.
Definition special_functions.hpp:1437
special::i
constexpr complex_type< T > i(static_cast< T >(0), static_cast< T >(1))
I constant.
special::sinh_taylor
T sinh_taylor(T x)
sinh(x)
Definition special_functions.hpp:1131
complex_type
std::complex< T > complex_type
Type alias for complex types.
Definition special_functions.hpp:36