Merge pull request #1624 from flintlib/catalan

Use new formulas from Jorge Zuniga to compute Catalan and zeta(3) faster

doc/source/constants.rst

@@ -72,9 +72,42 @@ Catalan's constant is computed using the hypergeometric series
 
 .. math ::
 
-    C = \frac{1}{64} \sum_{k=1}^{\infty} \frac{256^k (580k^2-184k+15)}{k^3(2k-1){6k\choose 3k}{6k\choose 4k}{4k\choose 2k}}
+    C = \frac{1}{768} \sum_{k=1}^{\infty} \frac{(-4096)^k P(k)}
+        {k^3 (2k-1)(3k-1)(3k-2)(6k-1)(6k-5) {5k \choose k} {10k \choose 5k} {12k \choose 6k}}
 
-given in [PP2010]_.
+where
+
+.. math ::
+
+    \begin{align}
+        P(k) & = -43203456k^6 + 92809152k^5 - 76613904k^4 \\
+             & + 30494304k^3 - 6004944k^2 + 536620^k - 17325,
+    \end{align}
+
+discovered by Zuniga [Zun2023]_.
+It was previously computed using a series given in [PP2010]_.
+
+Apery's constant
+-------------------------------------------------------------------------------
+
+Apery's constant `\zeta(3)` is computed using the hypergeometric series
+
+.. math ::
+
+    \zeta(3) = \frac{1}{48} \sum_{k=1}^{\infty} \frac{(-1)^{k-1} P(k)}{k^5 (2k-1)^3(3k-1)(3k-2)(4k-1)(4k-3)(6k-1)(6k-5){5k \choose k}{5k \choose 2k}{9k \choose 4k}{10k \choose 5k}{12k \choose 6k}}
+
+where
+
+.. math ::
+
+    \begin{align}
+        P(k) & = 1565994397644288k^{11} - 6719460725627136k^{10} + 12632254526031264k^9 \\
+             & - 13684352515879536k^8 + 9451223531851808k^7 - 4348596587040104k^6 \\
+             & + 1352700034136826k^5 - 282805786014979k^4 + 38721705264979k^3 \\
+             & - 3292502315430k^2 + 156286859400k - 3143448000,
+    \end{align}
+
+discovered by Zuniga [Zun2023]_.
 
 Khinchin's constant
 -------------------------------------------------------------------------------
@@ -114,16 +147,6 @@ from this formula. We don't use the reflection formula for the zeta function,
 as the arithmetic in Euler-Maclaurin summation is faster at `s = -1`
 than at `s = 2`.
 
-Apery's constant
--------------------------------------------------------------------------------
-
-Apery's constant `\zeta(3)` is computed using the hypergeometric series
-
-.. math ::
-
-    \zeta(3) = \frac{1}{64} \sum_{k=0}^\infty
-        (-1)^k (205k^2 + 250k + 77) \frac{(k!)^{10}}{[(2k+1)!]^5}.
-
 Reciprocal Fibonacci constant
 -------------------------------------------------------------------------------
 

doc/source/references.rst

@@ -295,10 +295,12 @@ References
 
 .. [Zip1985] \R. Zippel. "Simplification of expressions involving radicals". Journal of Symbolic Computation (1985) 1, 189-210. https://doi.org/10.1016/S0747-7171(85)80014-6
 
+.. [Zun2023] \J. Zuniga, "Catalan's constant fast convergent series", https://mathoverflow.net/q/424055
+
 .. [vHP2012] \M. van Hoeij and V. Pal. "Isomorphisms of algebraic number fields". Journal de Théorie des Nombres de Bordeaux, Vol. 24, No. 2 (2012), pp. 293-305. https://doi.org/10.2307/43973105
 
 .. [vdH1995] \J. van der Hoeven, "Automatic numerical expansions". Proc. of the conference Real numbers and computers (1995), 261-274. https://www.texmacs.org/joris/ane/ane-abs.html
 
 .. [vdH2006] \J. van der Hoeven, "Computations with effective real numbers". Theoretical Computer Science, Volume 351, Issue 1, 14 February 2006, Pages 52-60. https://doi.org/10.1016/j.tcs.2005.09.060
 
-All referenced works: [AbbottBronsteinMulders1999]_, [Apostol1997]_, [Ari2011]_, [Ari2012]_, [Arn2010]_, [Arn2012]_, [ArnoldMonagan2011]_, [BBC1997]_, [BBC2000]_, [BBK2014]_, [BD1992]_, [BF2020]_, [BFSS2006]_, [BJ2013]_, [BM1980]_, [BZ1992]_, [BZ2011]_, [BaiWag1980]_, [BerTas2010]_, [Blo2009]_, [Bodrato2010]_, [Boe2020]_, [Bog2012]_, [Bor1987]_, [Bor2000]_, [Bre1978]_, [Bre1979]_, [Bre2010]_, [BrentKung1978]_, [BuhlerCrandallSompolski1992]_, [CGHJK1996]_, [CP2005]_, [Car1995]_, [Car2004]_, [Chen2003]_, [Cho1999]_, [Coh1996]_, [Coh2000]_, [Col1971]_, [CraPom2005]_, [DYF1999]_, [DelegliseNicolasZimmermann2009]_, [DomKanTro1987]_, [Dup2006]_, [Dus1999]_, [EHJ2016]_, [EM2004]_, [Fie2007]_, [FieHof2014]_, [Fil1992]_, [GCL1992]_, [GG2003]_, [GS2003]_, [GVL1996]_, [Gas2018]_, [Gos1974]_, [GowWag2008]_, [GraMol2010]_, [HM2017]_, [HS1967]_, [HZ2004]_, [HanZim2004]_, [Har2010]_, [Har2012]_, [Har2015]_, [Har2018]_, [Hart2010]_, [Hen1956]_, [Hoe2001]_, [Hoe2009]_, [Hor1972]_, [Iliopoulos1989]_, [JB2018]_, [JM2018]_, [JR1999]_, [Joh2012]_, [Joh2013]_, [Joh2014a]_, [Joh2014b]_, [Joh2014c]_, [Joh2015]_, [Joh2015b]_, [Joh2016]_, [Joh2017]_, [Joh2017a]_, [Joh2017b]_, [Joh2018a]_, [Joh2018b]_, [JvdP2002]_, [Kahan1991]_, [KanBac1979]_, [Kar1998]_, [Knu1997]_, [Kob2010]_, [Kri2013]_, [Leh1970]_, [LukPatWil1996]_, [MP2006]_, [MPFR2012]_, [MasRob1996]_, [Mic2007]_, [Miy2010]_, [Mos1971]_, [Mul2000]_, [NIST2012]_, [NakTurWil1997]_, [Olv1997]_, [PP2010]_, [PS1973]_, [PS1991]_, [Paterson1973]_, [PernetStein2010]_, [Pet1999]_, [Pla2011]_, [Pla2017]_, [RF1994]_, [Rad1973]_, [Rademacher1937]_, [Ric1992]_, [Ric1995]_, [Ric1997]_, [Ric2007]_, [Ric2009]_, [RosSch1962]_, [Rum2010]_, [Smi2001]_, [SorWeb2016]_, [Ste2002]_, [Ste2010]_, [Stehle2010]_, [Stein2007]_, [Sut2007]_, [StoMul1998]_, [Str1997]_, [Str2012]_, [Tak2000]_, [ThullYap1990]_, [Tre2008]_, [Tru2011]_, [Tru2014]_, [Tur1953]_, [Villard2007]_, [WaktinsZeitlin1993]_, [Wei2000]_, [Whiteman1956]_, [Zip1985]_, [vHP2012]_, [vdH1995]_, [vdH2006]_
+All referenced works: [AbbottBronsteinMulders1999]_, [Apostol1997]_, [Ari2011]_, [Ari2012]_, [Arn2010]_, [Arn2012]_, [ArnoldMonagan2011]_, [BBC1997]_, [BBC2000]_, [BBK2014]_, [BD1992]_, [BF2020]_, [BFSS2006]_, [BJ2013]_, [BM1980]_, [BZ1992]_, [BZ2011]_, [BaiWag1980]_, [BerTas2010]_, [Blo2009]_, [Bodrato2010]_, [Boe2020]_, [Bog2012]_, [Bor1987]_, [Bor2000]_, [Bre1978]_, [Bre1979]_, [Bre2010]_, [BrentKung1978]_, [BuhlerCrandallSompolski1992]_, [CGHJK1996]_, [CP2005]_, [Car1995]_, [Car2004]_, [Chen2003]_, [Cho1999]_, [Coh1996]_, [Coh2000]_, [Col1971]_, [CraPom2005]_, [DYF1999]_, [DelegliseNicolasZimmermann2009]_, [DomKanTro1987]_, [Dup2006]_, [Dus1999]_, [EHJ2016]_, [EM2004]_, [Fie2007]_, [FieHof2014]_, [Fil1992]_, [GCL1992]_, [GG2003]_, [GS2003]_, [GVL1996]_, [Gas2018]_, [Gos1974]_, [GowWag2008]_, [GraMol2010]_, [HM2017]_, [HS1967]_, [HZ2004]_, [HanZim2004]_, [Har2010]_, [Har2012]_, [Har2015]_, [Har2018]_, [Hart2010]_, [Hen1956]_, [Hoe2001]_, [Hoe2009]_, [Hor1972]_, [Iliopoulos1989]_, [JB2018]_, [JM2018]_, [JR1999]_, [Joh2012]_, [Joh2013]_, [Joh2014a]_, [Joh2014b]_, [Joh2014c]_, [Joh2015]_, [Joh2015b]_, [Joh2016]_, [Joh2017]_, [Joh2017a]_, [Joh2017b]_, [Joh2018a]_, [Joh2018b]_, [JvdP2002]_, [Kahan1991]_, [KanBac1979]_, [Kar1998]_, [Knu1997]_, [Kob2010]_, [Kri2013]_, [Leh1970]_, [LukPatWil1996]_, [MP2006]_, [MPFR2012]_, [MasRob1996]_, [Mic2007]_, [Miy2010]_, [Mos1971]_, [Mul2000]_, [NIST2012]_, [NakTurWil1997]_, [Olv1997]_, [PP2010]_, [PS1973]_, [PS1991]_, [Paterson1973]_, [PernetStein2010]_, [Pet1999]_, [Pla2011]_, [Pla2017]_, [RF1994]_, [Rad1973]_, [Rademacher1937]_, [Ric1992]_, [Ric1995]_, [Ric1997]_, [Ric2007]_, [Ric2009]_, [RosSch1962]_, [Rum2010]_, [Smi2001]_, [SorWeb2016]_, [Ste2002]_, [Ste2010]_, [Stehle2010]_, [Stein2007]_, [Sut2007]_, [StoMul1998]_, [Str1997]_, [Str2012]_, [Tak2000]_, [ThullYap1990]_, [Tre2008]_, [Tru2011]_, [Tru2014]_, [Tur1953]_, [Villard2007]_, [WaktinsZeitlin1993]_, [Wei2000]_, [Whiteman1956]_, [Zip1985]_, [Zun2023]_, [vHP2012]_, [vdH1995]_, [vdH2006]_

src/arb/const_apery.c

@@ -21,14 +21,15 @@ arb_const_apery_eval(arb_t s, slong prec)
     arb_init(t);
     hypgeom_init(series);
 
-    fmpz_poly_set_str(series->A, "3  77 250 205");
+    fmpz_poly_set_str(series->A, "12  2539221463380 55800350430619 543035311757517 3094818610007883 11495362203169095 29262452925092202 52160137207884216 65180430489299744 56019281176488240 31567339140195744 10506477648460032 1565994397644288");
     fmpz_poly_set_str(series->B, "1  1");
-    fmpz_poly_set_str(series->P, "6  0 0 0 0 0 -1");
-    fmpz_poly_set_str(series->Q, "6  32 320 1280 2560 2560 1024");
+    fmpz_poly_set_str(series->P, "15  0 0 0 0 0 -30 691 -6781 37374 -127976 283232 -406224 364896 -186624 41472");
+    fmpz_poly_set_str(series->Q, "15  -44008272000 -2334151436400 -53522442803340 -703273183134030 -5931859745397870 -34140867105175650 -139058868850409430 -409481300311614720 -880500176512163280 -1382139595517666400 -1565294958171053280 -1244539247650560000 -658690593528960000 -208277254886400000 -29753893555200000");
 
     prec += 4 + FLINT_CLOG2(prec);
     arb_hypgeom_infsum(s, t, series, prec, prec);
-    arb_mul_ui(t ,t, 64, prec);
+    arb_mul_ui(t, t, 1031443875, prec);
+    arb_mul_2exp_si(t, t, 11);
     arb_div(s, s, t, prec);
 
     hypgeom_clear(series);

src/arb/const_catalan.c

@@ -22,14 +22,14 @@ arb_const_catalan_eval(arb_t s, slong prec)
     arb_init(t);
     hypgeom_init(series);
 
-    fmpz_poly_set_str(series->P, "5  32 160 288 224 64");
-    fmpz_poly_set_str(series->Q, "5  225 3240 14904 23328 11664");
-    fmpz_poly_set_str(series->A, "3  411 976 580");
-    fmpz_poly_set_str(series->B, "5  1 5 9 7 2");
+    fmpz_poly_set_str(series->A, "7  1999553 21620948 94165776 211938912 260619984 166411584 43203456");
+    fmpz_poly_set_str(series->B, "1  1");
+    fmpz_poly_set_str(series->P, "9  0 0 0 1280 -17536 86400 -195840 207360 -82944");
+    fmpz_poly_set_str(series->Q, "9  363825 12034680 150240200 918651040 3101725520 6073920000 6863040000 4147200000 1036800000");
 
-    prec += FLINT_CLOG2(prec);
+    prec += 4 + FLINT_CLOG2(prec);
     arb_hypgeom_infsum(s, t, series, prec, prec);
-    arb_mul_ui(t, t, 450, prec);
+    arb_mul_ui(t, t, 2182950, prec);
     arb_div(s, s, t, prec);
 
     hypgeom_clear(series);