reduce number of digits to machine precision (16)

Expander · Oct 11, 2018 · 78bfd7c · 78bfd7c
1 parent 8167134
commit 78bfd7c
Show file tree

Hide file tree

Showing 5 changed files with 111 additions and 110 deletions.
diff --git a/src/Li.cpp b/src/Li.cpp
@@ -23,41 +23,41 @@ namespace {
    /// Bernoulli numbers B0, ..., B49
    /// Table[BernoulliB[n], {n,0,49}]
    const double bernoulli[N] = {
-      1, -0.5               , 0.16666666666666666     , 0,
-      -0.03333333333333333  , 0, 0.023809523809523808 , 0,
-      -0.03333333333333333  , 0, 0.07575757575757576  , 0,
-      -0.2531135531135531   , 0, 1.1666666666666667   , 0,
-      -7.092156862745098    , 0, 54.971177944862156   , 0,
-      -529.1242424242424    , 0, 6192.123188405797    , 0,
-      -86580.25311355312    , 0, 1.4255171666666667e6 , 0,
-      -2.7298231067816094e7 , 0, 6.015808739006424e8  , 0,
-      -1.5116315767092157e10, 0, 4.296146430611667e11 , 0,
-      -1.3711655205088334e13, 0, 4.883323189735932e14 , 0,
-      -1.9296579341940068e16, 0, 8.416930475736827e17 , 0,
-      -4.0338071854059454e19, 0, 2.1150748638081993e21, 0,
-      -1.2086626522296526e23, 0
+      1, -0.5               , 1./6.                   , 0,
+      -3.333333333333333e-02, 0, 2.380952380952381e-02, 0,
+      -3.333333333333333e-02, 0, 7.575757575757576e-02, 0,
+      -2.531135531135531e-01, 0, 1.166666666666667e+00, 0,
+      -7.092156862745098e+00, 0, 5.497117794486215e+01, 0,
+      -5.291242424242424e+02, 0, 6.192123188405797e+03, 0,
+      -8.658025311355312e+04, 0, 1.425517166666667e+06, 0,
+      -2.729823106781609e+07, 0, 6.015808739006424e+08, 0,
+      -1.511631576709215e+10, 0, 4.296146430611667e+11, 0,
+      -1.371165520508833e+13, 0, 4.883323189735932e+14, 0,
+      -1.929657934194006e+16, 0, 8.416930475736827e+17, 0,
+      -4.033807185405945e+19, 0, 2.115074863808199e+21, 0,
+      -1.208662652229652e+23, 0
    };
 
    /// 1/n! for n = 1, ..., 50
    /// Table[1/Factorial[n], {n,1,50}]
    const double fac_inv[N] = {
-      1                     , 0.5                   , 0.16666666666666666    ,
-      0.041666666666666664  , 0.008333333333333333  , 0.001388888888888889   ,
-      0.0001984126984126984 , 0.0000248015873015873 , 2.7557319223985893e-6  ,
-      2.755731922398589e-7  , 2.505210838544172e-8  , 2.08767569878681e-9    ,
-      1.6059043836821613e-10, 1.1470745597729725e-11, 7.647163731819816e-13  ,
-      4.779477332387385e-14 , 2.8114572543455206e-15, 1.5619206968586225e-16 ,
-      8.22063524662433e-18  , 4.110317623312165e-19 , 1.9572941063391263e-20 ,
-      8.896791392450574e-22 , 3.8681701706306835e-23, 1.6117375710961184e-24 ,
-      6.446950284384474e-26 , 2.4795962632247972e-27, 9.183689863795546e-29  ,
-      3.279889237069838e-30 , 1.1309962886447718e-31, 3.7699876288159054e-33 ,
-      1.2161250415535181e-34, 3.800390754854744e-36 , 1.151633562077195e-37  ,
-      3.387157535521162e-39 , 9.67759295863189e-41  , 2.688220266286636e-42  ,
-      7.265460179153071e-44 , 1.911963205040282e-45 , 4.902469756513544e-47  ,
-      1.225617439128386e-48 , 2.9893108271424046e-50, 7.117406731291439e-52  ,
-      1.6552108677421951e-53, 3.7618428812322616e-55, 8.359650847182804e-57  ,
-      1.817315401561479e-58 , 3.866628513960594e-60 , 8.055476070751238e-62  ,
-      1.643974708316579e-63 , 3.2879494166331584e-65
+      1                    , 0.5                  , 1./6.                ,
+      4.166666666666666e-02, 8.333333333333333e-03, 1.388888888888889e-03,
+      1.984126984126984e-04, 2.480158730158730e-05, 2.755731922398589e-06,
+      2.755731922398589e-07, 2.505210838544172e-08, 2.087675698786810e-09,
+      1.605904383682161e-10, 1.147074559772973e-11, 7.647163731819816e-13,
+      4.779477332387385e-14, 2.811457254345521e-15, 1.561920696858623e-16,
+      8.220635246624330e-18, 4.110317623312165e-19, 1.957294106339126e-20,
+      8.896791392450574e-22, 3.868170170630684e-23, 1.611737571096118e-24,
+      6.446950284384474e-26, 2.479596263224797e-27, 9.183689863795546e-29,
+      3.279889237069838e-30, 1.130996288644772e-31, 3.769987628815905e-33,
+      1.216125041553518e-34, 3.800390754854744e-36, 1.151633562077195e-37,
+      3.387157535521162e-39, 9.677592958631890e-41, 2.688220266286636e-42,
+      7.265460179153071e-44, 1.911963205040282e-45, 4.902469756513544e-47,
+      1.225617439128386e-48, 2.989310827142405e-50, 7.117406731291439e-52,
+      1.655210867742195e-53, 3.761842881232262e-55, 8.359650847182804e-57,
+      1.817315401561479e-58, 3.866628513960594e-60, 8.055476070751238e-62,
+      1.643974708316579e-63, 3.287949416633158e-65
    };
 
    bool is_close(double a, double b, double eps = epsilon)
@@ -280,7 +280,7 @@ namespace {
 double Cl(long n, double x)
 {
    using std::exp;
-   const double PI = 3.1415926535897932384626433832795;
+   const double PI = 3.141592653589793;
    const std::complex<double> i(0.,1.);
    const std::complex<double> li = Li(n, exp(i*x));
 
@@ -333,7 +333,7 @@ std::complex<double> Li(long n, const std::complex<double>& z)
    if (std::abs(z) <= 1.) {
       u = -clog(1. - z);
    } else { // az > 1.
-      const double PI = 3.1415926535897932384626433832795;
+      const double PI = 3.141592653589793;
       const std::complex<double> IPI2(0.,2*PI);
       const std::complex<double> lnz = clog(-z);
       u = -clog(1. - 1./z);

diff --git a/src/Li3.cpp b/src/Li3.cpp
@@ -39,7 +39,7 @@ namespace {
 double Cl3(double x)
 {
    using std::exp;
-   const double PI = 3.1415926535897932384626433832795;
+   const double PI = 3.141592653589793;
    const std::complex<double> i(0.,1.);
 
    while (x >= 2*PI)
@@ -65,13 +65,13 @@ std::complex<double> Li3(const std::complex<double>& z)
 
    const double bf[18] = {
       1., -3./8., 17./216., -5./576.,
-      0.00012962962962962962962962962962963,  0.000081018518518518518518518518518519,
-     -3.4193571608537594932152755282007e-6 , -1.3286564625850340136054421768707e-6,
-      8.6608717561098513479465860418241e-8 ,  2.5260875955320399764844209288654e-8,
-     -2.1446944683640647609338850757365e-9 , -5.14011062201297891533581769272e-10,
-      5.2495821146008294363940888085581e-11,  1.0887754406636318375372971570425e-11,
-     -1.2779396094493695305581831754072e-12, -2.3698241773087452099797778810124e-13,
-      3.1043578879654622942847532704656e-14,  5.2617586299125060841318392511225e-15
+      1.296296296296296e-04,  8.101851851851851e-05,
+     -3.419357160853759e-06, -1.328656462585034e-06,
+      8.660871756109851e-08,  2.526087595532039e-08,
+     -2.144694468364064e-09, -5.140110622012978e-10,
+      5.249582114600829e-11,  1.088775440663631e-11,
+     -1.277939609449369e-12, -2.369824177308745e-13,
+      3.104357887965462e-14,  5.261758629912506e-15
    };
 
    if (is_close(z, 0.))

diff --git a/src/Li4.cpp b/src/Li4.cpp
@@ -40,7 +40,7 @@ namespace {
 double Cl4(double x)
 {
    using std::exp;
-   const double PI = 3.1415926535897932384626433832795;
+   const double PI = 3.141592653589793;
    const std::complex<double> i(0.,1.);
 
    while (x >= 2*PI)
@@ -64,32 +64,33 @@ double Cl4(double x)
  */
 std::complex<double> Li4(const std::complex<double>& z)
 {
-   const double PI    = 3.1415926535897932384626433832795;
+   const double PI    = 3.141592653589793;
    const double PI2   = PI*PI;
    const double PI4   = PI2*PI2;
-   const double zeta4 = 1.0823232337111381915160036965412;
+   const double zeta4 = 1.082323233711138;
    static const int N = 40;
 
    const double bf[N] = {
-      1., -7./16., 0.11651234567901234567901234567901, -0.019820601851851851851851851851852,
-      0.0019279320987654320987654320987654  , -0.000031057098765432098765432098765432,
-     -0.000015624009114857835298392473643526,  8.4851235467732066371522153835079e-7  ,
-      2.2909616603189711445359383547004e-7  , -2.1832614218526916939615352313765e-8  ,
-     -3.8828248791720155722806620380777e-9  ,  5.4462921032203321182579858808232e-10 ,
-      6.9608052106827254078772334134121e-11 , -1.3375737686445215199578072203635e-11 ,
-     -1.2784852685266571604146246361574e-12 ,  3.2605628580248922428788418178217e-13 ,
-      2.3647571168618257362309504812439e-14 , -7.9231351220311617024299900711372e-15 ,
-     -4.3452915709984187250497371626475e-16 ,  1.9236270062535920116126875526753e-16 ,
-      7.8124143331959546707222938968737e-18 , -4.6718038448036555203176282428722e-18 ,
-     -1.3435344329812847856260222675894e-19 ,  1.1356826851347343244764698375938e-19 ,
-      2.1152756202432586847505983414192e-21 , -2.7642026334746517388281729253731e-21 ,
-     -2.7068176608240064256109059558195e-23 ,  6.7372044828628572143267161265626e-23 ,
-      1.3287265456683822975818009045013e-25 , -1.6443773056367826467816763114889e-24 ,
-      8.2836058999339341109829673400349e-27 ,  4.0190848495069350699709315007621e-26 ,
-     -4.5757138444848790382359734346537e-28 , -9.8364109094615127758320974982117e-28 ,
-      1.6900339556037851067729523121903e-29 ,  2.4104805563059808504664904164902e-29 ,
-     -5.4266127056714182501325034058929e-31 , -5.9142429588741767864337599966928e-31 ,
-      1.6232110901087370772711176143968e-32 ,  1.4527595437740275946132587316158e-32
+      1., -7./16.,
+      1.165123456790123e-01, -1.982060185185185e-02,
+      1.927932098765432e-03, -3.105709876543209e-05,
+     -1.562400911485783e-05,  8.485123546773206e-07,
+      2.290961660318971e-07, -2.183261421852691e-08,
+     -3.882824879172015e-09,  5.446292103220332e-10,
+      6.960805210682725e-11, -1.337573768644521e-11,
+     -1.278485268526657e-12,  3.260562858024892e-13,
+      2.364757116861825e-14, -7.923135122031161e-15,
+     -4.345291570998418e-16,  1.923627006253592e-16,
+      7.812414333195954e-18, -4.671803844803655e-18,
+     -1.343534432981284e-19,  1.135682685134734e-19,
+      2.115275620243258e-21, -2.764202633474651e-21,
+     -2.706817660824006e-23,  6.737204482862857e-23,
+      1.328726545668382e-25, -1.644377305636782e-24,
+      8.283605899933934e-27,  4.019084849506935e-26,
+     -4.575713844484879e-28, -9.836410909461512e-28,
+      1.690033955603785e-29,  2.410480556305980e-29,
+     -5.426612705671418e-31, -5.914242958874176e-31,
+      1.623211090108737e-32,  1.452759543774027e-32
    };
 
    if (is_close(z, 0.))
@@ -101,7 +102,7 @@ std::complex<double> Li4(const std::complex<double>& z)
       const std::complex<double> IPI(0.,PI);
       const std::complex<double> zm1 = z - 1.;
       const std::complex<double> lzm1 = clog(zm1);
-      const double zeta3 = 1.2020569031595942853997381615114;
+      const double zeta3 = 1.202056903159594;
       const double ceil = std::arg(zm1) > 0. ? 1. : 0.;
       const std::complex<double> cs[8] = {
          zeta3,

diff --git a/src/Li5.cpp b/src/Li5.cpp
@@ -40,7 +40,7 @@ namespace {
 double Cl5(double x)
 {
    using std::exp;
-   const double PI = 3.1415926535897932384626433832795;
+   const double PI = 3.141592653589793;
    const std::complex<double> i(0.,1.);
 
    while (x >= 2*PI)
@@ -59,33 +59,33 @@ double Cl5(double x)
  */
 std::complex<double> Li5(const std::complex<double>& z)
 {
-   const double PI    = 3.1415926535897932384626433832795;
+   const double PI    = 3.141592653589793;
    const double PI2   = PI*PI;
    const double PI4   = PI2*PI2;
-   const double zeta5 = 1.0369277551433699263313654864570;
+   const double zeta5 = 1.036927755143369;
    static const int N = 40;
 
    const double bf[N] = {
       1., -15./32.,
-      0.13953189300411522633744855967078   , -0.028633777006172839506172839506173  ,
-      0.0040317412551440329218106995884774 , -0.00033985018004115226337448559670782,
-      4.5445184621617666490944600374313e-6 ,  2.3916808048569011882908870275245e-6 ,
-     -1.2762692600122746588544351898175e-7 , -3.1628984306505932440256787279547e-8 ,
-      3.2848118445335191621518574238472e-9 ,  4.7613713995660579048319132897732e-10,
-     -8.0846898171909830256460260362332e-11, -7.238764858773720694682921583759e-12 ,
-      1.943976011517396849305564920936e-12 ,  1.025697840597723597188135594332e-13 ,
-     -4.6180551009884830180582086241066e-14, -1.1535857196470580036842511483419e-15,
-      1.0903545401333393987977088380966e-15,  2.3148136317292526394079710319009e-18,
-     -2.5669917043265292194334891993397e-17,  4.5708620607314969014495962686014e-19,
-      6.0366779613205705882356103311411e-19, -2.167762494406241295879417172184e-20 ,
-     -1.4194096615600165298332266882011e-20,  7.5020009506413862553237752161953e-22,
-      3.3387045395078397164371515925447e-22, -2.3060040442620347682521515135259e-23,
-     -7.8581732456894818904499064631503e-24,  6.6683453043738808548651370461306e-25,
-      1.8509156540925297189464979688365e-25, -1.8591529445174085584103184057636e-26,
-     -4.362974648034589044726608174371e-27 ,  5.0611076099529284482263489587811e-28,
-      1.0291918249756878203788897930001e-28, -1.3551391221018316615687785328304e-29,
-     -2.4294059612957382655924154095657e-30,  3.5851973966503705211516473805307e-31,
-      5.7379658161039720640084653867328e-32, -9.4003593624568734535277452038948e-33
+      1.395318930041152e-01, -2.863377700617283e-02,
+      4.031741255144032e-03, -3.398501800411522e-04,
+      4.544518462161766e-06,  2.391680804856901e-06,
+     -1.276269260012274e-07, -3.162898430650593e-08,
+      3.284811844533519e-09,  4.761371399566057e-10,
+     -8.084689817190983e-11, -7.238764858773720e-12,
+      1.943976011517396e-12,  1.025697840597723e-13,
+     -4.618055100988483e-14, -1.153585719647058e-15,
+      1.090354540133339e-15,  2.314813631729252e-18,
+     -2.566991704326529e-17,  4.570862060731496e-19,
+      6.036677961320570e-19, -2.167762494406241e-20,
+     -1.419409661560016e-20,  7.502000950641386e-22,
+      3.338704539507839e-22, -2.306004044262034e-23,
+     -7.858173245689481e-24,  6.668345304373880e-25,
+      1.850915654092529e-25, -1.859152944517408e-26,
+     -4.362974648034589e-27,  5.061107609952928e-28,
+      1.029191824975687e-28, -1.355139122101831e-29,
+     -2.429405961295738e-30,  3.585197396650370e-31,
+      5.737965816103972e-32, -9.400359362456873e-33
    };
 
    if (is_close(z, 0.))
@@ -97,7 +97,7 @@ std::complex<double> Li5(const std::complex<double>& z)
       const std::complex<double> IPI(0.,PI);
       const std::complex<double> zm1 = z - 1.;
       const std::complex<double> lzm1 = clog(zm1);
-      const double zeta3 = 1.2020569031595942853997381615114;
+      const double zeta3 = 1.202056903159594;
       const double ceil = std::arg(zm1) > 0. ? 1. : 0.;
       const std::complex<double> cs[8] = {
          PI4/90.,

diff --git a/src/Li6.cpp b/src/Li6.cpp
@@ -40,7 +40,7 @@ namespace {
 double Cl6(double x)
 {
    using std::exp;
-   const double PI = 3.1415926535897932384626433832795;
+   const double PI = 3.141592653589793;
    const std::complex<double> i(0.,1.);
 
    while (x >= 2*PI)
@@ -64,34 +64,34 @@ double Cl6(double x)
  */
 std::complex<double> Li6(const std::complex<double>& z)
 {
-   const double PI    = 3.1415926535897932384626433832795;
+   const double PI    = 3.141592653589793;
    const double PI2   = PI*PI;
    const double PI4   = PI2*PI2;
    const double PI6   = PI2*PI4;
-   const double zeta6 = 1.0173430619844491397145179297909;
+   const double zeta6 = 1.017343061984449;
    static const int N = 40;
 
    const double bf[N] = {
       1., -31./64.,
-      0.15241340877914951989026063100137    , -0.034365555877057613168724279835391  ,
-      0.0057174797239368998628257887517147  , -0.00068180453746570644718792866941015,
-      0.000049960361948734493114517016962133, -4.9166051196039047720253082216462e-7 ,
-     -3.063297516130216378535303043104e-7   ,  1.4414599270849095361253742144801e-8 ,
-      3.727243823092410657685992456646e-9   , -3.7300867345487607207797722915926e-10,
-     -5.1246526816085832434008764611572e-11 ,  9.0541930956636682886879744762089e-12,
-      6.7381882615512517077610754493801e-13 , -2.1215831150303135318527968929661e-13,
-     -6.8408811719011697663920451878191e-15 ,  4.8691178462005581320638758651292e-15,
-     -4.8439878499872504165088964747316e-18 , -1.1027104849107490937043030257605e-16,
-      3.3353796916939381662488941184074e-18 ,  2.4735307488641352979154098748714e-18,
-     -1.4370616434232492021688368713474e-19 , -5.5047110335098118061482643505979e-20,
-      4.7467713917327224979130984066262e-21 ,  1.2158387178068105224373981741629e-21,
-     -1.4107552403561850041424030907801e-22 , -2.6638831253268346596585643767712e-23,
-      3.9667657428631007976790022608178e-24 ,  5.7821697358543615311236619348113e-25,
-     -1.0787778063164257317299887699592e-25 , -1.2407397086756909899014773697759e-26,
-      2.8704117917893601704260952468409e-27 ,  2.6235553563029330616574752038361e-28,
-     -7.5229485465754127261588121413434e-29 , -5.4401788379624696182029193072231e-30,
-      1.9502579532510166379386222338166e-30 ,  1.0978494282205187996117821359701e-31,
-     -5.0149583574163009207446958541576e-32 , -2.1286737504392761053563380691739e-33
+      1.524134087791495e-01, -3.436555587705761e-02,
+      5.717479723936899e-03, -6.818045374657064e-04,
+      4.996036194873449e-05, -4.916605119603904e-07,
+     -3.063297516130216e-07,  1.441459927084909e-08,
+      3.727243823092410e-09, -3.730086734548760e-10,
+     -5.124652681608583e-11,  9.054193095663668e-12,
+      6.738188261551251e-13, -2.121583115030313e-13,
+     -6.840881171901169e-15,  4.869117846200558e-15,
+     -4.843987849987250e-18, -1.102710484910749e-16,
+      3.335379691693938e-18,  2.473530748864135e-18,
+     -1.437061643423249e-19, -5.504711033509811e-20,
+      4.746771391732722e-21,  1.215838717806810e-21,
+     -1.410755240356185e-22, -2.663883125326834e-23,
+      3.966765742863100e-24,  5.782169735854361e-25,
+     -1.078777806316425e-25, -1.240739708675690e-26,
+      2.870411791789360e-27,  2.623555356302933e-28,
+     -7.522948546575412e-29, -5.440178837962469e-30,
+      1.950257953251016e-30,  1.097849428220518e-31,
+     -5.014958357416300e-32, -2.128673750439276e-33
    };
 
    if (is_close(z, 0.))
@@ -103,8 +103,8 @@ std::complex<double> Li6(const std::complex<double>& z)
       const std::complex<double> IPI(0.,PI);
       const std::complex<double> zm1 = z - 1.;
       const std::complex<double> lzm1 = clog(zm1);
-      const double zeta3 = 1.2020569031595942853997381615114;
-      const double zeta5 = 1.0369277551433699263313654864570;
+      const double zeta3 = 1.202056903159594;
+      const double zeta5 = 1.036927755143369;
       const double ceil = std::arg(zm1) > 0. ? 1. : 0.;
       const std::complex<double> cs[8] = {
          zeta5,