$ \newcommand{\setC}{\mathbb{C}} \newcommand{\setN}{\mathbb{N}} \newcommand{\setQ}{\mathbb{Q}} \newcommand{\setR}{\mathbb{R}} \newcommand{\setZ}{\mathbb{Z}} \newcommand{\SL}[2]{\mathrm{SL}_{#1}(\mathbb{#2})} \newcommand{\mt}[4]{\begin{bmatrix}#1\\#3\end{bmatrix}} \newcommand{\abs}[1]{\left\lvert#1\right\rvert} $ Via jupytext this file can be shown as a jupyter notebook.
Demo for derivation of $\frac{1}{\pi}$ formula for $N=29$¶
This file computes a Ramanujan-Sato series for $\frac{1}{\pi}$ of Gosper type and is connected to an article (RISC report 25-01) by Ralf Hemmecke, Peter Paule, and Cristian-Silviu Radu.
The $\frac{1}{\pi}$ formula used by the Gosper in 1985 to compute a world record of $17 \cdot 10^6$ digits can be found by setting $N=29$.
Ramanujan-Sato_series: \begin{align*} \frac{1}{\pi} &= \frac{2\sqrt{2}}{99^2} \sum_{n=0}^\infty \frac{(4n)!}{n!^4} \frac{26390n+1103}{396^{4n}} \end{align*}
The version that Gosper used is constructed with $N=29$ and a modular form \begin{align*} g(\tau)=(1+\lambda(\tau))\theta_3(\tau) \end{align*} of weight 2 for $\Gamma(2)$ and a modular function for $\Gamma(2)$ \begin{align*} h(\tau) &= \frac{\lambda(\tau) \cdot (1-\lambda(\tau))^2} {16 \cdot (1+\lambda(\tau))^4}. \end{align*}
Note that in this file we only demonstrate how to compute the numbers $\tau_N$, $h(\tau_N)$, $p_1(h(\tau_N))$, and $p_2(h(\tau_N))$. We assume that the factor $\frac{(4n)!}{n!}$ that appears in all $\frac{1}{\pi}$ formulas of the Gosper family has already been determined.
)cd ..
)read input/jfricas-test-support.input )quiet
The current FriCAS default directory is /home/hemmecke/backup/git/qeta All user variables and function definitions have been cleared. All )browse facility databases have been cleared. Internally cached functions and constructors have been cleared. )clear completely is finished. The current FriCAS default directory is /home/hemmecke/backup/git/qeta/tmp
AN ==> AlgebraicNumber
moebiusAN(gamma, x) ==> moebiusTransform(gamma, x)$QEtaOneOverPi(AN)
outputSpacing 0; -- no spacing in the output
outputGeneral 20; -- only output 20 digits
)set output formatted off
)set output algebra on
Computation of $\frac{1}{\pi}$ formula used by Gosper (condensed version)¶
In the QEta package we chose to collect all the steps for computing the respective coefficients of a Ramanujan-Sato $\frac{1}{\pi}$ formula by a single function from the Sato triple and the value $z_N=h(\tau_N)$.
Since the computation of the modular polynomials usually takes a long time (about 1 hour for $N=29$), we compute those polynomials only once and store them on disk.
)read projectdir.input )quiet
basedir := PROJECTDIR "/data/oneoverpi/gosper";
Let $N=29$ and let us define $W=\left\langle\Gamma(2)\cup\{g_1,g_2\}\right\rangle$ with \begin{alignat*}{2} g_1 &= \mt{1}{0}{1}{1}, &\qquad g_2 &= \mt{0}{-2}{1}{0}. \end{alignat*}
nn := 29;
g1:=matrix [[1,0],[1,1]],g2:=matrix [[0,-2],[1,0]]
In order to compute the $\frac{1}{\pi}$ formula used by Gosper, we need to find a Sato triple for $N=29$. We can actually compute several such triples and then select the one that corresponds to the Sato-Ramanujan series used by Gosper.
satotris := satoTriples(nn,2);
[x.fgamma::OF*x.ftau = moebiusAN(x.fgamma,x.ftau) for x in entries satotris]
With other parameters we also find the triple that we use in the paper.
satotri := gammaTauQCandidatesWurm(29, 59, -2).1;
gammaN := satotri.fgamma;
tauN := satotri.ftau;
gammaN::OF*tauN = moebiusAN(gammaN,tauN)
The $S$ and $T$ matrices are connected to transforming $\tau_N$ into $\tau_\infty= T \tau_N$, $S=T^{-1}$.
[gammafd := toFundamentalDomain(tauN,2), g2*g1^(-29)]
[g2^(-2), tmat := g2^(-2)*gammafd, smat:=tmat^(-1)]
tauInf := moebiusAN(tmat,tauN)
With these values, we can compute $z_N:=h(\tau_N)$ in a purely algebraic manner.
zN := 1/gosperIotaAt(nn,1,tauInf,basedir)
With that value, we compute easily the other coefficients of the $\frac{1}{\pi}$ formula, details below.
[aa,bb,ff] := gosperABFCoefficients(nn,gammaN,tauN,zN,basedir)
gosperOneOverPiFormula(nn,gammaN,tauN,zN,basedir)
Computation of $\frac{1}{\pi}$ formula used by Gosper (long version)¶
In this section, we explain in more detail what the functions from the previous sections actually do.
Setup¶
We introduce a few macros and functions for convenience.
C ==> QQ
MZZ ==> Matrix ZZ
SUP ==> SparseUnivariatePolynomial
RFZ ==> Fraction SUP ZZ -- rational function over C
SYS ==> systemCommand $ MoreSystemCommands
Q1PIQ ==> QEtaOneOverPiQ
Q1PIL ==> QEtaOneOverPiLambda
We collect expansion data of a function involving the modular lambda function and Jacobi theta functions (constants, i.e. $z=0$). This is used to find the cusps that are actually relevant for a samba computation.
CUSPDELTAMATRICES ==> _
Record(fdelta: PP, fgamma: SL2Z, fred: SL2Z, ftriang: MZZ)
CUSPEXPANSION C ==> Record(fmats: CUSPDELTAMATRICES, fex: L1 C)
CUSPEXPANSIONS C ==> XHashTable(SL2Z, CUSPEXPANSION C)
FUN C ==> (CUSPDELTAMATRICES, QQ) -> L1 C
$g(\tau)$, $h(\tau)$, and $\iota(\tau)$ expansion at $\infty$¶
We use the inverse of the following modular function (instead of the Klein $j$-invariant as in the Chudnovsky case).
In fact, we use $h(\tau)=\mathrm{hrf}(\lambda(\tau))$ where $\lambda$ is the modular $\lambda$-function.
We always expand the series in terms of $q=\exp(\pi i \tau)$ and precompute at least 6 terms of the series expansions for display.
)set stream calc 6
We define
$\iota(\tau) := \frac{1}{h(\tau)}$
and take the inverse of the hrf expression for $h$ namely
hrf := ((x * (1-x)^2) / (16 * (1+x)^4))::RFZ;
iotarf := 1/hrf;
ml := modularLambdaPower(1, 1)$CachedQFunctions(C, L1 C)
h := evalAtModularLambda(hrf)$QEtaModularLambdaTools(C)
Note that FriCAS uses a lazy Laurent series implementation that computes a coefficient at the time it is needed.
coefficient(h,50)
We define the function $g(\tau)=(1+\lambda(\tau))\theta_3(\tau)^4$. $g$ is a modular form for $\Gamma(2)$ of weight 2.
th3p4 := thetaPower4(3)$CachedQFunctions(C, L1 C)
g := (1 + ml) * th3p4
Guessing $Y$ with $g(\tau)=Y(h(\tau))$¶
We get $Y$ by computing $Y := \tilde{g}\circ \tilde{h}^{-1}$ where $\tilde{g}(q)=g(\tau)$ and $\tilde{h}(q)=h(\tau)$.
yseries := g(revert h)
We can roughtly determine its convergent radius by the root criterion. $1/r=\limsup_{n\to\infty}\sqrt[n]{\abs{y_n}}$. It seems to have a very small convergent radius.
1/(nthRoot(coefficient(yseries,200)::Float,200))
For making it precise, we need a recursion formula for the coefficients. It can be found algorithmically and even proven algorithmically. We only show here the finding (guessing) process.
)clear prop n
)expose RECOP
RecurrenceOperator is now explicitly exposed in frame initial
ycoeffs := [coefficient(yseries, n) for n in 0..22]
gr := guessRec(ycoeffs,functionName=='co).1
co := operator 'co;
eqy := getEq(gr);
gh := guessHolo(ycoeffs).1
res := guess(ycoeffs).1
From the last result, we can easily extract $\frac{(4n)!}{(n!)^4}$ by replacing the running variable by $k$ and forming the product from $k=1$ to $n$.
a:=argument(kernels(res).1);
factorFraction(eval(a.1,a.2=k-1)::Fraction(Polynomial ZZ))
Here we assume that this formula is already given. In fact, for our chosen $g$ and $h$, we get \begin{align} Y(x) &= \sum_{n=0}^\infty \frac{(4n)!}{(n!)^4} x^n. \end{align}
gosperFactor(n) ==> factorial(4*n)/((factorial(n)^4))
assertEquals(removeDuplicates [eval(gr,n=x) - coefficient(yseries,x) for x in 0..200],[0])
assertEquals(removeDuplicates [gosperFactor n - coefficient(yseries,n) for n in 0..200],[0])
assertTrue(zero? normalize eval(eqy,[co(n+1)=gosperFactor(n+1), co(n)=gosperFactor(n)]))
It is clear that the expansion for $g(\tau)$ and $h(\tau)$ converge for $\abs{q}<\frac{1}{2}$ where $q=\exp(\pi i \tau)$. It is easy to conclude that $Y(x)$ converges for all $x\in\setC$ with $\abs{x}<\frac{1}{256}$. Together, we can assure that $\tilde{g}(q)=Y(\tilde{h}(q))$ holds if $\abs{q}$ and $\tilde{h}(q)=h(\tau)$ are less than $\frac{1}{256}$.
Computation of $h(\tau_N)$¶
In fact, we compute $\iota(\tau_N)=\frac{1}{h(\tau_N)}$.
Let us recall the values from the Sato triple.
nn, gammaN, tauN, moebiusAN(gammaN,tauN), nn*tauN,smat,moebiusAN(tmat,tauN)
Since $\iota$ is invariant under group actions of $W=\left\langle\Gamma(2)\cup\{g_1,g_2\}\right\rangle$, we can instead evaluate $\iota$ at $\sqrt{-58}$. However, we do not do any numerical evaluation, but rather determine the value $\iota(\tau_N)$ as a root of a univariate polynomial.
Cusps of $W'$¶
Since we aim at the computation of the modular polynomial of $\iota(\tau)$ and $\iota(N\tau)$, we must consider the group \begin{align*} W' &= \left\{\mt{a}{b}{c}{d} \in W: \mt{a}{bN}{c/N}{d}\in W\right\}. \end{align*}
We have not explicitly determined the cusps for the group $W´$. Instead we note that $\iota(\tau)$ and $\iota(N\tau)$ are modular functions for $\Gamma(2N)$ and we can compute with the $q$-expansions at the cusps of $\Gamma(2N)$.
Unfortunately, there are 1260 cusps.
# cusps()$CongruenceSubgroupGammaN(58)
$\iota(\tau)$ and $\iota(N\tau)$, expansion at the reduced cusps¶
Fortunately, we can reduce the number of expansions drastically. Suppose that we have two transformations $\gamma_1,\gamma_2\in\SL2Z$. Furthermore, assume that $\iota(\gamma_i\tau)$ corresponds to the expansions $f_{A,i}(x)\in \setC[[x]]$, whereas $\iota(N \gamma_i\tau)$ corresponds to the expansions $f_{B,i}(x)\in \setC[[x]]$ in the canonical variable $x=\exp\left(\frac{\pi i \tau}{2N}\right)$ of $\Gamma(2N)$.
Clearly, if $f_{A,1}(x)=f_{A,2}(x)$ and $f_{B,1}(x)=f_{B,2}(x)$, then there is no need to compute with the expansions at both cusps, i.e. one of the cusps can be avoided. The same holds for the case $f_{A,1}(x^u)=f_{A,2}(\xi x^v)$ and $f_{B,1}(x^u)=f_{B,2}(\xi x^v)$ for two natural numbers $u, v$ and a root of unity $\xi$.
By such a reduction, we can reduce for $\iota(\gamma_i\tau)$ and $\iota(N \gamma_i\tau)$ to just the cusps $0$ and $\infty$ of $\Gamma(58)$.
We do the expansion in terms of $q=\exp\left(\frac{\pi i \tau}{2N}\right)$, since the cusp width of $\Gamma(2N)$ is $2N$.
seriesdir := basedir "/" string(nn);
modpoldir := basedir "/" string(nn);
fnbase := seriesdir "/hl";
SYS("system mkdir -p " seriesdir);
nn2: PP := 2*nn;
oo := infinity()$Cusp;
spitzs := cons(oo, remove(oo, cusps()$GAMMAN(nn2)));
cmats := [cuspDeltaMatrices(nn, nn2, x)$QEtaModularLambdaAuxiliaryPackage for x in spitzs];
# cmats
To demonstrate, we show here the expansions of $\iota(\tau)$ and $\iota(N\tau)$ at the cusps $\infty$, and $\frac{5}{58}$.
Note that cmats.i describes the cups, the corresponding
transformation matrix $\gamma$ and a splitting
$29\gamma= \gamma_R \cdot \gamma_T$ where $\gamma_R$ corresponds
to the matrix in field red of the data structure and
$\gamma_T$ corresponds to the triangular matrix in field
triang, i.e. we have
$\iota(29\gamma\tau)=\iota(\gamma_R \tau')$ for
$\tau':= \gamma_T \tau$. (see below).
cm := cmats.1
)set stream calc 200
iota1Inf := exIota1(cm, nn2)
)set stream calc 6000
iotaNInf := exIotaN(cm, nn2)
cm := cmats.25
)set stream calc 200
iota1_5_58 := exIota1(cm, nn2)
)set stream calc 6000
iotaN_5_58 := exIotaN(cm, nn2)
Since $\iota$ is a rational function in the modular $\lambda$
function, and the gamma fields of cmats.i for $i=1$ and
$i=25$ agree modulo 2, the expansions of
$\iota(\tau)$ at $\infty$ must agree with the
expansion of $\iota$ at the cusp $\frac{5}{58}$ according to
fungrim:099301. The same is
true for $\iota(29\tau)$, since the corresponding $\gamma_R$
matrices are congruent modulo 2 and the respective $\gamma_T$
matrices are equal. That means any computation done with
$\iota$ at $\infty$ is exactly the same to the computation
with the expansion at $\frac{5}{58}$. Thus,
the cusp $\frac{5}{58}$ can be excluded.
In a similar fashion we exclude many other cusps, so that for the case $\iota(\tau)$ and $\iota(29\tau)$ we only have to expand at the cusps $\infty$ and $0$.
xiord := 2 * lcm [(x.ftriang)(2,2) for x in cmats];
EXTENDEDCOEFFICIENTRING(C, xiord, CX, xi);
scred ==> reduceCuspsByScaling $ QETASCAL(C, xiord, CX, xi, GAMMAN nn2)
cmatsreduced1N := scred(exIota1, exIotaN, cmats)
assertEquals([cusp(cm.fgamma) for cm in cmatsreduced1N], [cusp(0,1),oo])
Note that we have a safety net here. After having derived a polynomial $p(x,y)$ just from appying MultiSamba to Laurent series expansions at two cusps $\infty$ and $0$, we can check easily whether it is such that $p(\iota(\tau),\iota(N\tau))=0$ by plugging in each of the 1260 cusp expansions and verify that the principle part (including the constant) of the Laurent series is 0. This is a finite check with elements in $\setQ[\xi]$ for a certain root of unity $\xi$ with $\xi^{58}=1$.
rcmats := reverse cmatsreduced1N; -- only cusps \infty and 0
cm := rcmats.1
)set stream calc 200
iota1Inf := exIota1(cm, nn2)
)set stream calc 6000
iotaNInf := exIotaN(cm, nn2)
cm := rcmats.2
)set stream calc 400
iota10 := exIota1(cm, nn2)
)set stream calc 12
iotaN0 := exIotaN(cm, nn2)
Modular polynomial of $\iota(\tau)$ and $\iota(N\tau)$¶
fIota1: FUN QQ := exIota1 $ Q1PIL(QQ);
fIotaN: FUN QQ := exIotaN $ Q1PIL(QQ);
The following computation takes less than 20 seconds. Nevertheless, we store the resulting polynomial $\phi_N(x,y)$ on disk and then read it back for other Sato triples for $N=29$.
phiN := gosperModularPolynomial("1N", fIota1, fIotaN, [x,y], nn2, rcmats, modpoldir);
The function gosperModularPolynomial works in the same way as is
demonstrated in
ModularEquationJLambda
to find the Relation between the Klein $j$ invariant
and the modular $\lambda$ function.
Rational bounds for $\iota(\tau_N)$¶
In order to find rational upper and lower bounds for the value of $\iota(\tau_N)=\iota\!\left(\frac{i\sqrt{58}+58}{1711}\right)=\iota(i\sqrt{58})$, we determine first a rational interval for $\exp(\pi i(i\sqrt{58}))=\exp(-\pi\sqrt{58})$. We use the facts $\frac{3141}{1000}<\pi<\frac{3142}{1000}$ and $\frac{7615}{1000}<\sqrt{58}<\frac{7616}{1000}$ and expand the Taylor series for $\exp$ up to 167 and 168 terms for the lower and upper bound. To make the numerator and denominator of the rational numbers a bit smaller, we take the 4th and 5th contined fraction convergent.
qi := expPiIInterval(sqrt(-58))
Next, we determine from the above bounds and a truncation of the theta series for $\lambda(\tau)=\frac{\theta_2(\tau)^4}{\theta_3(\tau)^4}$ a rational interval such that $\lambda_l < \lambda(i\sqrt(58)) < \lambda_u$. Finally, we determine a rational interval for $\iota_l < \iota(i\sqrt(58)) < \iota_u$.
iotai := gosperIotaInterval(1, tauInf)
[inf(iotai)::Float, sup(iotai)::Float]
Exact evaluation¶
We can use the fact that $\iota(\tau_N)=\iota(N\tau_N)$ in the modular polynomial $\phi_N(x,y)$ for which $\phi_N(\iota(\tau),\iota(N\tau))=0$. Then, clearly, $\iota(\tau_N)$ must be a root of $\phi_N(x,x)$.
We determine the actual value, by considering every factor (without its multiplicity) of $\phi_N(x,x)$ and check that in the interval from above there is exactly one real root.
phiNxx := eval(phiN, 'y='x)::UnivariatePolynomial('x,ZZ);
psiN := squareFreePart phiNxx;
psiNfactors := [x.factor for x in factors factor psiN]
We find that only the first factor contains a root in the rational interval $(\iota_l, \iota_u)$ from above.
nroots := [countRoots(p,iotai) for p in psiNfactors]
Since $\iota_l < \iota(\tau_N) < \iota_u$, $\iota(\tau_N)$ is a root of $\phi_N(x,x)$ and there is only one factor of $\phi_N(x,x)$ that has a root in the interval $[\iota_l, \iota_u]$, $\iota(\tau_N)$ must be equal to the root of this factor.
psiNfactors.1
rr := radicalRoots(psiNfactors.1,'x)
assertEquals(rr.1, 1/zN)
Modular polynomial for $L(\tau)$¶
We have (with $w=2$ and $q=\exp(\pi i \tau)$) $$ A(\tau) := \frac{w}{2 \pi i} \, \frac{h'(\tau)}{g(\tau)} = \frac{q \tilde{h}'(q)}{\tilde{g}(q)}. $$
)set stream calc 6
qq := monomial(1,1)$L1(C)
sera := qq*D(h)/g
This series can also be obtained via the expansion of a rational function of $\lambda(\tau)$. \begin{align} h'(\tau) &= \frac{d}{d\tau}r(\lambda(\tau)) = \pi i \, r'(\lambda(\tau)) \lambda'(\tau) = \pi i \, r'(\lambda(\tau)) \, \lambda(\tau) \, \theta_4(0,\tau)^4 \\ A(\tau) &= \frac{w}{2 \pi i} \frac{Z'(\tau)}{F(\tau)} = \frac{w r'(\lambda(\tau)) \, \lambda(\tau) \, \theta_4(0,\tau)^4} {32 (1+\lambda(\tau)) \theta_3(0,\tau)^4} = \frac{w \lambda(\tau) (1-\lambda(\tau))^2 (\lambda(\tau)^2-6\lambda(\tau)+1)} {32 (1+\lambda(\tau))^6} \end{align} if we take fungrim:04d3a6 into account, i.e. \begin{align} 1-\lambda(\tau) &= \frac{\theta_4(0,\tau)^4}{\theta_3(0,\tau)^4}. \end{align}
Note that in the following the question mark is the name of the
variable in $\setQ[?]$.
The series sera and the series that results from evaluating
the rational function arf at the the $q$-series expansion
of $\lambda(\tau)$ agree in the first 100 coefficients.
xx := monomial(1,1)$SUP(ZZ)::RFZ;
arf := xx * D(hrf) * (1-xx)/(1+xx); factorFraction arf
ser := sera - evalAtModularLambda(arf)$QEtaModularLambdaTools(C);
assertEquals(removeDuplicates [coefficient(ser,n) for n in 0..100], [0])
So we obtain $L(\tau)=\rho(\lambda(\tau))$ with $\rho$ given by the
variable lrf.
lrf := iotarf^2*arf;
factorFraction(lrf)
We can easily check that $p_L(x,y)=y^2 + x (256 - x)$ vanishes, if
we replace $x$ by iotarf, i.e. $\frac{1}{r}$, and $y$ by lrf,
i.e. $\rho$. That means that $p_L(\iota(\tau),L(\tau))=0$.
pL := y^2 + x * (256 - x)
eval(pL, [x=iotarf,y=lrf])
Another way to find that polynomial is by computing the modular polynomial for $\iota(\tau)$ and $L(\tau)$ as demonstrated in ModularEquationJLambda.
oo := infinity()$Cusp;
spitz2s := cons(oo, remove(oo, cusps()$GAMMAN(2)));
cmat2s := [cuspDeltaMatrices(1, 1, x)$QEtaModularLambdaAuxiliaryPackage for x in spitz2s];
[f.fgamma for f in cmat2s]
assertEquals(#cmat2s, 3)
modpol1L := gosperModularPolynomial("1L", exIota1, expandL, ['x,'y], 2, cmat2s, modpoldir)$Q1PIQ
Let us check that the actual expansions at the cusps also fulfil the polynomial.
check1L(n) ==> (_
cm := cmat2s.n; gamma:= cm.fgamma;_
iota := exIota1(cm,1);_
ll := expandL(cm,1);_
tracePrint("gamma", gamma);_
tracePrint(" iota", iota);_
tracePrint(" L", ll);_
eval(modpol1L, [x=iota, y=ll]))
check1L(1)
check1L(2)
check1L(3)
-- gamma:=matrix[[1, 0], [0, 1]] -- iota:=q^(-1)+104+4372*q+96256*q^2+1240002*q^3+10698752*q^4+74428120*q^5+O(q^6) -- L:=q^(-1)-24-3820*q-100352*q^2-1217598*q^3-10797056*q^4-74060072*q^5+O(q^6)
-- gamma:=matrix[[0, -1], [1, 0]] -- iota:=q^(-2)+104+4372*q^2+96256*q^4+O(q^5) -- L:=-q^(-2)+24+3820*q^2+100352*q^4+O(q^5)
-- gamma:=matrix[[1, -1], [1, 0]] -- iota:=-q^(-1)+104-4372*q+96256*q^2-1240002*q^3+10698752*q^4-74428120*q^5+O(q^6) -- L:=-q^(-1)-24+3820*q-100352*q^2+1217598*q^3-10797056*q^4+74060072*q^5+O(q^6)
In this polynomial $x$ corresponds to $\frac{1}{z_N}$, thus we replace $x$ by $\frac{1}{z_N}$ and compute the candidates for $A(\tau_N)$.
ataucands := [zN^2*(x::AN) for x in radicalRoots(eval(modpol1L, 'x=1/zN), 'y)]
Modular polynomials for of $K(\tau)$¶
Let $B(\tau) := \frac{H(\tau)}{g(\tau)}$ and $w=2$.
For the transformation of $B$ by some $\gamma=\mt{a}{b}{c}{d}\in\SL2Z$, we find $\bar{\gamma}:=\mt{\bar{a}}{\bar{b}}{\bar{c}}{\bar{d}}$ and $\bar{\tau}:=\frac{\bar{u}\tau+\bar{v}}{\bar{w}}$ such that $N\gamma=\mt{aN}{bN}{c}{d}=\mt{\bar{a}}{\bar{b}}{\bar{c}}{\bar{d}} \mt{\bar{u}}{\bar{v}}{0}{\bar{w}}$, where $\bar{a}=\frac{Na}{\gcd(N,c)}$, $\bar{c}=\frac{c}{\gcd(N,c)}$, and $\bar{b}$ and $\bar{d}$ are chosen such that $\bar{a}\bar{d}-\bar{b}\bar{c}=1$.
Then \begin{align*} B(\gamma\tau) &= \frac{\frac{w}{2 \pi i} \left( \frac{g'( \gamma \tau)}{g( \gamma \tau)} -N\frac{g'(\bar\gamma\bar\tau)}{g(\bar\gamma\bar\tau)} \right)} {g(\gamma\tau)} = \frac{ \phi(\gamma,\tau) -\frac{\gcd(N,c)^2}{N} \, \phi(\bar \gamma, \bar\tau) }{ (-1)^{ad+1} (1+r_{\gamma}(\lambda(\tau))) \theta_{S_3(\gamma)}(\tau)^4} \end{align*} where $S_3(\gamma)$ is given by fungrim:9bda2f and \begin{align*} \phi(\gamma,\tau) &= \frac{r_{\gamma}'(\lambda(\tau))\, \frac{w}{2 \pi i}\lambda'(\tau)}{1+r_\gamma(\lambda(\tau))} + 4 \frac{\frac{w}{2 \pi i}\theta_{S_3(\gamma)}'(\tau)} {\theta_{S_3(\gamma)}(\tau)}, \end{align*} follows from the transformation formulas of the Jacobi theta functions and the modular lambda function.
The expansion of
$K(\tau) = \frac{29}{24}\iota(\tau) \iota(N\tau) B(\tau)$ is
given by the function expandK.
spitzs := cons(oo, remove(oo, cusps()$GAMMAN(nn2)));
cmats := [cuspDeltaMatrices(nn, nn2, s)$QEtaModularLambdaAuxiliaryPackage for s in spitzs];
xiord := 2 * lcm [(cm.ftriang)(2,2) for cm in cmats]
EXTENDEDCOEFFICIENTRING(C, xiord, CX, xi);
fexpandK := expandK $ Q1PIL(QQ);
scred ==> reduceCuspsByScaling $ QETASCAL(QQ, xiord, CX, xi, GAMMAN nn2)
cmatsreduced1K := scred(exIota1, expandK, cmats);
cmatsreducedNK := scred(exIotaN, expandK, cmats);
assertEquals(cmatsreduced1K, cmatsreducedNK)
rcmats := cmatsreduced1K;
#rcmats
jtN := (simplifyRadicals(1/zN)$AlgebraicNumberTool)::Pol(AN)
modpol1K := gosperModularPolynomial("1K", fexpandK, fIota1, ['y,'x], nn2, rcmats, modpoldir);
modpolNK := gosperModularPolynomial("NK", fexpandK, fIotaN, ['y,'x], nn2, rcmats, modpoldir);
For the actual computation (for N=29), we did not use the
function gosperModularPolynomial in this way.
Here it is just loading an already
precomputed polynomial stored on disk.
The actual computation with rational coefficients of the
Laurent series lead to huge numerators and denominators
and would not have finished in reasonable time.
By experimentation with smaller $N$ we guessed a scaling factor of $\frac{29}{24}$ that we used in the definition of $K(\tau)$. That leads to a modular polynomial with integer coefficients. To reduce coefficient growth, we computed modulo the prime $10^{696}+61$ and eventually checked that the expansions at the cusps (with rational coefficients) make that polynomial (interpreted as a polynomial over $\setZ$) vanish.
checkK(n,m, mp,fexpandx) ==> (_
cm := rcmats.n; gamma:= cm.fgamma;_
iota := fexpandx(cm,m);_
kk := fexpandK(cm,m);_
tracePrint("gamma", gamma);_
tracePrint(" iota", iota);_
tracePrint(" K", kk);_
ser := eval(mp, [x=iota, y=kk])::L1(C);_
grdser := -degree(ser);_
tracePrint(" grd", grdser);_
removeZeroes(grdser, ser))
checkK(1,1/2, modpol1K,fIota1) -- 1000 sec
checkK(2,29, modpol1K,fIota1) -- 1100 sec
checkK(3,29/2, modpol1K,fIota1) -- 1200 sec
checkK(4,1, modpol1K,fIota1) -- 700 sec
-- gamma:=matrix[[2, -1], [29, -14]] -- iota:=q^(-1)+104+4372*q+96256*q^2+1240002*q^3+10698752*q^4+74428120*q^5+O(q^6) -- K:=-29*q^(-29)-1682*q^(-28)-34336*q^(-27)-344636*q^(-26)-2696768*q^(-25)-16658412*q^(-24)-89603968*q^(-23)+O(q^(-22)) -- grd:=1740
-- gamma:=matrix[[1, -28], [28, -783]] -- iota:=q^(-29)+O(q^(-22)) -- K:=-q^(-29)-82*q^(-28)-2600*q^(-27)-41788*q^(-26)-412216*q^(-25)-3206652*q^(-24)-20285600*q^(-23)+O(q^(-22)) -- grd:=3364
-- gamma:=matrix[[0, -1], [1, 0]] -- iota:=q^(-29)+O(q^(-22)) -- K:=q^(-29)+82*q^(-28)+2600*q^(-27)+41788*q^(-26)+412216*q^(-25)+3206652*q^(-24)+20285600*q^(-23)+O(q^(-22)) -- grd:=3364
-- gamma:=matrix[[1, 0], [0, 1]] -- iota:=q^(-1)+104+4372*q+96256*q^2+1240002*q^3+10698752*q^4+74428120*q^5+O(q^6) -- K:=29*q^(-29)+1682*q^(-28)+34336*q^(-27)+344636*q^(-26)+2696768*q^(-25)+16658412*q^(-24)+89603968*q^(-23)+O(q^(-22)) -- grd:=1740
checkK(1,1/2, modpolNK,fIotaN) -- 1050 sec
checkK(2,29, modpolNK,fIotaN) -- 650 sec
checkK(3,29/2, modpolNK,fIotaN) -- 910 sec
checkK(4,1, modpolNK,fIotaN) -- 1000 sec
-- gamma:=matrix[[2, -1], [29, -14]] -- iota:=q^(-29)+O(q^(-22)) -- K:=-29*q^(-29)-1682*q^(-28)-34336*q^(-27)-344636*q^(-26)-2696768*q^(-25)-16658412*q^(-24)-89603968*q^(-23)+O(q^(-22)) -- grd:=3364
-- gamma:=matrix[[1, -28], [28, -783]] -- iota:=q^(-1)+104+4372*q+96256*q^2+1240002*q^3+10698752*q^4+74428120*q^5+O(q^6) -- K:=-q^(-29)-82*q^(-28)-2600*q^(-27)-41788*q^(-26)-412216*q^(-25)-3206652*q^(-24)-20285600*q^(-23)+O(q^(-22)) -- grd:=1740
-- gamma:=matrix[[0, -1], [1, 0]] -- iota:=q^(-1)+104+4372*q+96256*q^2+1240002*q^3+10698752*q^4+74428120*q^5+O(q^6) -- K:=q^(-29)+82*q^(-28)+2600*q^(-27)+41788*q^(-26)+412216*q^(-25)+3206652*q^(-24)+20285600*q^(-23)+O(q^(-22)) -- grd:=1740
-- gamma:=matrix[[1, 0], [0, 1]] -- iota:=q^(-29)+O(q^(-22)) -- K:=29*q^(-29)+1682*q^(-28)+34336*q^(-27)+344636*q^(-26)+2696768*q^(-25)+16658412*q^(-24)+89603968*q^(-23)+O(q^(-22)) -- grd:=3364
toPolAN(pol) ==> map(coerce, pol)$PolynomialFunctions2(ZZ, AN)
ev1K := eval(toPolAN modpol1K, x, jtN);
evNK := eval(toPolAN modpolNK, x, jtN);
polK := gcd(ev1K, evNK) -- K(tau) is a root of polK
kcands := radicalRoots(polK/1, y)$RadicalSolvePackage(AN);
btaucands := [simplifyRadicals(24/nn*retract(x*zN^2)@AN) for x in kcands]
Select the right candidate¶
We have now two candidates for $p_1(h(\tau_N))=B(\tau_N)$ and two candidates for $p_2(h(\tau_N))=A(\tau_N)$ and must select the right one that actually yields the correct $\frac{1}{\pi}$ formula.
[ataucands,btaucands]
In our case it is rather easy. We get the following candidate formulas for $\frac{1}{\pi}$.
\begin{gather*} \frac{2\sqrt{2}}{99^2} \sum_{n=0}^\infty \frac{(4n)!}{n!^4} \frac{26390n+1103}{396^{4n}} \qquad\qquad(1) \end{gather*}
\begin{gather*} \frac{2\sqrt{2}}{99^2} \sum_{n=0}^\infty \frac{(4n)!}{n!^4} \frac{26390n-1103}{396^{4n}} \qquad\qquad(2) \end{gather*}
\begin{gather*} \frac{2\sqrt{2}}{99^2} \sum_{n=0}^\infty \frac{(4n)!}{n!^4} \frac{-26390n+1103}{396^{4n}} \qquad\qquad(3) \end{gather*}
\begin{gather*} \frac{2\sqrt{2}}{99^2} \sum_{n=0}^\infty \frac{(4n)!}{n!^4} \frac{-26390n-1103}{396^{4n}} \qquad\qquad(4) \end{gather*}
We know already that one pair $(a,b)$ of candidates yields a $\frac{1}{\pi}$ formula and thus we must exclude the three other cases by algebraic means.
Let us state some simple fact.
Lemma (Bounds): Assume that $(y_n)_{n\ge0}$ is a sequence of positive rational numbers, $a,b,t\in\setR$, $a_l \le a \le a_u$, $b_l \le b \le b_u$, $t_l \le t \le t_u$, $a_l,a_u,b_l,b_u,t_l,t_u\in\setQ$. Assume that $\sum_{n=0}^\infty n y_n t^n$ is absolutely convergent. Let $0<m\in\setN$ and define \begin{alignat*}{2} A_m &:= \sum_{n=m+1}^\infty n y_n t^n, &\qquad B_m &:= \sum_{n=m+1}^\infty y_n t^n. \end{alignat*}
Let $A_l \le A_m \le A_u$, $B_l \le B_m \le B_u$, and $S_l < S < S_u$ with $A_l,A_u, B_l,B_u,S_l,S_u\in\setQ$. Assume that for each interval the lower and upper bound have the same sign. Depending on the signs of these bounds one can construct rational bounds $\Sigma_l$ and $\Sigma_u$ such that \begin{align*} \Sigma_l &\le \sum_{n=0}^\infty (a n + b) y_n t^n = a \left(\sum_{n=0}^m n y_n t^n + A_m\right) + b \left(\sum_{n=0}^m y_n t^n + B_m\right) \le \Sigma_u \end{align*}
If $S_u< \Sigma_l$ or $\Sigma_u <S_l$, then $S \neq \sum_{n=0}^\infty (a n + b) y_n t^n$.
The selection is done with an approximation of $\pi$ that is already known to be exact up to a certain amount of digits. The amount of pre-knowledge depends on the convergence rate of $y_n h(\tau_N)^n$ as $n$ tends to infinity.
We bring the factor $\frac{2\sqrt{2}}{99^2}$ to the other side. Thus, in our case we must check, \begin{align*} S := \frac{99^2}{2\sqrt{2}\pi} &\stackrel{?}{=} \sum_{n=0}^\infty \frac{(4n)!}{n!^4} \frac{\pm 26390n \pm 1103}{396^{4n}} \end{align*}
Clearly, the following inequation holds. \begin{align*} 24^n &= \le 256^n \prod_{k=1}^n \left(1-\frac{1}{2\,k}\right)\, \left(1-\frac{3}{4\,k}\right)\, \left(1-\frac{1}{4\,k}\right) = \frac{(4n)!}{(n!)^4}. \end{align*}
Furthermore, we have \begin{align*} ny_n &= n \frac{(4n)!}{(n!)^4} = n \cdot 256^n \prod_{k=1}^n \left(1-\frac{1}{2\,k}\right)\, \left(1-\frac{3}{4\,k}\right)\, \left(1-\frac{1}{4\,k}\right). \end{align*} We can easily check that for $n\le 10$ the inequation $n y_n < 320^ n$ holds. For $n>10$ we have $n \le \left(\frac{5}{4}\right)^n$, i.e. $n\cdot 256^n \le \left(\frac{5\cdot 256}{4}\right)^n=320^n$, and \begin{align*} n y_n &\le 320^n \prod_{k=1}^n \left(1-\frac{1}{2\,k}\right)\, \left(1-\frac{3}{4\,k}\right)\, \left(1-\frac{1}{4\,k}\right) \le 320^n. \end{align*} Therefore, for all $n\ge0$ \begin{align*} 24^n &\le y_n \le ny_n \le 320^n. \end{align*}
Assume $0 < 320 z_N < \frac{1}{2}$. We get \begin{align*} \frac{(24z_N)^{m+1}}{1 - 24 z_N} &< \sum_{n=m+1}^\infty y_n z_N^n = B_m \le A_m = \sum_{n=m+1}^\infty n y_n z_N^n < \frac{(320 z_N)^{m+1} }{1 - 320 z_N} \end{align*} and obtain rational bounds \begin{align*} A_l:=B_l &:= \left(\frac{24}{396^4}\right)^{m+1} \le B_m \le A_m \le 2 \left(\frac{320}{396^4}\right)^{m+1} =: A_u=:B_u. \end{align*}
Let us find rational bounds for $S$. Let us take a rational interval for $\pi$.
pi48 := 3141592653589793238462643383279502884197169399375/10^48; --correct digits of pi
pii := (pi48 - 5*10^(-48)) .. (pi48 + 4*10^(-48))
We can also create a rational interval for $\sqrt{2}$ that is as small as we want it.
dgts := 50;
digits(dgts+3);
truncatedSqrt2 := rationalApproximation(sqrt(2.0),dgts);
cf2 := continuedFraction truncatedSqrt2;
sq2i := convergents(cf2).51 .. convergents(cf2).52
[low(sq2i)^2 < 2, 2 < high(sq2i)^2]
Therefore, $S_l < S < S_u$ are given by
[Sl := 99^2/2/high(sq2i)/high(pii), Su := 99^2/2/low(sq2i)/low(pii)]
Since we have integer candidates for $a$ and $b$, we can choose $a_l=a=a_u$, $b_l=b=b_u$ and $t_l=t=t_u$. In our case we can select $m=1$.
Let us investigate the four cases for \begin{gather*} \sum_{n=0}^\infty \frac{(4n)!}{n!^4} \frac{a n+ b}{396^{4n}}. \end{gather*}
m:=1;
Al := Bl := (24*zN)^(m+1);
Au:=Bu:=2*(320*zN)^(m+1);
yy ==> gosperFactor
bound(a,b,A,B) ==> retract(a*(reduce(_+,[yy(n)*n*zN^n for n in 0..m])+A) _
+ b*(reduce(_+,[yy(n) *zN^n for n in 0..m])+B))@QQ
a := 26390; b := 1103;
Case 1:
Sigmal := bound(a,b,Al,Bl), Sigmau := bound(a,b,Au,Bu)
Su < Sigmal, Sigmau < Sl
Case 2:
Sigmal := bound(a,-b,Al,Bu), Sigmau := bound(a,-b,Au,Bl)
Su < Sigmal, Sigmau < Sl
Case 3:
Sigmal := bound(-a,b,Au,Bl), Sigmau := bound(-a,b,Al,Bu)
Su < Sigmal, Sigmau < Sl
Case 4:
Sigmal := bound(-a,-b,Au,Bu), Sigmau := bound(-a,-b,Al,Bl)
Su < Sigmal, Sigmau < Sl
In cases 2-4 there is at least one condition of the above Lemma fulfilled, so these candidates must be dismissed. The only remaining case is \begin{align*} \frac{99^2}{2\sqrt{2}\pi} &= \sum_{n=0}^\infty \frac{(4n)!}{n!^4} \frac{26390n + 1103}{396^{4n}} \end{align*} which therefore proves the formula for $\pi$.