(* ::Package:: *)

(*This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program.  If not, see <http://www.gnu.org/licenses/>.*)


BeginPackage["RaduRK`",{"math4ti2`"}];

RISC`Version`RequiredMathematicaVersion[10.4];

(*Needs["math4ti2`"];*)

SetVar1::usage = "Defining variable q"

SetVar2::usage = "Defining variable t"

isRstar::usage = "Determines whether an eta quotient is modular over a classical modular curve of level N"

Pmr::usage = "Determines the number of factors on the left-hand side of a prospective RK identity"

pord::usage = "Gives the order of a modular function at infinity"

etaGenerators::usage = "Given a positive integer N, calling etaGenerators[N] returns the set of generators for the monoid of eta quotient modular functions over Gamma_0(N) with a pole only at infinity."

AB::usage = "Given a set of functions S, calling AB[S] returns an algebra basis for S.  In particular, if S=etaGenerators[N], calling AB[S] returns an algebra basis for the set of eta quotients over Gamma_0(N) with a pole only at infinity."

MW::usage = "Determines whether some f is a member of a C[t] module Z"

repMat::usage = "Gives a set of matrix representatives for each cusp of the classical modular curve of level N"

repMat0::usage = "Gives a set of matrix representatives in matrix form for each cusp of the classical modular curve of level N"

etaOrder::usage = "Gives the order of an eta quotient over a modular curve of level N at a cusp represented by a given matrix"

Delta::usage = "Given a positive integer M, an integer-valued vector r indexed by the divisors of M, we can define an arithmetic function a(n) whose generating function is a quotient of q-rising factorials with q-exponents indexed by the divisors of M, and powers indexed by r.  Now given an arithmetic progression mn+j in which 0 <= j < m and m>1, and a positive integer N, calling Delta[N,M,r,m,j] will return True if the generating function for a(mn+j) can be made into a modular function over Gamma_0(N), and False otherwise.";

minN::usage = "Given a generating function defined by M and r, and an arithmetic sequence mn+j, determines the minimal level of a modular curve over which a modular function can be constructed"

RK::usage = "Given a positive integer M, an integer-valued vector r indexed by the divisors of M, an arithmetic progression mn+j in which 0 <= j < m and m>1, and a positive integer N such that Delta[N,M,r,m,j] returns True, calling RKDelta[N,M,r,m,j] will return a computed Ramanujan--Kolberg identity, if one exists over the given congruence subgroup.  We assume that a certain theoretical condition on the space of eta quotients for the sake of increased speed, see Theorem 5 of companion paper.";

RKMan::usage = "Identical to RKDelta except that we also use a previously computed algebra basis A={t,{1,g1,...gv}}.  Calling RKDelta[N,M,r,m,j,A] will return a computed Ramanujan--Kolberg identity, if one exists over the given congruence subgroup.";

Begin["`Private`"]



SetVar1[x_]:=q=x;
SetVar2[y_]:=t=y;

w[r_]:=
 Sum[
  r[[d]],
  {d,1,Length[r]}
  ];


sigI[r_,N_]:=
 Sum[
  Divisors[N][[d]] r[[d]],
  {d,1,Length[Divisors[N]]}
  ];


sig0[r_,N_]:=
 Sum[
  N/Divisors[N][[d]] r[[d]],
  {d,1,Length[Divisors[N]]}
  ];


PiP[r_,N_]:=
 Product[
  Divisors[N][[d]]^(Abs[r[[d]]]),
  {d,1,Length[Divisors[N]]}
  ];


isRstar[r_,N_]:=
 If[
  AllTrue[r,IntegerQ]==True
   &&w[r]==0
   &&Mod[sigI[r,N],24]==0
   &&Mod[sig0[r,N],24]==0
   &&IntegerQ[Sqrt[PiP[r,N]]]==True,
  True,
  False
  ];


Pmr[N_,M_,r_,m_,j_]:=
 DeleteCases[
  DeleteDuplicates[
   Table[
    If[
     GCD[li,24 m]==1,
     Mod[
      (li^2) j+(((li)^2-1)/24) sigI[r,M],
      m
      ]
     ],
    {li,0,24 m-1}
    ]
   ],
  Null
  ];


twoPower[d_]:=
 If[
  PrimePowerQ[d]==True
   &&Divisible[d,2]==True,
  Log[2,d],
  0
  ];


twoPowdiv[n_]:=
 Max[
  Table[
   twoPower[Divisors[n][[i]]],
   {i,1,Length[Divisors[n]]}
   ]
  ];


S[M_,r_]:=
 twoPowdiv[
  Product[
   (Divisors[M][[i]])^(Abs[r[[i]]]),
   {i,1,Length[r]}
   ]
  ];


jd[M_,r_]:=
 If[
  S[M,r]==0,
  Product[
   (Divisors[M][[i]])^(Abs[r[[i]]]),
   {i,1,Length[r]}
   ],
  Product[
    (Divisors[M][[i]])^(Abs[r[[i]]]),
    {i,1,Length[r]}
    ]/S[M,r]
  ];


primediv[N_]:=
 DeleteCases[
  Table[
   If[
    PrimeQ[Divisors[N][[i]]]==True,
    Divisors[N][[i]],
    0
    ],
   {i,1,Length[Divisors[N]]}
   ],
  0
  ];


Delta[N_,M_,r_,m_,t_]:=
 If[
  SubsetQ[primediv[N],primediv[m]]==True
   &&SubsetQ[Divisors[m N],Divisors[M]]==True
   &&Mod[GCD[m^2-1,24] (m N^2)/M sig0[r,M],24]==0
   &&Mod[GCD[m^2-1,24] N w[r],8]==0
   &&Divisible[N,(24 m)/GCD[GCD[m^2-1,24] (-24 t-sigI[r,M]),24 m]]==True
   &&Implies[Divisible[m,2]==True,(Mod[GCD[m^2-1,24] N,4]==0&&Mod[N S[M,r],8]==0)||(Mod[S[M,r],2]==0&&Mod[N-N jd[M,r],8]==0)]==True,
  True,
  False
  ];


DeltaN[M_,r_,m_,j_]:=
 Module[
  {n0=2},
  While[
   Delta[n0,M,r,m,j]==False,
   n0=n0+1
   ];
  n0
  ];


pord[f_]:=
 If[
  -Exponent[Series[f,{q,0,1}],q,Min]==-Infinity,
  0,
  -Exponent[Series[f,{q,0,1}],q,Min]
  ];


LC[f_]:=
 Coefficient[Series[f,{q,0,pord[f]}],q,-pord[f]];
J[T_,z_,f_]:=(pord[f]-pord[z])/pord[T];
P[T_,z_,f_]:=
 If[
  Mod[pord[f],pord[T]]==Mod[pord[z],pord[T]]
   &&J[T,z,f]>=0,
  T^(J[T,z,f]) LC[f]/LC[z T^(J[T,z,f])],
  0
  ];


pordList[Z_,T_]:=
 Table[
  Mod[pord[Z[[i]]],pord[T]],
  {i,1,Length[Z]}
  ];


pordList1[Z_,T_,g_]:=
 If[
  Length[
    DeleteCases[
     Table[
      If[
       Mod[pord[Z[[i]]],pord[T]]==Mod[pord[g],pord[T]],
       pord[Z[[i]]],
       0
       ],
      {i,1,Length[Z]}
      ],
     0
     ]
    ]>=1,
  Join[
    DeleteCases[
      Table[
       If[
        Mod[pord[Z[[i]]],pord[T]]==Mod[pord[g],pord[T]],
        pord[Z[[i]]],
        0
        ],
       {i,1,Length[Z]}
       ],
      0
      ][[1]],
    {pord[g]}
    ][[1]],
  pord[g]
  ];


u1[zi_,zj_,T_]:=
  zj-LC[zj]/(LC[zi] (LC[T])^((pord[zj]-pord[zi])/pord[T])) zi T^((pord[zj]-pord[zi])/pord[T]);


u[zi_,zj_,T_]:=
 If[
  pord[u1[zi,zj,T]]!=0,
  u1[zi,zj,T],
  0
  ];


iPart0[T_,Z_,g_,i_,j_]:=
 If[
  i==j,
  t^(J[T,Z[[i]],g]) LC[g]/LC[Z[[i]] T^(J[T,Z[[i]],g])],
  0
  ];
i2Part[i_,x_]:=
 If[
  i==1,
  x,
  0
  ];


MC[T_,Z_,f_]:=
 Module[
  {g=f},
  While[
   MemberQ[pordList[Z,T],Mod[pord[g],pord[T]]]==True&&pordList1[Z,T,g]<=pord[g]&&pord[g]>0, 
   Do[
    If[
     Mod[pord[g],pord[T]]==Mod[pord[Z[[i]]],pord[T]],
     If[
      J[T,Z[[i]],g]>=0,
      g=Series[g-Z[[i]] T^(J[T,Z[[i]],g]) (LC[g]/LC[Z[[i]] T^(J[T,Z[[i]],g])]),
        {q,0,1}
        ],
      Break[]
      ]
     ],
    {i,1,Length[Z]}
    ]
   ]
  ;
  If[
   pord[g]==0,
   True,
   False
   ]
  ];


MW[T_,Z_,f_]:=
 Module[
  {g=f,M=Table[0,{i,1,Length[Z]}]},
  While[
   MemberQ[pordList[Z,T],Mod[pord[g],pord[T]]]==True&&pordList1[Z,T,g]<=pord[g]&&pord[g]>0, 
   Do[
    If[
     Mod[pord[g],pord[T]]==Mod[pord[Z[[i]]],pord[T]],
     If[
      J[T,Z[[i]],g]>=0,M=M+Table[iPart0[T,Z,g,i,j],{j,1,Length[Z]}];
      g=
       Normal[
        Series[g-Z[[i]] T^(J[T,Z[[i]],g]) LC[g]/LC[Z[[i]] T^(J[T,Z[[i]],g])],{q,0,0}]
        ],
      Break[]
      ]
     ],
    {i,1,Length[Z]}
    ]
   ];
  If[
   pord[g]==0,
   M+Table[i2Part[i,Coefficient[Series[g,{q,0,2}],q,0]],{i,1,Length[Z]}],
   {}
   ]
  ];


VB0[m_]:=
 RowReduce[
  Table[
   Coefficient[Normal[Series[m[[i]],{q,0,pord[m[[i]]]}]],q,j],
   {i,1,Length[m]},
   {j,Min[Table[-pord[m[[l]]],{l,1,Length[m]}]],-1}
   ]
  ];


VB1[m_]:=
 Transpose[
  LinearSolve[Transpose[Table[Coefficient[Normal[Series[m[[i]],{q,0,pord[m[[i]]]}]],q,j],
     {i,1,Length[m]},
     {j,Min[Table[-pord[m[[l]]],{l,1,Length[m]}]],-1}
     ]
    ]
   ,Transpose[VB0[m]]
   ]
  ];


VB2[m_]:=
 If[
  Length[m]==Dimensions[VB1[m]][[1]],
  m,
  Join[
   m,
   Table[
    0,
    {j,1,Dimensions[VB1[m]][[1]]-Length[m]}
    ]
   ]
  ];


VB3[m_]:=
 DeleteCases[
  VB1[m] . VB2[m],
  0
  ];


VB[m_]:=
 Reverse[
  VB3[m]
  ];


MB[m_]:=
 Module[{L={},M=m},
  While[TrueQ[M==L]==False,
   Set[L,VB[M]];Set[M,VB[M]];
   Do[
    If[
     Mod[pord[L[[i]]],pord[L[[1]]]]==Mod[pord[L[[j]]],pord[L[[1]]]],
     M=
      DeleteCases[
       Join[
        Delete[L,j],{u[L[[i]],L[[j]],L[[1]]]}
        ],0
       ]
     ],
    {i,1,Length[L]},
    {j,i+1,Length[L]}
    ]
   ]
  ;
  M
  ];


MBt[m_]:=
 Delete[
  MB[m],
  1
  ];


AB[m_]:=
 Module[
  {L={},M=m,T=1},
  While[
   TrueQ[L==M]==False,
   Set[L,MBt[M]];
   Set[T,MB[M][[1]]];
   Set[M,MBt[M]];
   Do[
    If[
     MC[T,Join[{1},L],L[[i]] L[[j]]]==False,
     M=Union[L,{T,L[[i]] L[[j]]}]
     ],
    {i,1,Length[L]},
    {j,1,Length[L]}
    ]
   ]
  ;
  Join[
   {T},
   {Join[{1},M]}
   ]
  ];


div[N_]:=
 Delete[
  Divisors[N],
  1
  ];


Eta[d_]:=
 q^(d/24) QPochhammer[q^d,q^d];


A[N_]:=
 Table[
  Divisors[N][[i]]/(24 GCD[Divisors[N][[i]],N/Divisors[N][[i]]]) ((GCD[Divisors[N][[j]],N/Divisors[N][[i]]]^2-Divisors[N][[j]])/Divisors[N][[j]]),
  {i,2,Length[Divisors[N]]},
  {j,2,Length[Divisors[N]]}
  ];


cap[N_]:=
 LCM@@
  LCM@@Denominator[
    A[N]
    ];


etaSys0[N_]:=
 Table[
  Sum[
    cap[N] A[N][[i]][[j]] x[j],
    {j,1,Length[A[N][[1]]]}
    ]>=0,
  {i,1,Length[A[N][[1]]]}
  ];


etaSys1[N_]:=
 Join[
  {
   x[0]+Sum[x[j],{j,1,Length[div[N]]}]==0,
   x[0]+Sum[div[N][[j]] x[j],
      {j,1,Length[div[N]]}]+24 z1==0,
   N x[0]+Sum[N/div[N][[j]] x[j],{j,1,Length[div[N]]}
      ]+24 z2==0
   },
  Table[
   Sum[
      pp[N,i,div[N][[j]]] x[j],
      {j,1,Length[div[N]]}
      ]+2 yr[i]==0,
   {i,1,Length[primeF[N]]}
   ]
  ];


etaGenerators0[N_]:=
 zsolve[
   Join[etaSys1[N],etaSys0[N]],
   Join[{x[0]},
    Table[
     x[j],
     {j,1,Length[div[N]]}
     ],
    {z1,z2},
    Table[yr[i],{i,1,Length[primeF[N]]}]
    ]
   ][[2]];


etaGenerators[N_]:=
 Module[{c=etaGenerators0[N],d={}},
  d=Sort[
    Table[
     Product[
      (Eta[Divisors[N][[j]]]/Eta[1])^(c[[i]][[j]]),{j,1,Length[Divisors[N]]}
      ],
     {i,1,Length[c]}],pord[#1]<=pord[#2]&]
  ;
  d
  ];


sqTest[d_,N_]:=
 If[
  Divisible[N,d]==True&&SquareFreeQ[d]==True,
  d,
  0
  ];


maxSqFr[N_]:=
 TakeLargest[
   Table[
    sqTest[Divisors[N][[i]],N],
    {i,1,Length[Divisors[N]]}
    ],
   1
   ][[1]];


sRep1[N_,d_]:=
 DeleteCases[
  Table[
   If[
    GCD[x,d]==1&&
     GCD[x,GCD[d,N/d]]==1,
    x,
    0
    ],
   {x,1,d}
   ],
  0
  ];


sRep0[N_,d_]:=
 Module[
  {L=sRep1[N,d],L0={sRep1[N,d][[1]]}},
  Do[
   If[
    Product[
      Mod[L[[x]]-L0[[y]],GCD[d,N/d]],
      {y,1,Length[L0]}
      ]!=0,
    L0=Union[L0,{L[[x]]}]
    ],
   {x,2,Length[L]}
   ];
  L0
  ];


sRep[N_,d_]:=
 Table[
  sRep0[N,d][[x]]/d,
  {x,1,Length[sRep0[N,d]]}
  ];


rep[N_]:=
 Module[
  {M={},d=1},
  While[
   d<Length[Divisors[N]]+1,
   M=Union[M,sRep[N,Divisors[N][[d]]]];d++
   ]
  ;M
  ];


rep0[N_]:=
 Table[
  {Numerator[rep[N][[i1]]],Denominator[rep[N][[i1]]]},
  {i1,1,Length[rep[N]]}
  ];


bd0[q_]:=
 FindInstance[
  Numerator[q] d-b Denominator[q]==1,
  {b,d},
  Integers
  ];


bd[q_]:=
 {
  bd0[q][[1]][[1]][[2]],
  bd0[q][[1]][[2]][[2]]
  };


repMat[N_]:=
 Table[
  {
   {Numerator[rep[N][[i]]],bd[rep[N][[i]]][[1]]},
   {Denominator[rep[N][[i]]],bd[rep[N][[i]]][[2]]}
   },
  {i,1,Length[rep[N]]}
  ];


repMat0[N_]:=
 Table[
  {
    {Numerator[rep[N][[i]]],bd[rep[N][[i]]][[1]]},
    {Denominator[rep[N][[i]]],bd[rep[N][[i]]][[2]]}
    }//MatrixForm,
  {i,1,Length[rep[N]]}
  ];


IdTest[y_,N_]:=
 If[
  Divisible[y[[2]][[1]],N]==True,
  y,
  0
  ];


repMatId[N_]:=
 DeleteCases[
   Table[
    IdTest[repMat[N][[i]],N],
    {i,1,Length[repMat[N]]}
    ],
   0
   ][[1]];


repMatNid[N_]:=
 DeleteCases[
  repMat[N],
  repMatId[N]
  ];


repMatNid0[N_]:=
 DeleteCases[
  repMat0[N],
  repMatId[N]//MatrixForm
  ];


p[y_,M_,r_,m_]:=
 Min[
  Table[
   1/24 Sum[
     r[[i]] (
       GCD[
          Divisors[M][[i]] (y[[1]][[1]]+GCD[m^2-1,24] l y[[2]][[1]]),m y[[2]][[1]]
          ]^2/(Divisors[M][[i]] m)
       ),
     {i,1,Length[Divisors[M]]}
     ],
   {l,0,m-1}
   ]
  ];


etaOrder[y_,N_,s_]:=
 Sum[
  (N/(24 y[[2]][[1]] GCD[y[[2]][[1]],N/y[[2]][[1]]])) s[[i]] GCD[Divisors[N][[i]],y[[2]][[1]]]^2/Divisors[N][[i]],
  {i,1,Length[Divisors[N]]}
  ];


pStar[y_,N_,s_]:=
 1/24 Sum[
   s[[i]] GCD[Divisors[N][[i]],y[[2]][[1]]]^2/Divisors[N][[i]],
   {i,1,Length[Divisors[N]]}
   ];


cuspSys0[f_,g_,N_]:=
 {
  mc g[[1]]-f[[1]]+nc f[[2]]+N jc==0,
  g[[2]]-mc f[[2]]+N kc==0
  };


cuspIdF[f_,g_,N_]:=
 Module[
  {
   x=zsolve[
     cuspSys0[f,g,N],
     {mc,nc,jc,kc}
     ],y={}
   },
  y=Join[
    {GCD[x[[1]][[1]][[1]],N]},
    Table[GCD[x[[3]][[i]][[1]],N],{i,1,Length[x[[3]]]}]
    ]
  ;
  If[Min[y]==1,Print["Same Cusps"],Print["Distinct Cusps"]]
  ];


cuspIdF0[f_,g_,N_]:=
 Module[
  {
   x=zsolve[
     cuspSys0[f,g,N],
     {mc,nc,jc,kc}
     ],y={}
   },
  If[
   Length[x[[1]]]==0,
   y={},
   y=Join[
     {GCD[x[[1]][[1]][[1]],N]},
     Table[GCD[x[[3]][[i]][[1]],N],{i,1,Length[x[[3]]]}]
     ]
   ]
  ;
  If[
   Length[y]==0,0,
   If[
    Min[y]==1,
    1,
    0
    ]
   ]
  ];


cuspCheck[f_,N_]:=
 Module[
  {x1=1,j1=1,R=rep0[N]},
  While[
   j1<=Length[R],
   If[
    cuspIdF0[f,R[[j1]],N]==1,x1=rep[N][[j1]]
    ];j1++
   ];
  x1
  ];


v[N_,M_,r_,m_,t_,PM_]:=
 Sum[
  1/m (1-m^2) (24 PM[[j]]+sigI[r,M]),
  {j,1,Length[PM]}
  ];


primeF[N_]:=
 Table[
  FactorInteger[N][[j]][[1]],
  {j,1,Length[FactorInteger[N]]}
  ];


pp[N_,i_,d_]:=
 If[
  Mod[IntegerExponent[d,primeF[N][[i]]],2]==1,
  1,
  0
  ];


prodS[N_,M_,r_,m_,t_,PM_]:=
 Product[
   (m Divisors[M][[i]])^(Abs[r[[i]]]),{i,1,Length[Divisors[M]]}]^(Length[PM]);


sSys[N_,M_,r_,m_,t_,PM_]:=
 Join[
  {
   Sum[s[j],{j,1,Length[Divisors[N]]}]==-Length[PM] w[r],
   Sum[Divisors[N][[j]] s[j],{j,1,Length[Divisors[N]]}]+24 mu11==Mod[-v[N,M,r,m,t,PM]-Length[PM] m sigI[r,M],24],
   Sum[N/Divisors[N][[j]] s[j],{j,Length[Divisors[N]]}]+24 mu22==Mod[-Length[PM] (m N)/M sig0[r,M],24]
   },
  Table[Sum[pp[N M,i,Divisors[N][[j]]] s[j],{j,1,Length[Divisors[N]]}]+2 ep[i]==pp[N M,i,prodS[N,M,r,m,t,PM]],{i,1,Length[primeF[N M]]}]
  ];


capP[N_,j_]:=
 LCM@@
  Table[
   Denominator[
    (GCD[Divisors[N][[i]],repMatNid[N][[j]][[2]][[1]]]^2/(24 Divisors[N][[i]]))],{i,1,Length[Divisors[N]]}
   ];


denP[N_,M_,r_,m_]:=
 Table[
  Denominator[
   p[repMatNid[N][[j]],M,r,m]
   ],
  {j,1,Length[repMatNid[N]]}
  ];


cA[N_,r_,M_,m_,j_]:=LCM@@{denP[N,M,r,m][[j]],capP[N,j]};


cAP[N_,j_]:=LCM@@{capP[N,j]};


MuTest[N_,M_,r_,m_,t_,PM_]:=
 Module[{x=1,R=repMatNid[N],cA=Table[cA[N,r,M,m,i],{i,1,Length[repMatNid[N]]}]},
  x=Table[
    Sum[
      (cA[[j]]) s[i] (GCD[Divisors[N][[i]],R[[j]][[2]][[1]]]^2/(24 Divisors[N][[i]])),
      {i,1,Length[Divisors[N]]}
      ]>=-(cA[[j]]) Length[PM] p[R[[j]],M,r,m],
    {j,1,Length[R]}
    ];
  x
  ];


sCand00t[N_,M_,r_,m_,t_,PM_]:=
 Module[
  {x={},A=sSys[N,M,r,m,t,PM],B=MuTest[N,M,r,m,t,PM]},
  x=
   zsolve[
    Join[A,B],
    Join[
     Table[s[j],{j,1,Length[Divisors[N]]}],
     {mu11,mu22},
     Table[ep[i],{i,1,Length[primeF[N M]]}]
     ]
    ]
  ;
  x];


sCand0t[N_,M_,r_,m_,t_,PM_]:=
 Module[
  {x={},A=sSys[N,M,r,m,t,PM],B=MuTest[N,M,r,m,t,PM]},
  x=zsolve[
     Join[A,B],
     Join[
      Table[s[j],{j,1,Length[Divisors[N]]}],
      {mu11,mu22},
      Table[ep[i],{i,1,Length[primeF[N M]]}]
      ]
     ][[1]]
  ;
  x];


scand[N_,M_,r_,m_,t_,PM_]:=
 Module[
  {Z={},T=sCand0t[N,M,r,m,t,PM]},
  Z=Table[
    Table[
     T[[j]][[i]],{i,1,Length[Divisors[N]]}],
    {j,1,Length[T]}
    ];
  Z
  ];


sMinCand0[N_,M_,r_,m_,t_,PM_]:=
 Module[
  {X=scand[N,M,r,m,t,PM],f={}},
  f=X[[1]];
  If[
   X=={},
   f={0},
   Do[
    If[
     Abs[
       Floor[
        p[{{1,0},{0,1}},M,r,m]+pStar[{{1,0},{0,1}},N,X[[j]]]
        ]
       ]<Abs[
       Floor[
        p[{{1,0},{0,1}},M,r,m]+pStar[{{1,0},{0,1}},N,f]
        ]
       ],
     f=X[[j]],
     f=f+0
     ],
    {j,1,Length[X]}
    ]
   ]
  ;
  f
  ];


sMinCand[N_,M_,r_,m_,t_,PM_]:=
If[
scand[N,M,r,m,t,PM]=={},
{},
sMinCand0[N,M,r,m,t,PM]
];


ordE[s_,N_,M_,r_,m_]:=
 Abs[
  Floor[
   p[{{1,0},{0,1}},M,r,m]+pStar[{{1,0},{0,1}},N,s]
   ]
  ];


sSysK[N_]:=
 Join[
  {
   Sum[mu[j],{j,1,Length[Divisors[N]]}]==0,
   Sum[Divisors[N][[j]] mu[j],{j,1,Length[Divisors[N]]}]+24 mu33==0,
   Sum[N/Divisors[N][[j]] mu[j],{j,Length[Divisors[N]]}]+24 mu44==0
   },
  Table[Sum[pp[N,i,Divisors[N][[j]]] mu[j],{j,1,Length[Divisors[N]]}]+2 ep2[i]==0,{i,1,Length[primeF[N]]}
   ]
  ];


MuTestK[N_]:=
 Module[{x=1,R=repMatNid[N],cA=Table[cAP[N,i],{i,1,Length[repMatNid[N]]}]},
  x=Table[
    Sum[
      (cA[[j]]) mu[i] (GCD[Divisors[N][[i]],R[[j]][[2]][[1]]]^2/(24 Divisors[N][[i]])),
      {i,1,Length[Divisors[N]]}
      ]>=1,
    {j,1,Length[R]}
    ];
  x
  ];


muCand00t[N_]:=
 Module[
  {x={},A=sSysK[N],B=MuTestK[N]},
  x=
   zsolve[
    Join[A,B],
    Join[
     Table[mu[j],{j,1,Length[Divisors[N]]}],
     {mu33,mu44},
     Table[ep2[i],{i,1,Length[primeF[N]]}]
     ]
    ]
  ;
  x];


muCand0t[N_]:=
 Module[
  {x={},A=sSysK[N],B=MuTestK[N]},
  x=zsolve[
     Join[A,B],
     Join[
      Table[mu[j],{j,1,Length[Divisors[N]]}],
      {mu33,mu44},
      Table[ep2[i],{i,1,Length[primeF[N]]}]
      ]
     ][[1]]
  ;
  x];


mucand[N_]:=
 Module[
  {Zm={},Tm=muCand0t[N]},
  Zm=Table[
    Table[
     Tm[[j]][[i]],{i,1,Length[Divisors[N]]}],
    {j,1,Length[Tm]}
    ];
  Zm
  ];


muMinCand[N_]:=
 Module[
  {Xm=mucand[N],fm={}},
  fm=Xm[[1]];
  If[
   Xm=={},
   fm={0},
   Do[
    If[
     Abs[
       sigI[Xm[[j]],N]
       ]<Abs[
       sigI[fm,N]
       ],
     fm=Xm[[j]],
     fm=fm+0
     ],
    {j,1,Length[Xm]}
    ]
   ]
  ;
  fm
  ];


muE[N_]:=
 Module[
  {
   muv=muMinCand[N]
   },1;
  Product[
   Eta[Divisors[N][[i]]]^(muv[[i]]),
   {i,1,Length[muv]}
   ]
  ];


muAB[X_,h_]:=
 Module[
  {
   Yf=
    Join[{X[[1]]},Table[X[[2]][[l]],{l,2,Length[X[[2]]]}]],
   Zf={},
   Zm={},
   Zout={},
   hm=pord[h]
   },
  Zf=Table[pord[Yf[[j]]],{j,1,Length[Yf]}];
  Zm=Table[Mod[pord[Yf[[j]]],hm],{j,1,Length[Yf]}];
  While[
   Length[Zm]<hm-1,
   Do[
    If[
     MemberQ[Mod[l1 Zf[[1]]+Zf[[l2]],hm],Zm]==False,
     Yf=Join[Yf,{Yf[[1]]^(l1) Yf[[l2]]}];Zm=Join[Zm,{Mod[l1 Zf[[1]]+Zf[[l2]],hm]}]
     ],
    {l1,1,hm},
    {l2,1,Length[Zf]}
    ]
   ]
  ;
  Zout=
   Sort[
    Yf,
    pord[#1]<=pord[#2]&
    ];
  Join[{h},{Join[{1},DeleteCases[Table[If[pord[Zout[[k]]]==hm,-1,Zout[[k]]],{k,1,Min[Length[Zout],hm]}],-1]]}]
  ];


mutK[X_,N_]:=
 Module[
  {
   Yf=
    Join[{X[[1]]},Table[X[[2]][[l]],{l,2,Length[X[[2]]]}]],
   Zf={},
   Zm={},
   Zmout={},
   h=muE[N],
   hm=0,
   km=0
   },
  hm=pord[h];
  Zf=Table[pord[Yf[[j]]],{j,1,Length[Yf]}];
  Zm=Table[Mod[pord[Yf[[j]]],hm],{j,1,Length[Yf]}];
  While[
   Length[Zm]<hm-1,
   Do[
    If[
     MemberQ[Mod[l1 Zf[[1]]+Zf[[l2]],hm],Zm]==False,
     Zf=Join[Zf,{l1 Zf[[1]]+Zf[[l2]]}];Zm=Join[Zm,{Mod[l1 Zf[[1]]+Zf[[l2]],hm]}]
     ],
    {l1,1,hm},
    {l2,1,Length[Zf]}
    ]
   ]
  ;
  Zmout=
   Sort[
    Zf,
    #1<=#2&
    ];
  km=Sum[
    If[
     Zmout[[l]]==hm,
     0,
     Floor[Zmout[[l]]/hm]
     ],
    {l,1,Min[Length[Zmout],hm]}
    ];
  h^(km)
  ];


f1[s_,N_,M_,r_,m_,PM_]:=
Product[
Eta[Divisors[N][[i]]]^(s[[i]]),
{i,1,Length[Divisors[N]]}
] Product[
q^((24 PM[[i]]+sigI[r,M])/(24 m)),
{i,1,Length[PM]}
]


fLHS0[s_,N_,M_,r_,m_,PM_]:=
 Module[
  {
   L=Series[Product[QPochhammer[q^(Divisors[M][[d]]),q^(Divisors[M][[d]])]^(r[[d]]),{d,1,Length[Divisors[M]]}],{q,0,m (pord[f1[s,N,M,r,m,PM]]+1)}
     ],
   g=
    Normal[Series[
      f1[s,N,M,r,m,PM],{q,0,0}
      ]]
   },
  1;
  Normal[
   Series[
    g Product[
      Sum[
       Coefficient[L,q,m n+PM[[i]]] q^n,
       {n,0,pord[f1[s,N,M,r,m,PM]]}
       ],
      {i,1,Length[PM]}
      ],
    {q,0,0}
    ]
   ]
  ];


fLHS[s_,N_,M_,r_,m_,PM_]:=
If[s=={},"No Membership",fLHS0[s,N,M,r,m,PM]];


RKCase1[N_,M_,r_,m_,j_]:=
Module[
{
Gx=etaGenerators[N],
X={},
PM1=Pmr[N,M,r,m,j],
f=1,
s0={},
R={},
pref=1,
return1={},
return2={}
},
Print["\!\(\*UnderscriptBox[\(\[Product]\), \(\[Delta] | M\)]\)(\!\(\*SuperscriptBox[\(q\), \(\[Delta]\)]\);\!\(\*SuperscriptBox[\(q\), \(\[Delta]\)]\)\!\(\*SubsuperscriptBox[\()\), \(\[Infinity]\), SubscriptBox[\(r\), \(\[Delta]\)]]\) = \!\(\*UnderoverscriptBox[\(\[Sum]\), \(n = 0\), \(\[Infinity]\)]\)a(n)\!\(\*SuperscriptBox[\(q\), \(n\)]\)"];
Print["\!\(\*FrameBox[
 RowBox[{
  RowBox[{
   SubscriptBox[\"f\", \"1\"], 
   RowBox[{
    RowBox[{\"(\", \"q\", \")\"}], \"\[CenterDot]\", 
    RowBox[{
     UnderscriptBox[\"\[Product]\", 
      RowBox[{
       RowBox[{\"j\", \"'\"}], \"\[Element]\", 
       RowBox[{
        SubscriptBox[\"P\", 
         RowBox[{\"m\", \",\", \"r\"}]], 
        RowBox[{\"(\", \"j\", \")\"}]}]}]], 
     RowBox[{
      UnderoverscriptBox[
       RowBox[{\"  \", \"\[Sum]\"}], 
       RowBox[{\"   \", 
        RowBox[{\"n\", \"=\", \"0\"}]}], \"\[Infinity]\"], \"a\", 
      RowBox[{\"(\", 
       RowBox[{\"mn\", \"+\", 
        RowBox[{\"j\", \"'\"}]}], \")\"}], 
      SuperscriptBox[\"q\", \"n\"]}]}]}]}], \" \", \"=\", \" \", 
  RowBox[{
   UnderscriptBox[\"\[Sum]\", 
    RowBox[{\"g\", \"\[Element]\", \"AB\"}]], 
   RowBox[{
    RowBox[{\"g\", \"\[CenterDot]\", 
     SubscriptBox[\"p\", \"g\"]}], 
    RowBox[{\"(\", \"t\", \")\"}]}]}]}]]\)"];
Print["Modular Curve: \!\(\*SubscriptBox[\(X\), \(0\)]\)(N)"];
X=AB[Gx];
s0=sMinCand[N,M,r,m,j,PM1];
f=fLHS[s0,N,M,r,m,PM1];
R=MW[X[[1]],X[[2]],f];
return1=
{
"N:",
"{M,(\!\(\*SubscriptBox[\(r\), \(\[Delta]\)]\)\!\(\*SubscriptBox[\()\), \(\[Delta] | M\)]\)}:",
"m:",
"\!\(\*SubscriptBox[\(P\), \(m, r\)]\)(j):",
"\!\(\*SubscriptBox[\(f\), \(1\)]\)(q):",
"t:",
"AB:",
"{\!\(\*SubscriptBox[\(p\), \(g\)]\)(t): g\[Element]AB}",
"Common Factor:"
};
return2=
{
N,
{M,r},
m,
PM1,
If[
s0=={},
"Select Another N",
TraditionalForm[f1[s0,N,M,r,m,PM1]]],
TraditionalForm[X[[1]]],
TraditionalForm[X[[2]]],
If[
R=={}||s0=={},
"No Membership",
R
],
If[
GCD@@Flatten[CoefficientList[R,t]]==1||R=={}||s0=={},
"None",
GCD@@Flatten[CoefficientList[R,t]]
]
};
Grid[Transpose[{return1,return2}],Dividers->{{False,True},{False,True,True,True,True,True,True,True,True}},Spacings->{1,2}]
];


RKCase1Man[N_,M_,r_,m_,j_,Ab_]:=
Module[
{
X=Ab,
PM1=Pmr[N,M,r,m,j],
f=1,
s0={},
R={},
return1={},
return2={}
},
Print["\!\(\*UnderscriptBox[\(\[Product]\), \(\[Delta] | M\)]\)(\!\(\*SuperscriptBox[\(q\), \(\[Delta]\)]\);\!\(\*SuperscriptBox[\(q\), \(\[Delta]\)]\)\!\(\*SubsuperscriptBox[\()\), \(\[Infinity]\), SubscriptBox[\(r\), \(\[Delta]\)]]\) = \!\(\*UnderoverscriptBox[\(\[Sum]\), \(n = 0\), \(\[Infinity]\)]\)a(n)\!\(\*SuperscriptBox[\(q\), \(n\)]\)"];
Print["\!\(\*FrameBox[
 RowBox[{
  RowBox[{
   SubscriptBox[\"f\", \"1\"], 
   RowBox[{
    RowBox[{\"(\", \"q\", \")\"}], \"\[CenterDot]\", 
    RowBox[{
     UnderscriptBox[\"\[Product]\", 
      RowBox[{
       RowBox[{\"j\", \"'\"}], \"\[Element]\", 
       RowBox[{
        SubscriptBox[\"P\", 
         RowBox[{\"m\", \",\", \"r\"}]], 
        RowBox[{\"(\", \"j\", \")\"}]}]}]], 
     RowBox[{
      UnderoverscriptBox[
       RowBox[{\"  \", \"\[Sum]\"}], 
       RowBox[{\"   \", 
        RowBox[{\"n\", \"=\", \"0\"}]}], \"\[Infinity]\"], \"a\", 
      RowBox[{\"(\", 
       RowBox[{\"mn\", \"+\", 
        RowBox[{\"j\", \"'\"}]}], \")\"}], 
      SuperscriptBox[\"q\", \"n\"]}]}]}]}], \" \", \"=\", \" \", 
  RowBox[{
   UnderscriptBox[\"\[Sum]\", 
    RowBox[{\"g\", \"\[Element]\", \"AB\"}]], 
   RowBox[{
    RowBox[{\"g\", \"\[CenterDot]\", 
     SubscriptBox[\"p\", \"g\"]}], 
    RowBox[{\"(\", \"t\", \")\"}]}]}]}]]\)"];
Print["Modular Curve: \!\(\*SubscriptBox[\(X\), \(0\)]\)(N)"];
s0=sMinCand[N,M,r,m,j,PM1];
f=fLHS[s0,N,M,r,m,PM1];
R=MW[X[[1]],X[[2]],f];
return1=
{
"N:",
"{M,(\!\(\*SubscriptBox[\(r\), \(\[Delta]\)]\)\!\(\*SubscriptBox[\()\), \(\[Delta] | M\)]\)}:",
"m:",
"\!\(\*SubscriptBox[\(P\), \(m, r\)]\)(j):",
"\!\(\*SubscriptBox[\(f\), \(1\)]\)(q):",
"t:",
"AB:",
"{\!\(\*SubscriptBox[\(p\), \(g\)]\)(t): g\[Element]AB}",
"Common Factor:"
};
return2=
{
N,
{M,r},
m,
PM1,
If[
s0=={},
"Select Another N",
TraditionalForm[f1[s0,N,M,r,m,PM1]]],
TraditionalForm[X[[1]]],
TraditionalForm[X[[2]]],
If[
R=={}||s0=={},
"No Membership",
R
],
If[
GCD@@Flatten[CoefficientList[R,t]]==1||R=={}||s0=={},
"None",
GCD@@Flatten[CoefficientList[R,t]]
]
};
Grid[Transpose[{return1,return2}],Dividers->{{False,True},{False,True,True,True,True,True,True,True,True}},Spacings->{1,2}]
];


minN[M_,r_,m_,t_]:=
Module[
{n=1},
While[
Delta[n,M,r,m,t]==False,
n=n+1
];
n
];


divKron[N_,d_]:=
If[
Divisible[N,d]==True,
1,
0
];


repMatNidm[n_,N_]:=
Module[
{Rn=repMatNid[n],Rn1={}},
Rn1=
Table[
If[Divisible[Rn[[k]][[2]][[1]],N]==True,0,Rn[[k]]],
{k,1,Length[Rn]}
]
;
DeleteCases[Rn1,0]
];


denPm[n_,N_,M_,r_,m_]:=
Table[
Denominator[
p[repMatNidm[n,N][[j]],M,r,m]
],
{j,1,Length[repMatNidm[n,N]]}
];


capPm[n_,N_,j_]:=
LCM@@
Table[
Denominator[
(GCD[Divisors[n][[i]],repMatNidm[n,N][[j]][[2]][[1]]]^2/(24 Divisors[n][[i]]))],{i,1,Length[Divisors[n]]}
];


cAm[n_,N_,r_,M_,m_,j_]:=LCM@@{denPm[n,N,M,r,m][[j]],capPm[n,N,j]};


MuTesta[n_,N_,M_,r_,m_,t_,PM_]:=
Module[{x=1,R=repMatNidm[n,N],cAm=Table[cAm[n,N,r,M,m,i],{i,1,Length[repMatNidm[n,N]]}]},
x=Table[
Sum[
(cAm[[j]]) s[i] (GCD[Divisors[n][[i]],R[[j]][[2]][[1]]]^2/(24 Divisors[n][[i]])),
{i,1,Length[Divisors[n]]}
]>=-(cAm[[j]]) Length[PM] p[R[[j]],M,r,m],
{j,1,Length[R]}
];
x
];


sCand0taD[N_,n_,M_,r_,m_,t_,PM_,ja_]:=
Module[
{
x={},
A=sSys[n,M,r,m,t,PM],
B=MuTesta[n,N,M,r,m,t,PM],
C={
Sum[divKron[N,Divisors[n][[j]]] N/Divisors[n][[j]] s[j],{j,Length[Divisors[n]]}]+24 mu00==ja
},
D=
DeleteCases[
Table[
Ceiling[(N/Divisors[n][[k]]-Floor[N/Divisors[n][[k]]])] s[k]==0,
{k,1,Length[Divisors[n]]}
],
True
]
},
x=zsolve[
Join[A,B,C,D],
Join[
Table[s[j],{j,1,Length[Divisors[n]]}],
{mu00,mu11,mu22},
Table[ep[i],{i,1,Length[primeF[n M]]}]
]
]
;
If[
x[[1]]=={},
{},
x[[1]][[1]]
]
];


sCand0taD[N_,n_,M_,r_,m_,t_,PM_,ja_]:=
Module[
{
x={},
A=sSys[n,M,r,m,t,PM],
B=MuTesta[n,N,M,r,m,t,PM],
C={
Sum[divKron[N,Divisors[n][[j]]] N/Divisors[n][[j]] s[j],{j,Length[Divisors[n]]}]+24 mu00==ja
},
D=
DeleteCases[
Table[
Ceiling[(N/Divisors[n][[k]]-Floor[N/Divisors[n][[k]]])] s[k]==0,
{k,1,Length[Divisors[n]]}
],
True
]
},
x=zsolve[
Join[A,B,C,D],
Join[
Table[s[j],{j,1,Length[Divisors[n]]}],
{mu00,mu11,mu22},
Table[ep[i],{i,1,Length[primeF[n M]]}]
]
]
;
If[
x[[1]]=={},
{},
x[[1]][[1]]
]
];


sCand0taD1[N_,n_,M_,r_,m_,j_,PM_,ja_,Ab_]:=
Module[
{
v=sCand0taD[N,n,M,r,m,j,PM,ja],
fv1=1,
fv=1,
Rv={}
},
fv1=If[v=={},0,f1[v,n,M,r,m,PM]];
fv=If[v=={},0,fLHS[v,n,M,r,m,PM]];
Rv=MW[Ab[[1]],Ab[[2]],fv];
If[
v=={},
{},
If[
Rv=={},
{},
{
If[
Length[v]==Length[Divisors[n]],
v,
Drop[v,{Length[Divisors[n]]+1,Length[v]}]
],
Rv . Ab[[2]]/.t->Ab[[1]],
Rv,
fv,
fv1,
Ab
}
]
]
];


fLHS20[s_,N_,n_,M_,r_,m_,PM_,km_]:=
Module[
{
L=1,
g=1,
mu=muE[n]
},
L=Series[Product[QPochhammer[q^(Divisors[M][[d]]),q^(Divisors[M][[d]])]^(r[[d]]),{d,1,Length[Divisors[M]]}],{q,0,m (pord[f1[s,N,M,r,m,PM] mu^(km)]+1)}
];
g=
Normal[
Series[
f1[s,N,M,r,m,PM] mu^(km),
{q,0,0}
]
];
If[
s=={},
1,
Normal[
Series[
g Product[
Sum[
Coefficient[L,q,m n0+PM[[i]]] q^n0,
{n0,0,pord[f1[s,N,M,r,m,PM] mu^(km)]}
],
{i,1,Length[PM]}
],
{q,0,0}
]
]
]
];


fLHS2[s_,N_,n_,M_,r_,m_,PM_,km_]:=
If[s=={}||s=={0},"No Membership",fLHS20[s,N,n,M,r,m,PM,km]];


sCand0taD2[N_,n_,M_,r_,m_,j_,PM_,ja_,Ab_,mu_,muv_]:=
If[
sCand0taD1[N,n,M,r,m,j,PM,ja,Ab]=={},
{},
Module[
{
SC=sCand0taD1[N,n,M,r,m,j,PM,ja,Ab],
ord=0,
orda=0,
LH=1,
RH=1
},
ord=
Min[
Table[
(n/GCD[repMat[n][[kn]][[2]][[1]]^2,n]) (pStar[repMat[n][[kn]],n,SC[[1]]]+Length[PM] p[repMat[n][[kn]],M,r,m]),{kn,1,Length[repMat[n]]}
]
];
orda=
If[
ord<=-1,
Ceiling[Abs[ord]],
Ceiling[ord]
];
LH=
fLHS2[SC[[1]],N,n,M,r,m,PM,orda];
RH=
mu^(orda) SC[[2]];
If[
pord[Normal[Series[LH-RH,{q,0,0}]]]==0,
SC,
{}
]
]
];


RKCase2A[N_,Sc_,M_,r_,m_,j_,PM_]:=
Module[
{
return1={},
return2={},
nd=LCM[N,minN[M,r,m,j]]
},
Print["\!\(\*UnderscriptBox[\(\[Product]\), \(\[Delta] | M\)]\)(\!\(\*SuperscriptBox[\(q\), \(\[Delta]\)]\);\!\(\*SuperscriptBox[\(q\), \(\[Delta]\)]\)\!\(\*SubsuperscriptBox[\()\), \(\[Infinity]\), SubscriptBox[\(r\), \(\[Delta]\)]]\) = \!\(\*UnderoverscriptBox[\(\[Sum]\), \(n = 0\), \(\[Infinity]\)]\)a(n)\!\(\*SuperscriptBox[\(q\), \(n\)]\)"];
Print["\!\(\*FrameBox[RowBox[{RowBox[{
   SubscriptBox[\"f\", \"1\"], RowBox[{RowBox[{\"(\", \"q\", \")\"}], \"\[CenterDot]\", RowBox[{
     UnderscriptBox[\"\[Product]\", RowBox[{RowBox[{\"j\", \"'\"}], \"\[Element]\", RowBox[{
        SubscriptBox[\"P\", RowBox[{\"m\", \",\", \"r\"}]], RowBox[{\"(\", \"j\", \")\"}]}]}]], RowBox[{
      UnderoverscriptBox[RowBox[{\"  \", \"\[Sum]\"}], RowBox[{\"   \", RowBox[{\"n\", \"=\", \"0\"}]}], \"\[Infinity]\"], \"a\", RowBox[{\"(\", RowBox[{\"mn\", \"+\", RowBox[{\"j\", \"'\"}]}], \")\"}], 
      SuperscriptBox[\"q\", \"n\"]}]}]}]}], \" \", \"=\", \" \", RowBox[{
   UnderscriptBox[\"\[Sum]\", RowBox[{\"g\", \"\[Element]\", \"AB\"}]], RowBox[{RowBox[{\"g\", \"\[CenterDot]\", 
     SubscriptBox[\"p\", \"g\"]}], RowBox[{\"(\", \"t\", \")\"}]}]}]}]]\)"];
Print["Modular Curve: \!\(\*SubscriptBox[\(X\), \(0\)]\)(n)"];
return1=
{
"n:",
"{M,(\!\(\*SubscriptBox[\(r\), \(\[Delta]\)]\)\!\(\*SubscriptBox[\()\), \(\[Delta] | M\)]\)}:",
"m:",
"\!\(\*SubscriptBox[\(P\), \(m, r\)]\)(j):",
"\!\(\*SubscriptBox[\(f\), \(1\)]\)(q):",
"t:",
"AB:",
"{\!\(\*SubscriptBox[\(p\), \(g\)]\)(t): g\[Element]AB}",
"Common Factor:"
};
return2=
{
nd,
{M,r},
m,
PM,
TraditionalForm[Sc[[5]]],
TraditionalForm[Sc[[6]][[1]]],
TraditionalForm[Sc[[6]][[2]]],
Sc[[3]],
If[
GCD@@Flatten[CoefficientList[Sc[[3]],t]]==1||Sc[[3]]=={}||Sc[[1]]=={},
"None",
GCD@@Flatten[CoefficientList[Sc[[3]],t]]
]
};
Grid[Transpose[{return1,return2}],Dividers->{{False,True},{False,True,True,True,True,True,True,True,True}},Spacings->{1,2}]
];


RKCase2[N_,M_,r_,m_,j_]:=
Module[
{
ja=0,
y={},
z={},
n=LCM[N,minN[M,r,m,j]],
PM={},
mu=1,
muv={},
G=etaGenerators[N],
X={}
},
PM=Pmr[n,M,r,m,j];
mu=muE[n];
muv=muMinCand[n];
X=AB[G];
y=sCand0taD2[N,n,M,r,m,j,PM,ja,X,mu,muv];
While[
ja<=23&&sCand0taD2[N,n,M,r,m,j,PM,ja,X,mu,muv]=={},
ja++
]
;
y=sCand0taD2[N,n,M,r,m,j,PM,ja,X,mu,muv];
z=
If[
y=={},
Print["Select Another N"],
RKCase2A[N,y,M,r,m,j,PM]
];
z
];


RKCase2Man[N_,M_,r_,m_,j_,Ab_]:=
Module[
{
ja=0,
y={},
z={},
n=LCM[N,minN[M,r,m,j]],
PM={},
mu=1,
muv={}
},
PM=Pmr[n,M,r,m,j];
mu=muE[n];
muv=muMinCand[n];
y=sCand0taD2[N,n,M,r,m,j,PM,ja,Ab,mu,muv];
While[
ja<=23&&sCand0taD2[N,n,M,r,m,j,PM,ja,Ab,mu,muv]=={},
ja++
]
;
y=sCand0taD2[N,n,M,r,m,j,PM,ja,Ab,mu,muv];
z=
If[
y=={},
Print["Select Another N"],
RKCase2A[N,y,M,r,m,j,PM]
];
z
];


RK[N_,M_,r_,m_,j_]:=
If[
Delta[N,M,r,m,j]==True,
RKCase1[N,M,r,m,j],
RKCase2[N,M,r,m,j]
];


RKMan[N_,M_,r_,m_,j_,Ab_]:=
If[
Delta[N,M,r,m,j]==True,
RKCase1Man[N,M,r,m,j,Ab],
RKCase2Man[N,M,r,m,j,Ab]
];


copyright[line_, lines___] := Module[
    {txt = line <> StringJoin[Map[("\n"<>#)&,{lines}]]},
    If[$Notebooks,
        CellPrint[Cell[txt, "Text",
            FontColor -> RGBColor[0, 0, 0],
            CellFrame -> 0.5,
            Background -> RGBColor[0.796887, 0.789075, 0.871107],
            ShowAutoSpellCheck -> False (* Needed from MMA11 on *)
        ]],
        Print[txt]
    ]];

copyright[
"RaduRK: Ramanujan--Kolberg Program Version 3.4",
    "Copyright (C) 2021, Nicolas Allen Smoot <nicolas.smoot@risc.jku.at>",
    "Research Institute for Symbolic Computation",
    "Johannes Kepler Universit\[ADoubleDot]t, Linz"];
End[]
EndPackage[]





