{

if (std::real(z) == 0.0)

return std::complex<double>(erfcx(std::imag(z)), std::real(z));

else if (std::imag(z) == 0.0)

return std::complex<double>(std::exp(-std::real(z) * std::real(z)), scaled_w_im(std::real(z)));

double a, a2, c;

if (relerr <= DBL_EPSILON)

{

relerr = DBL_EPSILON;

a = 0.518321480430085929872; // pi / sqrt(-log(eps*0.5))

c = 0.329973702884629072537; // (2/pi) * a;

a2 = 0.268657157075235951582; // a^2

}

else

{

if (relerr > 0.1)

relerr = 0.1; // not sensible to compute < 1 digit

a = ModuleBase::PI / std::sqrt(-std::log(relerr * 0.5));

c = (2 / ModuleBase::PI) * a;

a2 = a * a;

}

const double x = std::fabs(std::real(z));

const double y = std::imag(z), ya = std::fabs(y);

std::complex<double> ret = 0.; // return value

double sum1 = 0, sum2 = 0, sum3 = 0, sum4 = 0, sum5 = 0;

if (ya > 7 || (x > 6 && (ya > 0.1 || (x > 8 && ya > 1e-10) || x > 28)))

{

const double ispi = 1 / std::sqrt(ModuleBase::PI); // 1 / sqrt(pi)

double xs = y < 0 ? -std::real(z) : std::real(z); // compute for -z if y < 0

if (x + ya > 4000)

{ // nu <= 2

if (x + ya > 1e7)

{ // nu == 1, w(z) = i/sqrt(pi) / z

// scale to avoid overflow

if (x > ya)

{

double yax = ya / xs;

double denom = ispi / (xs + yax * ya);

ret = std::complex<double>(denom * yax, denom);

}

else if (std::isinf(ya))

return ((std::isnan(x) || y < 0) ? std::complex<double>(NaN, NaN) : std::complex<double>(0, 0));

else

{

double xya = xs / ya;

double denom = ispi / (xya * xs + ya);

ret = std::complex<double>(denom, denom * xya);

}

}

else

{ // nu == 2, w(z) = i/sqrt(pi) * z / (z*z - 0.5)

double dr = xs * xs - ya * ya - 0.5, di = 2 * xs * ya;

double denom = ispi / (dr * dr + di * di);

ret = std::complex<double>(denom * (xs * di - ya * dr), denom * (xs * dr + ya * di));

}

}

else

{ // compute nu(z) estimate and do general continued fraction

const double c0 = 3.9, c1 = 11.398, c2 = 0.08254, c3 = 0.1421, c4 = 0.2023; // fit

double nu = std::floor(c0 + c1 / (c2 * x + c3 * ya + c4));

double wr = xs, wi = ya;

for (nu = 0.5 * (nu - 1); nu > 0.4; nu -= 0.5)

{

// w <- z - nu/w:

double denom = nu / (wr * wr + wi * wi);

wr = xs - wr * denom;

wi = ya + wi * denom;

}

{ // w(z) = i/sqrt(pi) / w:

double denom = ispi / (wr * wr + wi * wi);

ret = std::complex<double>(denom * wi, denom * wr);

}

}

if (y < 0)

{

// use w(z) = 2.0*exp(-z*z) - w(-z),

// but be careful of overflow in exp(-z*z)

// = exp(-(xs*xs-ya*ya) -2*i*xs*ya)

return 2.0 * std::exp(std::complex<double>((ya - xs) * (xs + ya), 2 * xs * y)) - ret;

}

else

return ret;

}

else if (x < 10)

{

double prod2ax = 1, prodm2ax = 1;

double expx2;

if (std::isnan(y))

return std::complex<double>(y, y);

if (relerr == DBL_EPSILON)

{ // use precomputed exp(-a2*(n*n)) table

if (x < 5e-4)

{ // compute sum4 and sum5 together as sum5-sum4

const double x2 = x * x;

expx2 = 1 - x2 * (1 - 0.5 * x2); // exp(-x*x) via Taylor

// compute exp(2*a*x) and exp(-2*a*x) via Taylor, to double precision

const double ax2 = 1.036642960860171859744 * x; // 2*a*x

const double exp2ax = 1 + ax2 * (1 + ax2 * (0.5 + 0.166666666666666666667 * ax2));

const double expm2ax = 1 - ax2 * (1 - ax2 * (0.5 - 0.166666666666666666667 * ax2));

for (int n = 1; 1; ++n)

{

const double coef = expa2n2[n - 1] * expx2 / (a2 * (n * n) + y * y);

prod2ax *= exp2ax;

prodm2ax *= expm2ax;

sum1 += coef;

sum2 += coef * prodm2ax;

sum3 += coef * prod2ax;

// really = sum5 - sum4

sum5 += coef * (2 * a) * n * std::sinh((2 * a) * n * x);

// test convergence via sum3

if (coef * prod2ax < relerr * sum3)

break;

}

}

else

{ // x > 5e-4, compute sum4 and sum5 separately

expx2 = std::exp(-x * x);

const double exp2ax = std::exp((2 * a) * x), expm2ax = 1 / exp2ax;

for (int n = 1; 1; ++n)

{

const double coef = expa2n2[n - 1] * expx2 / (a2 * (n * n) + y * y);

prod2ax *= exp2ax;

prodm2ax *= expm2ax;

sum1 += coef;

sum2 += coef * prodm2ax;

sum4 += (coef * prodm2ax) * (a * n);

sum3 += coef * prod2ax;

sum5 += (coef * prod2ax) * (a * n);

// test convergence via sum5, since this sum has the slowest decay

if ((coef * prod2ax) * (a * n) < relerr * sum5)

break;

}

}

}

else

{ // relerr != DBL_EPSILON, compute exp(-a2*(n*n)) on the fly

const double exp2ax = std::exp((2 * a) * x), expm2ax = 1 / exp2ax;

if (x < 5e-4)

{ // compute sum4 and sum5 together as sum5-sum4

const double x2 = x * x;

expx2 = 1 - x2 * (1 - 0.5 * x2); // exp(-x*x) via Taylor

for (int n = 1; 1; ++n)

{

const double coef = exp(-a2 * (n * n)) * expx2 / (a2 * (n * n) + y * y);

prod2ax *= exp2ax;

prodm2ax *= expm2ax;

sum1 += coef;

sum2 += coef * prodm2ax;

sum3 += coef * prod2ax;

// really = sum5 - sum4

sum5 += coef * (2 * a) * n * std::sinh((2 * a) * n * x);

// test convergence via sum3

if (coef * prod2ax < relerr * sum3)

break;

}

}

else

{ // x > 5e-4, compute sum4 and sum5 separately

expx2 = std::exp(-x * x);

for (int n = 1; 1; ++n)

{

const double coef = std::exp(-a2 * (n * n)) * expx2 / (a2 * (n * n) + y * y);

prod2ax *= exp2ax;

prodm2ax *= expm2ax;

sum1 += coef;

sum2 += coef * prodm2ax;

sum4 += (coef * prodm2ax) * (a * n);

sum3 += coef * prod2ax;

sum5 += (coef * prod2ax) * (a * n);

// test convergence via sum5, since this sum has the slowest decay

if ((coef * prod2ax) * (a * n) < relerr * sum5)

break;

}

}

}

const double expx2erfcxy = // avoid spurious overflow for large negative y

y > -6 // for y < -6, erfcx(y) = 2*exp(y*y) to double precision

? expx2 * erfcx(y)

: 2 * std::exp(y * y - x * x);

if (y > 5)

{ // imaginary terms cancel

const double sinxy = std::sin(x * y);

ret = (expx2erfcxy - c * y * sum1) * cos(2 * x * y) + (c * x * expx2) * sinxy * sinc(x * y, sinxy);

}

else

{

double xs = std::real(z);

const double sinxy = std::sin(xs * y);

const double sin2xy = std::sin(2 * xs * y), cos2xy = std::cos(2 * xs * y);

const double coef1 = expx2erfcxy - c * y * sum1;

const double coef2 = c * xs * expx2;

ret = std::complex<double>(coef1 * cos2xy + coef2 * sinxy * sinc(xs * y, sinxy),

coef2 * sinc(2 * xs * y, sin2xy) - coef1 * sin2xy);

}

}

else

{ // x large: only sum3 & sum5 contribute (see above note)

if (std::isnan(x))

return std::complex<double>(x, x);

if (std::isnan(y))

return std::complex<double>(y, y);

ret = std::exp(-x * x); // |y| < 1e-10, so we only need exp(-x*x) term

// (round instead of ceil as in original paper; note that x/a > 1 here)

double n0 = std::floor(x / a + 0.5); // sum in both directions, starting at n0

double dx = a * n0 - x;

sum3 = std::exp(-dx * dx) / (a2 * (n0 * n0) + y * y);

sum5 = a * n0 * sum3;

double exp1 = std::exp(4 * a * dx), exp1dn = 1;

int dn;

for (dn = 1; n0 - dn > 0; ++dn)

{ // loop over n0-dn and n0+dn terms

double np = n0 + dn, nm = n0 - dn;

double tp = exp(-(a * dn + dx) * (a * dn + dx));

double tm = tp * (exp1dn *= exp1); // trick to get tm from tp

tp /= (a2 * (np * np) + y * y);

tm /= (a2 * (nm * nm) + y * y);

sum3 += tp + tm;

sum5 += a * (np * tp + nm * tm);

if (a * (np * tp + nm * tm) < relerr * sum5)

return ret

+ std::complex<double>((0.5 * c) * y * (sum2 + sum3),

(0.5 * c) * copysign(sum5 - sum4, std::real(z)));

;

}

while (1)

{ // loop over n0+dn terms only (since n0-dn <= 0)

double np = n0 + dn++;

double tp = std::exp(-(a * dn + dx) * (a * dn + dx)) / (a2 * (np * np) + y * y);

sum3 += tp;

sum5 += a * np * tp;

if (a * np * tp < relerr * sum5)

return ret

+ std::complex<double>((0.5 * c) * y * (sum2 + sum3),

(0.5 * c) * copysign(sum5 - sum4, std::real(z)));

;

}

}

return ret + std::complex<double>((0.5 * c) * y * (sum2 + sum3), (0.5 * c) * copysign(sum5 - sum4, std::real(z)));

;

{

if (x >= 0)

{

if (x > 45)

{ // continued-fraction expansion is faster

const double ispi = 0.56418958354775628694807945156; // 1 / sqrt(pi)

if (x > 5e7) // 1-term expansion, important to avoid overflow

return ispi / x;

/* 5-term expansion (rely on compiler for CSE), simplified from:

return ispi * ((x * x) * (x * x - 4.5) + 2) / (x * ((x * x) * (x * x - 5) + 3.75));

}

return w_im_y100(100 / (1 + x), x);

}

else

{ // = -FADDEEVA(w_im)(-x)

if (x < -45)

{ // continued-fraction expansion is faster

const double ispi = 0.56418958354775628694807945156; // 1 / sqrt(pi)

if (x < -5e7) // 1-term expansion, important to avoid overflow

return ispi / x;

/* 5-term expansion (rely on compiler for CSE), simplified from:

return ispi * ((x * x) * (x * x - 4.5) + 2) / (x * ((x * x) * (x * x - 5) + 3.75));

}

return -w_im_y100(100 / (1 - x), -x);

}

{

double x = std::real(z), y = std::imag(z);

if (y == 0)

return std::complex<double>(std::erf(x),

y); // preserve sign of 0

if (x == 0) // handle separately for speed & handling of y = Inf or NaN

return std::complex<double>(x, // preserve sign of 0

y * y > 720 ? (y > 0 ? Inf : -Inf) : std::exp(y * y) * scaled_w_im(y));

double mRe_z2 = (y - x) * (x + y); // Re(-z^2), being careful of overflow

double mIm_z2 = -2 * x * y; // Im(-z^2)

if (mRe_z2 < -750) // underflow

return (x >= 0 ? 1.0 : -1.0);

// Use Taylor series for small |z|, to avoid cancellation inaccuracy

// erf(z) = 2/sqrt(pi) * z * (1 - z^2/3 + z^4/10 - z^6/42 + z^8/216 + ...)

auto taylor = [&]() -> std::complex<double> {

std::complex<double> mz2 = std::complex<double>(mRe_z2, mIm_z2); // -z^2

return z

* (1.1283791670955125739

+ mz2

* (0.37612638903183752464

+ mz2

* (0.11283791670955125739

+ mz2 * (0.026866170645131251760 + mz2 * 0.0052239776254421878422))));

};

/* for small |x| and small |xy|,

auto taylor_erfi = [&]() -> std::complex<double> {

double x2 = x * x, y2 = y * y;

double expy2 = std::exp(y2);

return std::complex<double>(

expy2 * x

* (1.1283791670955125739 - x2 * (0.37612638903183752464 + 0.75225277806367504925 * y2)

+ x2 * x2 * (0.11283791670955125739 + y2 * (0.45135166683820502956 + 0.15045055561273500986 * y2))),

expy2

* (scaled_w_im(y)

- x2 * y * (1.1283791670955125739 - x2 * (0.56418958354775628695 + 0.37612638903183752464 * y2))));

};

/* Handle positive and negative x via different formulas,

if (x >= 0)

{

if (x < 8e-2)

{

if (std::fabs(y) < 1e-2)

return taylor();

else if (std::fabs(mIm_z2) < 5e-3 && x < 5e-3)

return taylor_erfi();

}

/* don't use complex exp function, since that will produce spurious NaN

return 1.0

- std::exp(mRe_z2)

* (std::complex<double>(std::cos(mIm_z2), std::sin(mIm_z2))

* scaled_w(std::complex<double>(-y, x), relerr));

}

else

{ // x < 0

if (x > -8e-2)

{ // duplicate from above to avoid fabs(x) call

if (std::fabs(y) < 1e-2)

return taylor();

else if (std::fabs(mIm_z2) < 5e-3 && x > -5e-3)

return taylor_erfi();

}

else if (std::isnan(x))

return std::complex<double>(NaN, y == 0 ? 0 : NaN);

/* don't use complex exp function, since that will produce spurious NaN

return std::exp(mRe_z2)

* (std::complex<double>(std::cos(mIm_z2), std::sin(mIm_z2))

* scaled_w(std::complex<double>(y, -x), relerr))

- 1.0;

}

{

std::complex<double> e = erf(std::complex<double>(-std::imag(z), std::real(z)), relerr);

return std::complex<double>(std::imag(e), -std::real(e));

{

double x = std::real(z), y = std::imag(z);

if (x == 0.)

return std::complex<double>(1,

/* handle y -> Inf limit manually, since

y * y > 720 ? (y > 0 ? -Inf : Inf) : -std::exp(y * y) * scaled_w_im(y));

if (y == 0.)

{

if (x * x > 750) // underflow

return std::complex<double>(x >= 0 ? 0.0 : 2.0,

-y); // preserve sign of 0

return std::complex<double>(x >= 0 ? std::exp(-x * x) * erfcx(x) : 2. - std::exp(-x * x) * erfcx(-x),

-y); // preserve sign of zero

}

double mRe_z2 = (y - x) * (x + y); // Re(-z^2), being careful of overflow

double mIm_z2 = -2 * x * y; // Im(-z^2)

if (mRe_z2 < -750) // underflow

return (x >= 0 ? 0.0 : 2.0);

if (x >= 0)

return std::exp(std::complex<double>(mRe_z2, mIm_z2)) * scaled_w(std::complex<double>(-y, x), relerr);

else

return 2.0 - std::exp(std::complex<double>(mRe_z2, mIm_z2)) * scaled_w(std::complex<double>(y, -x), relerr);

{

switch ((int)y100)

{

case 0: {

double t = 2 * y100 - 1;

return 0.28351593328822191546e-2

+ (0.28494783221378400759e-2

+ (0.14427470563276734183e-4

+ (0.10939723080231588129e-6

+ (0.92474307943275042045e-9

+ (0.89128907666450075245e-11 + 0.92974121935111111110e-13 * t) * t)

* t)

* t)

* t)

* t;

}

case 1: {

double t = 2 * y100 - 3;

return 0.85927161243940350562e-2

+ (0.29085312941641339862e-2

+ (0.15106783707725582090e-4

+ (0.11716709978531327367e-6

+ (0.10197387816021040024e-8

+ (0.10122678863073360769e-10 + 0.10917479678400000000e-12 * t) * t)

* t)

* t)

* t)

* t;

}

case 2: {

double t = 2 * y100 - 5;

return 0.14471159831187703054e-1

+ (0.29703978970263836210e-2

+ (0.15835096760173030976e-4

+ (0.12574803383199211596e-6

+ (0.11278672159518415848e-8

+ (0.11547462300333495797e-10 + 0.12894535335111111111e-12 * t) * t)

* t)

* t)

* t)

* t;

}

case 3: {

double t = 2 * y100 - 7;

return 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;

}

case 4: {

double t = 2 * y100 - 9;

return 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;

}

case 5: {

double t = 2 * y100 - 11;

return 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;

}

case 6: {

double t = 2 * y100 - 13;

return 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;

}

case 7: {

double t = 2 * y100 - 15;

return 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;

}

case 8: {

double t = 2 * y100 - 17;

return 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;

}

case 9: {

double t = 2 * y100 - 19;

return 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;

}

case 10: {

double t = 2 * y100 - 21;

return 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;

}

case 11: {

double t = 2 * y100 - 23;

return 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;

}

case 12: {

double t = 2 * y100 - 25;

return 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;

}

case 13: {

double t = 2 * y100 - 27;

return 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;

}

case 14: {

double t = 2 * y100 - 29;

return 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;

}

case 15: {

double t = 2 * y100 - 31;

return 0.10532778218603311137e0

+ (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;

}

case 16: {

double t = 2 * y100 - 33;

return 0.11380523107427108222e0

+ (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;

}

case 17: {

double t = 2 * y100 - 35;

return 0.12257529703447467345e0

+ (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;

}

case 18: {

double t = 2 * y100 - 37;

return 0.13166276955656699478e0

+ (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;

}

case 19: {

double t = 2 * y100 - 39;

return 0.14109647869803356475e0

+ (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;

}

case 20: {

double t = 2 * y100 - 41;

return 0.15091057940548936603e0

+ (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;

}

case 21: {

double t = 2 * y100 - 43;

return 0.16114648116017010770e0

+ (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;

}

case 22: {

double t = 2 * y100 - 45;

return 0.17185551279680451144e0

+ (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;

}

case 23: {

double t = 2 * y100 - 47;

return 0.18310194559815257381e0

+ (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;

}

case 24: {

double t = 2 * y100 - 49;

return 0.19496527191546630345e0

+ (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;

}

case 25: {

double t = 2 * y100 - 51;

return 0.20754006813966575720e0

+ (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;

}

case 26: {

double t = 2 * y100 - 53;

return 0.22093185554845172146e0

+ (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;

}

case 27: {

double t = 2 * y100 - 55;

return 0.23524827304057813918e0

+ (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;

}

case 28: {

double t = 2 * y100 - 57;

return 0.25058626331812744775e0

+ (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;

}

case 29: {

double t = 2 * y100 - 59;

return 0.26701724900280689785e0

+ (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;

}

case 30: {

double t = 2 * y100 - 61;

return 0.28457293586253654144e0

+ (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;

}

case 31: {

double t = 2 * y100 - 63;

return 0.30323425595617385705e0

+ (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;

}

case 32: {

double t = 2 * y100 - 65;

return 0.32292521181517384379e0

+ (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;

}

case 33: {

double t = 2 * y100 - 67;

return 0.34351233557911753862e0

+ (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;

}

case 34: {

double t = 2 * y100 - 69;

return 0.36480946642143669093e0

+ (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;

}

case 35: {

double t = 2 * y100 - 71;

return 0.38658679935694939199e0

+ (0.10966119158288804999e-1

+ (0.27521612041849561426e-4

+ (-0.57132774537670953638e-5