(*Programme accompanying the article "A Riccati differential 
  equation and free subgroup numbers for lifts of $PSL_2(\Z)$ 
  modulo powers of primes" by Christian Krattenthaler and
  Thomas M\"uller*)

(*The function FreeGF[p,e] below computes the rational function which
  provides the generating function for the number of free subgroups
  of index 4lambda in H(4) when the numbers are reduced modulo p^e*)

(*The polynomials P_n(z) and Q_n(z) in Theorem 1*)
PP[n_] := PP[n] = Factor[Sum[p2[n, k] z^(n - k), {k, 0, n}]];
QQ[n_] := QQ[n] = Factor[Sum[q2[n, k] z^(n - k), {k, 0, n}]];

pp = Pochhammer;

q2[n_, k_] := Module[{A,B,CC,DD,EE},
 A=EE=2;B=4;CC=1;DD=0;
 If[k < 0 || k > n, 
  0, (-1)^n B^(n - k)/k!/
    pp[EE/B - k, 
     2 k + 1] (Sum[(k + j)!/(k - j)!/
         j! (n - j)!/(n - k)! pp[-EE/B + j + 1, 
         k - j] pp[(A + 2 CC - EE)/2/B, j], {j, 0, 
        k}] pp[(A + 2 CC + EE)/2/B, n + 1] - 
     Sum[(k + j)!/(k - j)!/j! (n - j)!/(n - k)! pp[EE/B + j + 1, 
         k - j] pp[(A + 2 CC + EE)/2/B, j], {j, 0, 
        k}] pp[(A + 2 CC - EE)/2/B, n + 1])]
    ]

p2[n_, k_] := Module[{A,B,CC,DD,EE},
 A=EE=2;B=4;CC=1;DD=0;
 If[k < 0 || k > n, 
  0, (-1)^(n + 1) B^(n - k)/k!/pp[EE/B - k, 2 k + 1]/2/
    CC (Sum[(k + j)!/(k - j)!/j! (n - j)!/(n - k)! pp[-EE/B + j + 1, 
         k - j] pp[(A + 2 CC - EE)/2/B, 
         j] (A + If[k == 0, -EE, (2 k j B - (k - j) EE)/(k + j)]), {j,
         0, k}] pp[(A + 2 CC + EE)/2/B, n + 1] - 
     Sum[(k + j)!/(k - j)!/j! (n - j)!/(n - k)! pp[EE/B + j + 1, 
         k - j] pp[(A + 2 CC + EE)/2/B, 
         j] (A + If[k == 0, EE, (2 k j B + (k - j) EE)/(k + j)]), {j, 
        0, k}] pp[(A + 2 CC - EE)/2/B, n + 1])]
   ]

(*Function which converts the rational number r, viewed modulo m,
  into an integer*)
Ratmod[r_, m_] := 
 Mod[LinearSolve[{{Denominator[r]}}, {1}, Modulus -> m]*Numerator[r], 
   m][[1]]

(*Function which transforms a rational function Expr=P(z)/Q(z),
  with Q(z) in factored form, into a cleaner form 
  (this concerns the numerator)*)
Putzen[Expr_,m_]:=Module[{ZZ,NN,Erg,const},
   const=(Denominator[Expr]/.z->0);
   ZZ=Numerator[Expr]/const;
   ZZ=Sum[Ratmod[Coefficient[ZZ,z,ii],m]z^ii,
                          {ii, 0, Exponent[ZZ,z]}];
   NN=Denominator[Expr]/const;
   ZZ/NN
   ]

(*The product Pi_+*Pi_-*)
Prod[n_]:=Product[(4l+1)(4l+3),{l,0,n}]

(*The function FreeGF[p,e] computes the rational function which
  provides the generating function for the number of free subgroups
  of index 6lambda in PSL_2(Z) when the numbers are reduced modulo p^e.
  The result is output as partial fraction expansion of the form
  predicted by Theorem 11.*)
FreeGF[p_,e_]:=Module[{Erg,Var,P,Q,P1,Q1,P2,Q2,nn,kk,const},
   If[!PrimeQ[p],Print[p," is no prime number"];Return[]];
   If[p==2,Print["This function is not able to handle the prime p=2"];
           Return[]];
   If[p==3,
      P1=1;
      Q1=1;
      kk=0;
      While[!IntegerQ[Prod[kk*p]/p^e],
            kk++;
         ];
      P2=PolynomialMod[PP[kk*p],p^e];
      Q2=PolynomialMod[QQ[kk*p],p^e];
      P={};Q={};
      P={P1};Q={Q1};
      Erg=P1/Q1;
      For[kk=1,kk<=e,kk++,
          Var=Together[(P2/Q2-Erg)Q1^(kk+1)/p^kk];
          Var=Factor[Var,Modulus->p];
          Var=PolynomialMod[Var,p];
          Erg=Erg+Var*p^kk/Q1^(kk+1);
         ];      
     ];
   If[Mod[p,4]==1,
      P1=PolynomialMod[PP[(p-1)/4],p];
      Q1=Factor[QQ[(p-1)/4],Modulus->p];
      kk=0;
      While[!IntegerQ[Prod[(p-1)/4+kk*p]/p^e],
            kk++;
         ];
      P2=PolynomialMod[PP[(p-1)/4+kk*p],p^e];
      Q2=PolynomialMod[QQ[(p-1)/4+kk*p],p^e];
      P={};Q={};
      If[(Q1/.z->0)=!=1,
         For[kk=1,kk<=Length[Q1],kk++,
             If[NumberQ[Q1[[kk]]],
                Erg=1,
                If[Head[Q1[[kk]]]===Power,
                   Erg=Q1[[kk,1]];
                   Var=Q1[[kk,2]],
                   Erg=Q1[[kk]];
                   Var=1;
                  ];
                const=Erg/.z->0;
                Q=Join[Q,
                       {(Sum[Ratmod[Coefficient[Erg,z,ii]/const,p]z^ii,
                             {ii, 0, Exponent[Erg,z]}])^Var}];
               ];
            ];
         Q1=Product[Q[[kk]],{kk,1,Length[Q]}];
        ];
      P={P1};Q={Q1};
      Erg=P1/Q1;
      For[kk=1,kk<=e,kk++,
          Var=Together[(P2/Q2-Erg)Q1^(kk+1)/p^kk];
          Var=Factor[Var,Modulus->p];
          Var=PolynomialMod[Var,p];
          Erg=Erg+Var*p^kk/Q1^(kk+1);
         ];      
     ];
   If[Mod[p,4]==3&&p=!=3,
      P1=PolynomialMod[PP[(p-3)/4],p];
      Q1=Factor[QQ[(p-3)/4],Modulus->p];
      kk=0;
      While[!IntegerQ[Prod[(p-3)/4+kk*p]/p^e],
            kk++;
         ];
      P2=PolynomialMod[PP[(p-3)/4+kk*p],p^e];
      Q2=PolynomialMod[QQ[(p-3)/4+kk*p],p^e];
      P={};Q={};
      If[(Q1/.z->0)=!=1,
         For[kk=1,kk<=Length[Q1],kk++,
             If[NumberQ[Q1[[kk]]],
                Erg=1,
                If[Head[Q1[[kk]]]===Power,
                   Erg=Q1[[kk,1]];
                   Var=Q1[[kk,2]],
                   Erg=Q1[[kk]];
                   Var=1;
                  ];
                const=Erg/.z->0;
                Q=Join[Q,
                       {(Sum[Ratmod[Coefficient[Erg,z,ii]/const,p]z^ii,
                             {ii, 0, Exponent[Erg,z]}])^Var}];
               ];
            ];
         Q1=Product[Q[[kk]],{kk,1,Length[Q]}];
        ];
      P={P1};Q={Q1};
      Erg=P1/Q1;
      For[kk=1,kk<=e,kk++,
          Var=Together[(P2/Q2-Erg)Q1^(kk+1)/p^kk];
          Var=Factor[Var,Modulus->p];
          Var=PolynomialMod[Var,p];
          Erg=Erg+Var*p^kk/Q1^(kk+1);
         ];      
     ];
   Erg=Apart[Erg];
   If[Head[Erg]===Plus,
      Erg=Sum[Putzen[Erg[[i]], p^e], {i, 1, Length[Erg]}],
      Erg=Putzen[Erg,p^e]
     ];
   Erg
   ]







