diffusionFactor2.cc

#include <assert.h>
#include <float.h>
#include <math.h>

#include "diffusionFactor2.h"

#ifdef __cplusplus
extern "C" {
#endif

double diffusionFactor2(double x)
{
    static const double a_m160[] = {
        0.03043351465358116,
        0.0052408158879836035,
        0.00067603436947200957,
        9.6906016427047165e-05,
        1.4589693269296478e-05,
        2.2599709015760758e-06,
        3.5665284682141206e-07,
        5.7029364769451708e-08,
        9.2088036202491301e-09,
        1.4983042248000036e-09,
        2.4525808290580496e-10,
        4.0345255474062071e-11,
        6.6642115253560996e-12,
        1.1046315810533649e-12,
        1.8364675844364025e-13,
        3.0610775248382879e-14,
        5.1138925270347952e-15,
        8.5606756200932375e-16,
        1.4354550385902103e-16,
        2.4088043958724545e-17,
        3.9385015302212454e-18
    };
    static const double a_m80[] = {
        0.042875904689467288,
        0.007209426563241781,
        0.00090927673972108591,
        0.00012753758747126521,
        1.8798771299551956e-05,
        2.8520962036009132e-06,
        4.4099489945542269e-07,
        6.9109593230610286e-08,
        1.0939649685577586e-08,
        1.7452440700629315e-09,
        2.8016949352042695e-10,
        4.5207291357665451e-11,
        7.3257761543250868e-12,
        1.1914500075596096e-12,
        1.9438149126458507e-13,
        3.1799097292178623e-14,
        5.21450648864145e-15,
        8.5692233925781008e-16,
        1.410682551675202e-16,
        2.3241777374612997e-17,
        3.7334187197300791e-18
    };
    static const double a_m40[] = {
        0.059963175073054216,
        0.0098993381250291525,
        0.00122888698356201,
        0.00016990303381926204,
        2.4711138774925709e-05,
        3.7023702466966373e-06,
        5.6570403392577033e-07,
        8.7655090280002247e-08,
        1.3725736055262137e-08,
        2.1670407838754021e-09,
        3.4441010877238304e-10,
        5.503729028989188e-11,
        8.8355137250870296e-12,
        1.4239951345647648e-12,
        2.3028067267343693e-13,
        3.7350324223653444e-14,
        6.0739269297352238e-15,
        9.9008600229174862e-16,
        1.6170036436654494e-16,
        2.6440657802272236e-17,
        4.2240524707015702e-18
    };
    static const double a_m20[] = {
        0.08318426223480628,
        0.013313268374736007,
        0.0016098032418279938,
        0.0002174058941277162,
        3.0949072035199655e-05,
        4.5456853410547137e-06,
        6.817555154282486e-07,
        1.0380180476908638e-07,
        1.598661375927885e-08,
        2.4844869252713606e-09,
        3.8896562707642755e-10,
        6.1269379308596461e-11,
        9.7012709313062235e-12,
        1.5429541447780714e-12,
        2.4635886347221512e-13,
        3.9470588217430402e-14,
        6.3431632089533856e-15,
        1.0222140426509087e-15,
        1.6510721155020624e-16,
        2.6702915336703069e-17,
        4.2228625910851895e-18
    };
    static const double a_m9[] = {
        0.11625317175453577,
        0.020230500432115689,
        0.0026841163757812584,
        0.000400018150379431,
        6.309982686849994e-05,
        1.0302721242097596e-05,
        1.7222546669302032e-06,
        2.9292305311565961e-07,
        5.0491094893621118e-08,
        8.7969204438940848e-09,
        1.5462655991544518e-09,
        2.738236970764213e-10,
        4.8801192788924355e-11,
        8.7458325364765392e-12,
        1.5750490556619191e-12,
        2.8488758332022055e-13,
        5.1730233646619079e-14,
        9.4263847564666446e-15,
        1.7230889540925401e-15,
        3.1555171809480902e-16,
        5.6096572002517297e-17
    };
    static const double a_m4[] = {
        0.16285593998638476,
        0.026057727371401426,
        0.0032444926850706221,
        0.00045948600696913881,
        6.9496192939853222e-05,
        1.0955503986269055e-05,
        1.7780828437199838e-06,
        2.949830852521079e-07,
        4.9791351199250328e-08,
        8.5237875106074325e-09,
        1.4764707938413207e-09,
        2.5832734813754341e-10,
        4.559091982556144e-11,
        8.1073373803106564e-12,
        1.451402124984054e-12,
        2.6139121010652892e-13,
        4.73285945069507e-14,
        8.6112199314614046e-15,
        1.5735626312627608e-15,
        2.8836567129895304e-16,
        5.1341893725263162e-17
    };
    static const double a_001[] = {
        0.21615320591176612,
        0.025667470909059423,
        0.0024591308682437364,
        0.00027403764849678023,
        3.3136150434335862e-05,
        4.2266141140533153e-06,
        5.6029476765066014e-07,
        7.6496075117207105e-08,
        1.0691490149984854e-08,
        1.523193162245219e-09,
        2.2050157405036199e-10,
        3.2355844878476651e-11,
        4.8033934459803503e-12,
        7.2033066644637672e-13,
        1.0898340864848178e-13,
        1.6618152037786782e-14,
        2.5516917551314589e-15,
        3.9436143824304815e-16,
        6.1212893013346769e-17,
        9.5198436290933311e-18,
        1.4732726395909797e-18
    };
    static const double a_01[] = {
        0.26911436547621914,
        0.026550990939845882,
        0.0022316852768348933,
        0.00022589510019385264,
        2.5417645370067592e-05,
        3.0698601195676094e-06,
        3.9028809370426237e-07,
        5.1592053940669445e-08,
        7.0318076195644763e-09,
        9.8227275994005179e-10,
        1.4000668464021816e-10,
        2.0293168920651209e-11,
        2.9832923150793548e-12,
        4.4390055092774891e-13,
        6.6741506423884432e-14,
        1.0125993313317228e-14,
        1.5485891958647422e-15,
        2.3864645069121559e-16,
        3.6898045899269329e-17,
        5.7191664598610359e-18,
        9.0253377275058221e-19
    };
    static const double a_03[] = {
        0.32012926467434183,
        0.023439827535842644,
        0.0015865814867846795,
        0.00013630377416242785,
        1.3469276418857162e-05,
        1.4606555057838224e-06,
        1.6914720156044868e-07,
        2.0557671755005003e-08,
        2.5920464655707826e-09,
        3.3633544080523081e-10,
        4.4653353036645457e-11,
        6.0400254998947521e-12,
        8.2972326202977265e-13,
        1.154690027869403e-13,
        1.6247759534536112e-14,
        2.3080560383715018e-15,
        3.305887309932944e-16,
        4.7907538138267221e-17,
        6.8349911957107006e-18,
        9.2221720504868209e-19,
        1.6085559017381665e-19
    };
    static const double a_06[] = {
        0.37178370164960722,
        0.026059013245511079,
        -0.0002864235290850683,
        0.00012237029300851102,
        -7.0224378260013625e-06,
        1.5905262620054674e-06,
        -1.4930143886968221e-07,
        2.8226054854609916e-08,
        -3.328856376657302e-09,
        5.8085035152287293e-10,
        -7.7509041182665165e-11,
        1.3007553389692023e-11,
        -1.8661611001232423e-12,
        3.0731892699252263e-13,
        -4.6112636278849983e-14,
        7.5299944747619651e-15,
        -1.1630752991000809e-15,
        1.8961405279899465e-16,
        -2.9954957165528079e-17,
        4.7469339760220999e-18,
        -7.5942553956685555e-19
    };
    static const double a_08[] = {
        0.41586503319497564,
        0.018567733719387337,
        0.00042592333826135716,
        5.3990025316010915e-05,
        4.3724286495273355e-06,
        5.5705473207608488e-07,
        6.539057613016484e-08,
        8.7200163110557916e-09,
        1.1682823821158203e-09,
        1.636448568110036e-10,
        2.3291521608022997e-11,
        3.3865645154038974e-12,
        4.9918820162545008e-13,
        7.455349349405021e-14,
        1.1250578789935639e-14,
        1.7134627876611151e-15,
        2.6308036643673566e-16,
        4.090604446606768e-17,
        6.2393253592796771e-18,
        9.076966402386084e-19,
        1.2729695150164627e-19
    };
    static const double a_09[] = {
        0.44619746616271777,
        0.011606122045773676,
        0.00035382736799969433,
        3.012684887954366e-05,
        2.9397219654356154e-06,
        3.3575439242470255e-07,
        4.1375904917077875e-08,
        5.3890684095233641e-09,
        7.2970735788821398e-10,
        1.0175269943753691e-10,
        1.4516987614077291e-11,
        2.1094516417608354e-12,
        3.1116140789980548e-13,
        4.6477701348211444e-14,
        7.0163459512457254e-15,
        1.0688336066905224e-15,
        1.6414562748056135e-16,
        2.5650416397386227e-17,
        3.7985797543152851e-18,
        5.040249385185589e-19,
        8.6155351206437401e-20
    };
    static const double a_095[] = {
        0.46545117238173106,
        0.007480253686070166,
        0.00023685523641974401,
        1.8559830233775195e-05,
        1.8375929615348941e-06,
        2.0900558506730081e-07,
        2.5786128743067678e-08,
        3.3597953874739569e-09,
        4.5518527834562034e-10,
        6.3503259789894595e-11,
        9.0640125845768299e-12,
        1.3176135039396327e-12,
        1.9442927086041339e-13,
        2.905118946744786e-14,
        4.3868910860096161e-15,
        6.6847388446404177e-16,
        1.0268685290309733e-16,
        1.6142995918185957e-17,
        2.2877956556316148e-18,
        3.284874439256677e-19,
        6.1631730638312899e-20,
    };
    static const double a_098[] = {
        0.47894250145124079,
        0.0059759053641499169,
        0.0002439042351159064,
        2.4342790265704914e-05,
        3.1670778786592217e-06,
        4.7262157934345584e-07,
        7.6542449086696733e-08,
        1.3091357122618394e-08,
        2.3281213799761949e-09,
        4.2632939718401429e-10,
        7.9870973243087215e-11,
        1.5239290024792436e-11,
        2.9514604477765125e-12,
        5.7879924914578889e-13,
        1.1471160754501311e-13,
        2.2941510164376882e-14,
        4.6243920254769254e-15,
        9.3884809194448653e-16,
        1.915352848629764e-16,
        3.913969842672671e-17,
        7.7262311958312254e-18
    };
    static const double a_0992[] = {
        0.48839701704992866,
        0.0033270028409703319,
        0.00013195084631676041,
        1.306268834728929e-05,
        1.7028554762349209e-06,
        2.5444557219946899e-07,
        4.1248103032871373e-08,
        7.0599914091042044e-09,
        1.2562440924468296e-09,
        2.3015181121949108e-10,
        4.3134409001373342e-11,
        8.2326416940641698e-12,
        1.5948908008379551e-12,
        3.1284246078681052e-13,
        6.2014919953991073e-14,
        1.2404837202083042e-14,
        2.5009175348204652e-15,
        5.0791515917078665e-16,
        1.0354453426454569e-16,
        2.114065164659533e-17,
        4.1964432301113054e-18
    };
    static const double a_0997[] = {
        0.49372640256336159,
        0.0019275134609139084,
        7.954094205087676e-05,
        8.4007488344470155e-06,
        1.1708253340623425e-06,
        1.8694426948753732e-07,
        3.237527785670156e-08,
        5.9189494781715684e-09,
        1.1248855336792873e-09,
        2.2009914017580867e-10,
        4.4053532061214941e-11,
        8.9791854920333388e-12,
        1.8576352556273172e-12,
        3.8911676511967414e-13,
        8.2370166511819811e-14,
        1.759455160104699e-14,
        3.7878100127352627e-15,
        8.2135542580260201e-16,
        1.7890045740290426e-16,
        3.9055801830046285e-17,
        8.2021830423836414e-18
    };
    static const double a_0999[] = {
        0.49679382230913094,
        0.0010948643551945106,
        4.9339775513703531e-05,
        5.8106148827030648e-06,
        9.0297933642160258e-07,
        1.6071611723980855e-07,
        3.1022996460780713e-08,
        6.3215527873357896e-09,
        1.3390319160991274e-09,
        2.9201310679876942e-10,
        6.5142659671346827e-11,
        1.479866018397699e-11,
        3.4122968226239452e-12,
        7.9664971347960851e-13,
        1.879570420951645e-13,
        4.474770189598618e-14,
        1.0737054175239111e-14,
        2.5943473663175315e-15,
        6.3022348996212765e-16,
        1.5338621163551473e-16,
        3.5530854052596389e-17
    };
    static const double a_09999[] = {
        1.0000002221483859,
        -1.3929527778011679e-07,
        -3.5807158709469604e-08,
        1.5121654311494091e-09,
        5.3970512437778956e-10,
        1.749186943964158e-10,
        6.0104289187874681e-11,
        2.1852250560483192e-11,
        8.3138373108038966e-12,
        3.2790311109965967e-12,
        1.3312071301023399e-12,
        5.5334276553985535e-13,
        2.3455336032391335e-13,
        1.0107390081969312e-13,
        4.4169493005249048e-14,
        1.9535843986653469e-14,
        8.7294821875382235e-15,
        3.9322541378815156e-15,
        1.7738812444397418e-15,
        7.8254356946136493e-16,
        2.9694490500757159e-16
    };
    static const double a_099999[] = {
        1.0000000253876808,
        -2.3380627410207544e-08,
        7.6175735905551203e-10,
        1.9722123521977164e-10,
        4.5532178143588751e-11,
        1.2850367703835401e-11,
        4.1082667969836677e-12,
        1.4271380972653781e-12,
        5.2590327616044593e-13,
        2.0252111309884293e-13,
        8.069314582432292e-14,
        3.3036837980890544e-14,
        1.3828425631761439e-14,
        5.8957520165539961e-15,
        2.5527953180570246e-15,
        1.1198233767191412e-15,
        4.9702409020479538e-16,
        2.2307975985638059e-16,
        9.959171384402562e-17,
        4.3930839455615036e-17,
        1.6866442724945631e-17
    };
    static const double a_09999996[] = {
        1.0000000014245567,
        -1.5114761508160367e-09,
        1.0896504763822493e-10,
        1.7720251693123193e-11,
        4.6966472388221927e-12,
        1.6006987608743089e-12,
        6.3030644175405775e-13,
        2.7258921186736259e-13,
        1.2589734975431724e-13,
        6.1043592262679418e-14,
        3.0724725219390323e-14,
        1.5930611038331021e-14,
        8.4606542588704426e-15,
        4.5831136091932846e-15,
        2.5239245630639817e-15,
        1.4093595668813952e-15,
        7.9412388123700772e-16,
        4.5025109219702187e-16,
        2.5134387791958008e-16,
        1.343533006409061e-16,
        5.8508841166493046e-17
    };
    static const double a_1[] = {
        1.0000000000211662,
        -2.4059848555846177e-11,
        2.281829848666116e-12,
        3.737551769401022e-13,
        1.1829146499963193e-13,
        5.011304284900719e-14,
        2.5138290919778954e-14,
        1.4098805002109372e-14,
        8.5589707444449261e-15,
        5.5099928835836536e-15,
        3.7094376842670246e-15,
        2.5861164648247456e-15,
        1.8543832557467907e-15,
        1.353781298372793e-15,
        1.0024481658126094e-15,
        7.4729151025312198e-16,
        5.559453154118945e-16,
        4.0537544799109808e-16,
        2.8139951100112463e-16,
        1.7792919295671934e-16,
        8.7158238215986494e-17
    };

    assert(x >= 0.0);
    assert(x <= 1.0);

    if (x == 0.0)
        return 0.0;
    if (x == 1.0)
        return 0.5;

    static const double minratio = 1.e-323;
    const double xin = x;
    double xmin, xmax;
    unsigned degree = 20, near1 = 0;
    const double* a = 0;

    if (x < 1.e-160)
    {
        a = a_m160;
        xmin = log(minratio);
        xmax = log(1.e-160);
        if (x < minratio)
            x = minratio;
        x = log(x);
    }
    else if (x < 1.e-80)
    {
        a = a_m80;
        xmin = log(1.e-160);
        xmax = log(1.e-80);
        x = log(x);
    }
    else if (x < 1.e-40)
    {
        a = a_m40;
        xmin = log(1.e-80);
        xmax = log(1.e-40);
        x = log(x);
    }
    else if (x < 1.e-20)
    {
        a = a_m20;
        xmin = log(1.e-40);
        xmax = log(1.e-20);
        x = log(x);
    }
    else if (x < 1.e-9)
    {
        a = a_m9;
        xmin = log(1.e-20);
        xmax = log(1.e-9);
        x = log(x);
    }
    else if (x < 0.0001)
    {
        a = a_m4;
        xmin = log(1.e-9);
        xmax = log(0.0001);
        x = log(x);
    }
    else if (x < 0.01)
    {
        a = a_001;
        xmin = log(0.0001);
        xmax = log(0.01);
        x = log(x);
    }
    else if (x < 0.1)
    {
        a = a_01;
        xmin = log(0.01);
        xmax = log(0.1);
        x = log(x);
    }
    else if (x < 0.3)
    {
        a = a_03;
        xmin = log(0.1);
        xmax = log(0.3);
        x = log(x);
    }
    else if (x < 0.6)
    {
        a = a_06;
        xmin = 0.3;
        xmax = 0.6;
    }
    else if (x < 0.8)
    {
        a = a_08;
        xmin = 0.6;
        xmax = 0.8;
    }
    else if (x < 0.9)
    {
        a = a_09;
        xmin = 0.8;
        xmax = 0.9;
    }
    else if (x < 0.95)
    {
        a = a_095;
        xmin = 0.9;
        xmax = 0.95;
    }
    else if (x < 0.98)
    {
        a = a_098;
        xmin = 0.95;
        xmax = 0.98;
    }
    else if (x < 0.992)
    {
        a = a_0992;
        xmin = 0.98;
        xmax = 0.992;
    }
    else if (x < 0.997)
    {
        a = a_0997;
        xmin = 0.992;
        xmax = 0.997;
    }
    else if (x < 0.999)
    {
        a = a_0999;
        xmin = 0.997;
        xmax = 0.999;
    }
    else if (x < 0.9999)
    {
        a = a_09999;
        xmin = 0.999;
        xmax = 0.9999;
        near1 = 1;
    }
    else if (x < 0.99999)
    {
        a = a_099999;
        xmin = 0.9999;
        xmax = 0.99999;
        near1 = 1;
    }
    else if (x < 0.9999996)
    {
        a = a_09999996;
        xmin = 0.99999;
        xmax = 0.9999996;
        near1 = 1;
    }
    else
    {
        a = a_1;
        xmin = 0.9999996;
        xmax = 1.0;
        near1 = 1;
    }

    x = 2.0*(x - xmin)/(xmax - xmin) - 1.0;
    const long double twox = 2.0L*x;
    long double rp2 = 0.0L, rp1 = 0.0L, r = 0.0L;
    unsigned k;
    for (k = degree; k > 0U; --k)
    {
        r = twox*rp1 - rp2 + a[k];
        rp2 = rp1;
        rp1 = r;
    }
    double cheb = x*rp1 - rp2 + a[0];
    if (near1)
        cheb /= (pow(3.0/4.0*(1.0 - xin), 2.0/3.0) + 2.0);
    return cheb;
}

#ifdef __cplusplus
}
#endif