diff --git a/src/Li.cpp b/src/Li.cpp index c9771ead..119fdbce 100644 --- 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 i(0.,1.); const std::complex li = Li(n, exp(i*x)); @@ -333,7 +333,7 @@ std::complex Li(long n, const std::complex& 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 IPI2(0.,2*PI); const std::complex lnz = clog(-z); u = -clog(1. - 1./z); diff --git a/src/Li3.cpp b/src/Li3.cpp index a5475d1d..36df037d 100644 --- 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 i(0.,1.); while (x >= 2*PI) @@ -65,13 +65,13 @@ std::complex Li3(const std::complex& 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 index 4dda5261..675d0321 100644 --- 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 i(0.,1.); while (x >= 2*PI) @@ -64,32 +64,33 @@ double Cl4(double x) */ std::complex Li4(const std::complex& 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 Li4(const std::complex& z) const std::complex IPI(0.,PI); const std::complex zm1 = z - 1.; const std::complex 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 cs[8] = { zeta3, diff --git a/src/Li5.cpp b/src/Li5.cpp index 299800bf..51b81b48 100644 --- 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 i(0.,1.); while (x >= 2*PI) @@ -59,33 +59,33 @@ double Cl5(double x) */ std::complex Li5(const std::complex& 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 Li5(const std::complex& z) const std::complex IPI(0.,PI); const std::complex zm1 = z - 1.; const std::complex 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 cs[8] = { PI4/90., diff --git a/src/Li6.cpp b/src/Li6.cpp index 66a2605b..d5364580 100644 --- 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 i(0.,1.); while (x >= 2*PI) @@ -64,34 +64,34 @@ double Cl6(double x) */ std::complex Li6(const std::complex& 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 Li6(const std::complex& z) const std::complex IPI(0.,PI); const std::complex zm1 = z - 1.; const std::complex 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 cs[8] = { zeta5,