DiscriminantOfQuaternionAlgebra = 10;
DiscriminantOfCMPoint = 52;

DB = DiscriminantOfQuaternionAlgebra;
{tk2c, Ds2, tk3c, Ds3, c2} = 
        Switch[DB, 10, {3, 3, -2, 2, 2^(2)}, 
          6, {-6, 3, 4, 1, 6^(-6)}, 
          14, {2, 2, -4, 1, 1}];

<< NumberTheory`NumberTheoryFunctions`
<< NumberTheory`AlgebraicNumberFields`
<< SmithForm`
Off[General::spell1]
Squarefree[xxx_] := 
      Module[{divx}, 
divx = 
          FactorInteger[xxx]; 
Product[
          divx[[i, 1]]^Mod[divx[[i, 2]], 2], {i, 1, 
            Length[divx]}]
]
d = 
        Abs[If[Mod[DiscriminantOfCMPoint, 4] == 0, 
            DiscriminantOfCMPoint/4, 
            DiscriminantOfCMPoint]];
DK = 
        If[Mod[Squarefree[d], 4] == 3, Squarefree[d], 
          4  Squarefree[d]];
hk = ClassNumber[-DK];
PrimeDivisors[xxx_] := 
      If[xxx == 
          1, {}, Transpose[FactorInteger[xxx]][[1]]]
TwoAdicReduce[xxx_] := 
      Module[{minpoly, X, rts, rt1, rt2, rat, 
          irrat}, 
minpoly = 
          MinimalPolynomial[xxx, X]; 
rts = 
          X /. Solve[minpoly == 0, X]; 
{rt1, rt2} = 
          If[Length[rts] == 2, rts, 
            Append[rts, rts[[1]]]]; 
rat = 
          FullSimplify[(rt1 + rt2)/2]; 
irrat = 
          FullSimplify[(rt1 - rt2)/2]; 
ratmult = 
          If[irrat == 0, 0, 
            FullSimplify[
              Sqrt[irrat^2/
                  Squarefree[
                    FullSimplify[
                      irrat^2]]]]]; 
FullSimplify[
          rat + ratmult]
]
\(ord[ppp_, ttt_] := 
      Module[{fac, t = 0, i, tt}, 
If[
          PrimeQ[ppp], 
If[ttt == 0, 
            t = Infinity, 
tt = 
              If[ppp == 2, \ TwoAdicReduce[ttt], 
                ttt]; 
fac = 
              FactorInteger[tt]; 
For[i = 1, 
              i < Length[fac] + 1, 
If[
                fac[\([i, 1]\)] == ppp, 
                t = fac[\([i, 2]\)]]; 
\(i++\)], 
            Print["First input ", ppp\ , " not Prime"]; 
            t = 
              Null
]]; 
t
]\)
RamifiedPrimes = PrimeDivisors[DK];
Divides[xxx_, yyy_] := 
      If[Mod[yyy, xxx] == 0, True, False]
stan[xxx_, ppp_] := {xxx/ppp^ord[ppp, xxx], ord[ppp, xxx]}
CharZp[xxx_, ppp_] := 
      If[ord[ppp, xxx] ≥ 0, 1, 0]
psi[xxx_] := E^(2 Pi I * 2^stan[xxx, 2][[2]])
SubsetOf[xxx_, yyy_] := 
      If[Intersection[xxx, yyy] == xxx, True, 
        False]
ModifiedJacobi[xxx_, yyy_] := 
      Module[{twoparty, resty}, 
If[
          Mod[xxx, 4] > 1  &&  Divides[2, yyy], 
          Print[{"UNDEFINED!!", xxx, yyy}], 
resty = 
            yyy/(2^ord[2, yyy]); 
twopartJac = 
            If[Divides[2, 
                yyy], (If[Mod[xxx, 4] == 0, 0, 
                    If[Mod[xxx, 8] == 1, 1, -1]])^ord[2, yyy], 
              1]; 
JacobiSymbol[xxx, resty]*
            twopartJac]]
reduce[xxx_, ppp_] := 
      Module[{denx, numx}, 
denx = 
          Denominator[xxx]; 
numx = 
          Numerator[xxx]; 
numx*
          PowerMod[denx, -1, ppp]
]
epsi[uuu_] := 
        If[Mod[uuu, 2] == 1, Mod[(uuu - 1)/2, 2], 
          Null];
omeg[uuu_] := 
        If[Mod[uuu, 2] == 1, Mod[(uuu^2 - 1)/8, 2], 
          Null];
Hilbert[xxx_, yyy_, ppp_] := 
      Module[{alpha, beta, u, v, result}, 
alpha = 
          ord[ppp, xxx]; 
beta = 
          ord[ppp, yyy]; 
If[
          ppp == 2, 
u = 
            reduce[xxx/(ppp^alpha), 8]; 
v = 
            reduce[yyy/(ppp^beta), 8], 
u = 
            reduce[xxx/(ppp^alpha), ppp]; 
v = 
            reduce[yyy/(ppp^beta), ppp]]; 
result = 
          If[ppp == 
              2, (-1)^(epsi[u]*epsi[v] + alpha*omeg[v] + 
                  beta*omeg[u]), (-1)^(alpha*beta*epsi[ppp])*
              ModifiedJacobi[u, ppp]^beta*
              ModifiedJacobi[v, ppp]^
                alpha]; 
result
]
rhop[ppp_, ttt_] := 
      If[ord[ppp, ttt] ≥ 0, 
If[
          ModifiedJacobi[-d, ppp] == 0, 
          1, 
Sum[(ModifiedJacobi[-d, ppp])^k, {k, 
              0, ord[ppp, ttt]}]], 0]
rho[ttt_] := 
      If[(ttt ∈  Integers), 
Module[{divt, 
            i}, 
divt = 
            PrimeDivisors[ttt]; 
Product[
            rhop[divt[[i]], ttt], {i, 1, 
              Length[divt]}]
], 
0]
WO[ddd_, lis_] := 
      Module[{res = ddd, i}, 
For[i = 1, 
          i ≤ Length[lis], 
res = 
            res/lis[[i]]^
                ord[lis[[i]], 
                  ddd]; 
i++]; 
res]
InL[mu_, ppp_] := 
      ord[ppp, mu[[1]]] ≥ 0  &&  
        ord[ppp, mu[[2]]] ≥ 0
chi[xxx_, ppp_] := 
      If[ord[ppp, xxx] == 0, Hilbert[xxx, ppp, ppp], 
        0]
\(Present[li_] := 
      Module[{i, top = "", bot = ""}, 
For[i = 1, 
          i ≤ Length[li], 
If[
            li[\([i, 2]\)] > 0, 
top = 
              top . \ ToString[li[\([i]\)]], 
bot = 
              bot . \ ToString[
                  li[\([i]\)]]]; 
\(i++\);]; 
Print[
          top]; 
Print["------------------------"]; 
Print[bot];
]\)
ExpandLog[xx_] := 
        Module[{fac = FactorInteger[xx]}, 
Sum[
            fac[[i, 2]]*Log[fac[[i, 1]]], {i, 1, 
              Length[fac]}]];
genus[Delt_, Dee_] := 
        2^(Length[PrimeDivisors[GCD[Delt, Dee]]] - 
              If[SubsetOf[PrimeDivisors[Delt], PrimeDivisors[Dee]], 1, 
                0]);
SizeOfZU = hk/genus[DK, DB];

q = 2;
For[i = 1, 
        i < 1000, 
Module[{pis, 
              cond = 1}, 
pis = 
              PrimeDivisors[DB]; 
If[
              Not[Mod[q, 8] == 5] || Mod[DB, q] == 0, 
              q = NextPrime[q], 
              For[j = 1, j ≤ Length[pis], 
                If[pis[[j]] == 2, 
                  If[Not[Mod[q, 8] == 5], cond = -1], 
                  If[JacobiSymbol[q, pis[[j]]] ≠ -1, 
                    cond = -1]]; 
j++]; 
If[cond == -1, q = NextPrime[q], 
                i = 10001];
];]
i++];
If[JacobiSymbol[DB, q] == 
        1, 
m = Position[Table[Mod[i^2, q], {i, 1, q}], 
            Mod[DB, q]][[1, 1]]; 
l = 
        Mod[PowerMod[2*DB, -1, q]*m, q]; 
mp = 
        2*DB*l;]

iota[{{aaa_, bbb_}, {ccc_, ddd_}}] := {{ddd, -bbb}, {-ccc, 
          aaa}}
inn[xxx_, yyy_] := FullSimplify[Tr[xxx . iota[yyy]]]
Q[xxx_] := FullSimplify[inn[xxx, xxx]/2]

o = {{1, 0}, {0, 1}};
b = {{Sqrt[q], 0}, {0, -Sqrt[q]}};
c = {{0, DB}, {1, 0}};
bc = b . c;
bcBasis[{{www_, xxx_}, {yyy_, zzz_}}] := 
      FullSimplify[{(www + zzz)/
            2, (www - zzz)/(2 Sqrt[q]), (xxx + DB yyy)/(2  
                DB), (xxx - DB yyy)/(2 DB Sqrt[q])}]
e1 = (o + b)/2;
e2 = (mp b + bc)/q;
e3 = e1 . e2; 
    L1 = FullSimplify[
          bcBasis[A o  +  (-2 A - D mp) e1 + C e2 +  
                D e3] /. {A -> -1, C -> 0, D -> 0}] . {o, 
          b, c, bc};
L2 = 
        FullSimplify[
            bcBasis[A o  +  (-2 A - D mp) e1 + C e2 +  
                  D e3] /. {A -> 0, C -> 1, D -> 0}] . {o, b, 
            c, bc};
L3 = 
        FullSimplify[
            bcBasis[A o  +  (-2 A - D mp) e1 + C e2 +  
                  D e3] /. {A -> 0, C -> 0, D -> 1}] . {o, b, 
            c, bc};
L3 = (q - 1) L2/2 + L3;
LBasis[xxx_] := 
      FullSimplify[
        bcBasis[xxx] . {{1, 0, 0, 0}, {0, 1, 0, 0}, {0, mp, -q, 
              2}, {0, -mp, q, 0}}]
eBasis[xxx_] := LBasis[xxx];

li = 80;
For[i = 0, 
        i < li, 
IT = 
          FullSimplify[
            Reduce[(Det[A L1 + B L2 + C L3] /. {C -> i}) == 
                d, {A, B}, Integers,  
              Backsubstitution -> True]]; 
If[
          IT =!= False, Ctemp = i; i = li + 999]; 
If[
          i ≤ 0, i = -i + 1, 
          If[i < li, 
            i = -i]];
];
Z = 
        A L1 + B L2 + Ctemp L3  /. 
          ToRules[Reduce[LogicalExpand[
                  IT /. {C[1] -> 0}][[1]], {A, B}, 
              Integers]];
{z0, z1, z2, z3} = LBasis[Z];

k =. ;
i =. ;
j =. ;
k = 
        k /. Solve[inn[Z, i*L1 + j*L2 + k*L3] == 0, 
              k][[1]];
U2 = 
        Coefficient[i*XX + j*YY + k*ZZ, i] /. {XX -> L1, YY -> L2, 
            ZZ -> L3};
U1 = 
        Coefficient[i*XX + j*YY + k*ZZ, j] /. {XX -> L1, YY -> L2, 
            ZZ -> L3};
U1 = 
        U1*DB q (2 z2 + (-1 + 
                    q) z3);
U2 = U2 - inn[U2, U1]/inn[U1, U1]*U1;
U2 = 
        FullSimplify[
          U2*(-4 mp^2 (-1 + q) q^2 z1^2 - 
                  8 mp (DB + mp^2 (-1 + q)) q z1 z2 + 
                  4 (DB - mp^2) (DB + 
                        mp^2 (-1 + q)) z2^2 + 
                  4 DB (DB + mp^2 (-1 + q)) q z2 z3 + 
                  DB (DB + mp^2 (-1 + q)) (-1 + 
                        q) q z3^2)/(2 q)];
k =. ;
i =. ;
j =. ;

alp1  =  
        DB q (2 z2 + (-1 + q) z3)/(2  
                q);
alp2 = -4 q (q z1 + mp z2)/(2  
                q);
alp3 = (-4 DB z2 + 4 mp (q z1 + mp z2) - 
              2 DB q z3)/(2  q);
g123 = GCD[alp1, alp2, alp3];
g12 = GCD[alp1, alp2]/g123;
g13 = GCD[alp1, alp3]/g123;
g23 = GCD[alp2, alp3]/g123;
If[g12 == 0, g12 = 1];
If[g13 == 0, g13 = 1];
If[g23 == 0, g23 = 1];
alp11 = alp1/(g12 g13 g123);
alp21 = alp2/(g12 g23 g123);
alp31 =  alp3/(g13 g23 g123);
{rr1, 
          ss1} = ExtendedGCD[alp11, 
            alp31][[2]];
{rr3, 
          ss3} = ExtendedGCD[
            alp11 g12,  -2 alp11 g12  -  
              alp31 g23][[2]];
Lm1 = alp11 g12 L2 + alp31 g23 L3;
Lm2  = 
        g13 L1  +  ss1 alp21 g12 L2 - 
          alp21 g23 rr1 L3;
If[
        GCD[Denominator[inn[Lm1, Lm2]/inn[Lm2, Lm2]], 
            Denominator[inn[Lm1, Lm2]/inn[Lm1, Lm1]]] > 2, 
        Lm1 = Lm1 + Lm2];

LmBasis[xxx_] := 
      Module[{A, B, Zero, i, j, k, lis}, 
Zero = 
          eBasis[A*Lm1 + B*Lm2 - xxx]; 
i = 
          If[eBasis[Lm1][[2]] == 0, 
            If[eBasis[Lm1][[3]] == 0, 4, 3], 
            2]; 
lis = 
          Complement[{2, 3, 4}, {i}]; 
j = 
          If[eBasis[Lm2][[lis[[1]]]] == 0, lis[[2]], 
            lis[[1]]]; 
k = Complement[{2, 3, 4}, {i, 
                j}][[1]]; 
A = 
          A /. Solve[Zero[[i]] == 0, 
                A][[1]]; 
B = 
          B /. Solve[Zero[[j]] == 0, 
                B][[1]]; 
If[
          FullSimplify[Zero[[k]]] == 0, {A, 
            B}, Print["Nope", {A, B}];]]

CreateOrtho[ppp_] := 
      Module[{num, den, Ax, Bx, IT, i, 
          j}, 
OL1 =. ; 
OL2 =. ; 
TOL =. ; 
den = 
          Denominator[inn[Lm1, Lm2]/inn[Lm1, Lm1]]; 
num = 
          Numerator[inn[Lm1, Lm2]/inn[Lm1, Lm1]]; 
If[
          MemberQ[PrimeDivisors[den], ppp], 
den = 
            Denominator[
              inn[Lm1, Lm2]/inn[Lm2, Lm2]]; 
num = 
            Numerator[inn[Lm1, Lm2]/inn[Lm2, Lm2]]; 
If[
            MemberQ[PrimeDivisors[den], 
              ppp], 
Ax =. ; 
Bx =. ; 
OL1 = Lm1; 
OL2 = 
              Ax*Lm1 + Bx*Lm2; 
i =. ; 
            j =. ; 
IT = 
              Q[i*OL1 + j*OL2]; 
Ax = 
              Ax /. Solve[
                    Coefficient[IT, i j, 1] == Coefficient[IT, i, 2], 
                    Ax][[1]]; 
Bx = 
              Bx /. Solve[
                    Coefficient[IT, j, 2] == Coefficient[IT, i, 2], 
                    Bx][[1]]; 
TOL  =  
              FullSimplify[{{1, -Ax/Bx}, {0, 
                    1/Bx}}]; 
diag = {Coefficient[IT, i, 
                  2], Coefficient[IT, i, 2]};, 
OL1 = 
              Lm2; 
OL2 = 
              den*Lm1 - 
                num*Lm2; 
TOL = {{num/den, 1}, {1/den, 
                  0}}; 
diag = 
              Module[{i, j}, {Coefficient[Q[i*OL1 + j*OL2], i, 2], 
                  Coefficient[Q[i*OL1 + j*OL2], j, 
                    2]}]];, 
OL1 = 
            Lm1; 
OL2 = 
            den*Lm2 - 
              num*Lm1; 
TOL = {{1, num/den}, {0, 
                1/den}}; 
diag = 
            Module[{i, j}, {Coefficient[Q[i*OL1 + j*OL2], i, 2], 
                Coefficient[Q[i*OL1 + j*OL2], j, 
                  2]}];
];
]
OLBasis[
          xxx_] := Transpose[
            TOL . Transpose[{LmBasis[
                    xxx]}]][[1]];
epiel[iii_, ppp_] := stan[diag[[iii]], ppp];
CreateOrtho[3];
tempdiag = diag;
If[Not[
        PrimeDivisors[
            Abs[Det[TOL]]] == {2}], 
pri = PrimeDivisors[Det[TOL]][[-1]]; 
CreateOrtho[
        pri];]
LatPrimes = 
        PrimeDivisors[
          Abs[tempdiag[[1]]*tempdiag[[2]]*diag[[1]]*
              diag[[2]]]];
diag =.

i =. 
K1 =. ;
K2 =. ;
Y =. ;
L =. ;
Te = K1*U1 + K2*U2;
Ite = Collect[inn[Te, i*Lm1 + j*Lm2], {i, j}];
If[
        Coefficient[Ite, i K1] ≠ 0, 
K1 = 
          K1 /. Solve[Coefficient[Ite, i] == L, 
                K1][[1]]; 
K2 = 
          K2 /. Solve[Coefficient[Ite, j] == Y, 
                K2][[1]], 
K1 = 
          K1 /. Solve[Coefficient[Ite, j] == L, 
                K1][[1]]; 
K2 = 
          K2 /. Solve[Coefficient[Ite, i] == Y, 
                K2][[1]]];
{K1, K2};
TeB = eBasis[Te];
TeBY = TeB /. {Y -> 1, L -> 0};
TeBL = TeB /. {L -> 1, Y -> 0};
mvBasis = Developer`HermiteNormalForm[{{TeBY[[2]], 
                TeBY[[3]], TeBY[[4]]}, {TeBL[[2]], TeBL[[3]], 
                TeBL[[4]]}}][[2]];
Lmv1 = (mvBasis[[1]] . {{L1}, {L2}, {L3}})[[1]];
Lmv2 = (mvBasis[[2]] . {{L1}, {L2}, {L3}})[[1]];
Zv = Z/(2  d);

RelMat = {LBasis[Z][[{2, 3, 4}]], LBasis[
              Lm1][[{2, 3, 4}]], LBasis[
              Lm2][[{2, 3, 4}]]};
CosetRep[RM_] := 
        Module[{SS, UU, VV}, 
{SS, {UU, VV}} = 
            IntegerSmithForm[RM]; 
VV = 
            Inverse[VV]; 
If[
            Abs[SS[[3, 3]]] == Abs[Det[RM]], Cycl = 1; 
            VV[[3, 1]]*L1 + VV[[3, 2]]*L2 + VV[[3, 3]]*L3, 
            Cycl = 0; {{VV[[3, 1]]*L1 + VV[[3, 2]]*L2 + 
                  VV[[3, 3]]*L3, 
                SS[[3, 3]]}, {VV[[2, 1]]*L1 + VV[[2, 2]]*L2 + 
                  VV[[2, 3]]*L3, 
                SS[[2, 2]]}, {VV[[1, 1]]*L1 + VV[[1, 2]]*L2 + 
                  VV[[1, 3]]*L3, 
                SS[[1, 
                    1]]}}]
];
CosRep = CosetRep[RelMat];
i =. 
inn[i*Zv, Z] == i;
LamLim = Abs[Det[RelMat]];
lam[nnn_] := 
        If[Length[CosRep] == 2, nnn*CosRep, 
          Floor[nnn/(CosRep[[1, 2]]*CosRep[[2, 2]])]*L1 + 
            Floor[Mod[
                    nnn, (CosRep[[1, 2]]*
                        CosRep[[2, 2]])]/(CosRep[[1, 2]])]*L2 + 
            Mod[nnn, CosRep[[1, 2]]]*L3];
plus[xxx_] := inn[xxx, Z]/inn[Z, Z]*Z
minus[xxx_] := xxx - plus[xxx]
lamplus[nnn_] := plus[lam[nnn]]
lamminus[nnn_] := minus[lam[nnn]]
FullSimplify[Det[lamplus[n]] ≥ 0, 
        n ∈ Intgers];
FullSimplify[Det[lamminus[n]] ≤ 0, 
        n ∈ Integers];

ModL[xxx_] := 
      Module[{bas}, 
bas = 
          Mod[eBasis[xxx], 1]; 
bas[[1]]*o + 
          bas[[2]]*c + bas[[3]]*bc + bas[[4]]*e]
EtaPlus[eta_] := plus[ModL[eta]]
EtaMinus[eta_] := minus[ModL[eta]]
kappalist[eta_, mmm_, lim_] := 
      Module[{QvalP, QvalM, top, X, last, etam, etap, lamp, lamm, n, i, 
          templist = {}}, 
etap = 
          EtaPlus[eta]; 
etam = 
          EtaMinus[eta]; 
For[n = 0, 
          n < LamLim, 
lamp = 
            lamplus[n]; 
lamm = 
            lamminus[n]; 
top = 
            Round[X /. Solve[D[Q[lamp + etap + X*Z], X] == 0, 
                    X][[1]]]; 
For[i = 0, 
            i ≤ lim, 
QvalP = 
              mmm - Q[lamp + 
                    etap + (top + (1 + i))*Z]; 
If[
              QvalP ≥ 
                0 , templist = 
                  Append[templist, {QvalP, (lamm + etam), 
                      n, (top + (1 + 
                              i))}];
]; 
QvalM = 
              mmm - Q[lamp + etap + (top - i)*Z]; 
If[
              QvalM ≥ 
                0 , 
templist = 
                  Append[templist, {QvalM, (lamm + etam), 
                      n, (top - 
                          i)}];
]; 
If[QvalM < 0  && QvalP < 0,  
              i = lim]; 
i++]; 
n++]; 
templist]

Hmu[mu_, ppp_] := 
      Module[{res = {}}, 
If[
          CharZp[mu[[1]], ppp] == 1, 
          res = Append[res, 1]]; 
If[
          CharZp[mu[[2]], ppp] == 1, 
          res = Append[res, 2]]; 
res]
Kmu[mu_, ppp_] := 
      If[CharZp[mu[[1]], ppp] == 1  &&  
          CharZp[mu[[2]], ppp] == 1, Infinity, 
        Module[{nothmu, i}, 
nothmu = 
            Complement[{1, 2}, Hmu[mu, ppp]]; 
Min[
            Table[epiel[nothmu[[i]], ppp][[2]] + 
                ord[ppp, mu[[nothmu[[i]]]]], {i, 1, 
                Length[nothmu]}]]]]
Lmu[mu_, kkk_, ppp_] := 
      Module[{it, h, i, res = {}}, 
h = 
          Hmu[mu, ppp]; 
For[i = 1, 
          i ≤ 
            Length[h], 
it = epiel[h[[i]], 
                  ppp][[2]] - kkk; 
If[
            it < 0  &&  Mod[it, 2] == 1, 
            res =  Append[res, 
                h[[i]]]]; 
i++]; 
res]
lmu[mu_, kkk_, ppp_] := Length[Lmu[mu, kkk, ppp]]
dmu[mu_, kkk_, ppp_] := 
      kkk + (1/2)*
          Sum[Min[{epiel[Hmu[mu, ppp][[i]], ppp][[2]] - kkk, 
                0}], {i, 1, Length[Hmu[mu, ppp]]}]
epimu[mu_, kkk_, ppp_] := 
      Hilbert[-1, ppp, ppp]^IntegerPart[lmu[mu, kkk, ppp]/2]*
        Module[{li, i}, 
li = 
            Lmu[mu, kkk, ppp]; 
Product[
            Hilbert[epiel[li[[i]], ppp][[1]], ppp, ppp], {i, 1, 
              Length[li]}]]
tmu[mu_, mmm_, ppp_] := 
      Module[{nothmu, res = 0, ep, i}, 
nothmu = 
          Complement[{1, 2}, Hmu[mu, ppp]]; 
For[i = 1, 
          i ≤ Length[nothmu], 
res = 
            res + diag[[nothmu[[i]]]]*
                mu[[nothmu[[i]]]]^2; 
i++]; 
mmm - res]
aye[mu_, mmm_, ppp_] := ord[ppp, tmu[mu, mmm, ppp]]
f1[mu_, mmm_, ppp_] := 
      If[Mod[lmu[mu, ord[ppp, mmm] + 1, ppp], 2] == 0, -1/ppp, 
        Hilbert[stan[mmm, ppp][[1]], ppp, ppp]/Sqrt[ppp]]
det2S[ppp_] := 
      ppp^(-ord[ppp, 4*diag[[1]]*diag[[2]]]/2)
Lp[ppp_, 
        XXX_] := (1 - 
            If[Divides[ppp, DK], 0, Hilbert[ppp, -DK, ppp]]*
              XXX/ppp)^(-1)
Wmp[mmm_, ppp_, mu_, XXX_] := 
      Module[{}, 
CreateOrtho[
          ppp]; 
nmu = OLBasis[mu]; 
If[
          ppp == 2, Wm2[mmm, nmu, XXX], 
          det2S[ppp]*Lp[ppp, XXX]*
            If[aye[nmu, mmm, ppp] ≥  Kmu[nmu, ppp], 
              Expand[(1 + (1 - 1/ppp)*
                      Sum[If[Mod[lmu[nmu, k, ppp], 2] == 0, 
                          epimu[nmu, k, ppp]*ppp^dmu[nmu, k, ppp]*XXX^k, 
                          0], {k, 1, Kmu[nmu, ppp]}])], 
              If[aye[nmu, mmm, ppp] ≥ 0, 
                Expand[(1 + (1 - 1/ppp)*
                        Sum[If[Mod[lmu[nmu, k, ppp], 2] == 0, 
                            epimu[nmu, k, ppp]*ppp^dmu[nmu, k, ppp]*XXX^k, 
                            0], {k, 1, aye[nmu, mmm, ppp]}] + 
                      epimu[nmu, aye[nmu, mmm, ppp] + 1, ppp]*
                        f1[nmu, tmu[nmu, mmm, ppp], ppp]*
                        ppp^dmu[nmu, aye[nmu, mmm, ppp] + 1, ppp]*
                        XXX^(aye[nmu, mmm, ppp] + 1))], 
                1]]]
]
Wmp0[mmm_, ppp_, mu_] := Wmp[mmm, ppp, mu, 1]
DWmp0[mmm_, ppp_, mu_] := 
      Module[{sss}, 
If[Wmp0[mmm, ppp, mu] == 0, 
          D[Wmp[mmm, ppp, mu, ppp^(-sss)], 
              sss] /. {sss -> 0}]]
ord2[mu_] := Min[{ord[2, mu[[2]]], ord[2, mu[[1]]]}]
Nmu2[mu_] := If[InL[mu, 2], {1}, {}]
dmu2[mu_, kkk_] := 
      kkk + Sum[
          Min[ord[2, diag[[1]]] - kkk, 0], {i, 1, 
            Length[Nmu[mu]]}]
pmu2[mu_, kkk_] := (-1)^
        Sum[ord[2, diag[[1]]] - kkk, {i, 1, Length[Nmu[mu]]}]
tmu2[mu_, mmm_] := If[InL[mu, 2], mmm, 0]
numu2[mu_, mmm_, kkk_] := tmu2[mu, mmm]*2^(3 - kkk)
aye2[mu_, mmm_] := ord[2, tmu2[mu, mmm]]
Kmu2[mu_, mmm_] := 
      If[InL[mu, 2], aye2[mu, mmm] + 4, 
        ord[2, diag[[1]]] + ord2[mu]]
dmu2d[mu_, kkk_] := 
      kkk + (1/2)*
          Sum[Min[{epiel[Hmu[mu, 2][[i]], 2][[2]] - kkk + 1, 
                0}], {i, 1, Length[Hmu[mu, 2]]}]
epimu2[mu_, kkk_] := 
      Module[{li, i}, 
li = 
          Lmu[mu, kkk - 1, 
            2]; 
Product[epiel[li[[i]], 
              2][[1]], {i, 1, Length[li]}]]
deltamu2[mu_, kkk_] := 
      Module[{tem}, 
tem = 
          Product[epiel[Hmu[mu, 2][[i]], 2][[2]] - kkk + 
              1, {i, 1, Length[Hmu[mu, 2]]}]; 
If[
          tem == 0, 0, 1]]
tmu2d[mu_, mmm_] := 
      Module[{nothmu, res = 0, ep, i}, 
nothmu = 
          Complement[{1, 2}, Hmu[mu, 2]]; 
For[i = 1, 
          i ≤ Length[nothmu], 
res = 
            res + diag[[nothmu[[i]]]]*
                mu[[nothmu[[i]]]]^2; 
i++]; 
mmm - res]
numu2d[mu_, mmm_, kkk_] := 
      tmu2d[mu, mmm]*2^(3 - kkk) - 
        Sum[If[epiel[Hmu[mu, 2][[i]], 2][[2]] < 
              kkk - 1, epiel[Hmu[mu, 2][[i]], 2][[1]], 
            0], {i, 1, Length[Hmu[mu, 2]]}]
aye2d[mu_, mmm_] := ord[2, tmu2d[mu, mmm]]
Kmu2d[mu_, mmm_] := 
      If[InL[mu, 2], aye2d[mu, mmm] + 30, 
        Module[{i, lis = {}}, 
For[i = 1, 
            i ≤ 2, 
If[
              ord[2, mu[[i]]] < -1, 
              lis = Append[
                  lis, epiel[i, 2][[2]] + ord[2, mu[[i]]] + 
                    1]]; 
If[
              ord[2, mu[[i]]] == -1, 
lis = 
                Append[lis, epiel[i, 2][[2]] + 
                    1]]; 
i++]; 
Min[
            lis]]]
Wm2[mmm_, mu_, XXX_] := 
      Module[{iji, jij}, 
If[
          Coefficient[Q[iji OL1 + jij OL2], iji jij, 1] == 
            0,   
Lp[2, XXX]*
            det2S[2]*(1 + 
                Sum[If[Mod[lmu[mu, k - 1, 2], 2] == 0, 
                    deltamu2[mu, k]*2^(dmu2d[mu, k] - 1)*
                      chi[epimu2[mu, k], 2]*psi[numu2d[mu, mmm, k]/8]*
                      CharZp[numu2d[mu, mmm, k]/4, 2]*XXX^k, 
                    deltamu2[mu, k]*2^(dmu2d[mu, k] - 3/2)*
                      chi[epimu2[mu, k]*numu2d[mu, mmm, k], 2]*XXX^k], {k, 1, 
                    Kmu2d[mu, 
                      mmm]}]), 
   
          Lp[2, XXX]*2^(-ord[2, diag[[1]]^2]/2)*(1 + 
                Sum[pmu2[mu, k]*2^(dmu2[mu, k] - 1)*
                    psi[numu2[mu, mmm, k]/8]*CharZp[numu2[mu, mmm, k]/4, 2]*
                    XXX^k, {k, 1, 
                    Kmu2[mu, mmm]}])
]]
Wm20[mmm_, mu_] := 
      Module[{XX}, Wm2[mmm, mu, XX] /. {XX -> 1}]
DWm20[mmm_, mu_] := 
      Module[{sss}, 
If[Wm20[mmm, mu] == 0, 
          D[Wm2[mmm, mu, 2^(-sss)], 
              sss] /. {sss -> 0}]]
ListKappa[kap_] := 
      For[iii = 1, iii ≤ Length[kap], 
        Print[{kap[[iii, 1]], eBasis[kap[[iii, 2]]], 
            OLBasis[kap[[iii, 2]]], {Wm20[kap[[iii, 1]], 
                OLBasis[kap[[iii, 2]]]], 
              DWm20[kap[[iii, 1]], OLBasis[kap[[iii, 2]]]]}, 
            Table[{Wmp0[kap[[iii, 1]], SpecialPrimes[[i]], 
                  OLBasis[kap[[iii, 2]]]], 
                DWmp0[kap[[iii, 1]], SpecialPrimes[[i]], 
                  OLBasis[kap[[iii, 2]]]]}, {i, 1, 
                Length[SpecialPrimes]}]}]; iii++]
NeededPrimes[mmm_] := 
      Union[Union[PrimeDivisors[Abs[mmm]], LatPrimes], {2}]
W[mmm_, mu_, sss_] := 
      Product[Wmp[mmm, NeededPrimes[mmm][[i]], 
          mu, NeededPrimes[mmm][[i]]^(-sss)], {i, 1, 
          Length[NeededPrimes[mmm]]}]
SumKappa[kap_] := 
      Expand[Sum[
          Module[{ss}, 
            If[kap[[i, 1]] == 0, 
              INFINITY, -2*DK^((ss + 1)/2)/hk*
                  D[W[kap[[i, 1]], kap[[i, 2]], ss], 
                    ss] /. {ss -> 0}]], {i, 1, 
            Length[kap]}]]

Print["d = ", d];
Print["Class Number = ", hk];
Print["Z = ", eBasis[Z]];
Print[""];
Print["Basis of Negative Plane:"];
Print["U1 = ", eBasis[U1]];
Print["U2 = ", eBasis[U2]];
X =. ; Y =. ;
Print["Quadratic Form: ",  Q[X U1  +  Y U2]]; 
    Print[""];
Print["Basis of Lattice over integers:"]
Print["Lm1 = ", eBasis[Lm1]];
Print["Lm2 = ", eBasis[Lm2]];
X =. ; Y =. ;
Print["Quadratic Form over this basis: ",  
        Q[X Lm1  +  Y Lm2]];
Print[""];
CreateOrtho[3];
\(If[Det[TOL] ≠ 1, 
CreateOrtho[
        2]; 
Print["Basis over the 2-adics:"]; 
Print["Lm1_2 = ", \(LmBasis[
            OL1]\)[\([1]\)], "*Lm1 + ", \(LmBasis[
            OL1]\)[\([2]\)], "*Lm2"]; 
Print["Lm2_2 = ", \(LmBasis[OL2]\)[\([1]\)], "*Lm1 + ", \(LmBasis[
            OL2]\)[\([2]\)], "*Lm2"]; 
Print["Quadratic Form: ", Q[X\ OL1 + Y\ OL2]]; 
      Print[""]; 
CreateOrtho[
        3]; 
fc = 
        PrimeDivisors[Det[TOL]]; 
If[
        fc == {2}, 
Print["Basis over p-adics for all other p:"]; 
Print["Lm1_p = ", \(LmBasis[
              OL1]\)[\([1]\)], "*Lm1 + ", \(LmBasis[
              OL1]\)[\([2]\)], "*Lm2"]; 
Print["Lm2_p = ", \(LmBasis[OL2]\)[\([1]\)], "*Lm1 + ", \(LmBasis[
              OL2]\)[\([2]\)], "*Lm2"]; 
Print["Quadratic Form: ", \ Q[X\ OL1 + Y\ OL2]]; 
        Print[""];, 
Print["Basis over p-adics for p NOT in ", 
          fc, ":"]; 
Print["Lm1_p = ", \(LmBasis[
              OL1]\)[\([1]\)], "*Lm1 + ", \(LmBasis[
              OL1]\)[\([2]\)], "*Lm2"]; 
Print["Lm2_p = ", \(LmBasis[OL2]\)[\([1]\)], "*Lm1 + ", \(LmBasis[
              OL2]\)[\([2]\)], "*Lm2"]; 
Print["Quadratic Form: ", \ Q[X\ OL1 + Y\ OL2]]; 
        Print[""]; 
CreateOrtho[
          fc[\([\(-1\)]\)]]; 
Print["Basis over p-adics for p in ", \ 
          Complement[
            fc, {2}], ":"]; 
Print["Lm1_p = ", \(LmBasis[OL1]\)[\([1]\)], "*Lm1 + ", \(LmBasis[
              OL1]\)[\([2]\)], "*Lm2"]; 
Print["Lm2_p = ", \(LmBasis[OL2]\)[\([1]\)], "*Lm1 + ", \(LmBasis[
              OL2]\)[\([2]\)], "*Lm2"]; 
Print["Quadratic Form: ", \ Q[X\ OL1 + Y\ OL2]]; 
        Print[""];]]\)
\(If[Length[CosRep] == 2, 
      Print["L/(L- + L+) generated by lam = ", eBasis[lam[1]]]; 

      Print["and is cyclic of order ", LamLim]; 

      Print["lam- = ", \(LmBasis[
            lamminus[1]]\)[\([1]\)], "*Lm1 + ", \(LmBasis[
            lamminus[1]]\)[\([2]\)], "*Lm2"]; 

      Print["lam+ = ", 
        inn[lamplus[1], Z]/inn[Z, Z], "*Z"];]\)
k2 = kappalist[0*o, Ds2, 100];
k3 = kappalist[0*o, Ds3, 100];
kp = k2;
Print["TABLE 1:"];
Print[
        Prepend[Table[{kp[[i, 
                    3]], (inn[lamplus[kp[[i, 3]]], Z]/inn[Z, Z] + 
                      kp[[i, 4]]) "Z", 
                kp[[i, 1]], -2*DK^((ss + 1)/2)/hk*
                    D[W[kp[[i, 1]], kp[[i, 2]], ss], 
                      ss] /. {ss -> 0}}, {i, 1, 
                Length[kp]}], {"i", "x", 
              ToString[Ds2] <> "-Q[x]", "k_u(m)"}] // 
          TableForm];
kp = k3;
Print["TABLE 2:"];
Print[
        Prepend[Table[{kp[[i, 
                    3]], (inn[lamplus[kp[[i, 3]]], Z]/inn[Z, Z] + 
                      kp[[i, 4]]) "Z", 
                kp[[i, 1]], -2*DK^((ss + 1)/2)/hk*
                    D[W[kp[[i, 1]], kp[[i, 2]], ss], 
                      ss] /. {ss -> 0}}, {i, 1, 
                Length[kp]}], {"i", "x", 
              ToString[Ds3] <> "-Q[x]", "k_u(m)"}] // 
          TableForm];
X =. ;
\(\(Leg = {"p", "u1", "u2", "W^*_{m,p}(s,u)", "W^*_{m,p}(0,u)", \ "d/ds W^*_{m,p}(0,u)"};\)\)
kp = k2;
Print["TABLE 3:"];
\(\(For[i1 = 1, 
        i1 ≤ Length[kp], 
pr = 
          NeededPrimes[
            kp[\([i1, 
                1]\)]]; 
Print["For u = ", \(LmBasis[
              kp[\([i1, 2]\)]]\)[\([1]\)], "*Lm1+", \(LmBasis[
              kp[\([i1, 
                  2]\)]]\)[\([2]\)], "*Lm2 = u1*Lm1_p+u2*Lm2_p  and m = ", kp[\([i1, 1]\)]]; 
Tb = 
          Table[{pr[\([j1]\)], 
              CreateOrtho[
                pr[\([j1]\)]]; \(OLBasis[
                  kp[\([i1, 2]\)]]\)[\([1]\)], \(OLBasis[
                  kp[\([i1, 2]\)]]\)[\([2]\)], 
              Simplify[
                Wmp[kp[\([i1, 1]\)], pr[\([j1]\)], kp[\([i1, 2]\)], X]], 
              Simplify[
                Wmp[kp[\([i1, 1]\)], pr[\([j1]\)], kp[\([i1, 2]\)], 
                    X] /. {X -> 1}], 
              If[Simplify[
                    Wmp[kp[\([i1, 1]\)], pr[\([j1]\)], kp[\([i1, 2]\)], 
                        X] /. {X -> 1}] == 0, 
                Simplify[
                  D[Wmp[kp[\([i1, 1]\)], pr[\([j1]\)], kp[\([i1, 2]\)], 
                        pr[\([j1]\)]^\((\(-X\))\)], 
                      X] /. {X -> 0}], ""]}, {j1, 1, 
              Length[pr]}]; 
Print[
          Prepend[Tb, Leg] // 
            TableForm]; 
Print["Thus k_u(", 
          kp[\([i1, 1]\)], ") = ", \ \(-2\)*DK^\((\((ss + 1)\)/2)\)/hk*
              D[W[kp[\([i1, 1]\)], kp[\([i1, 2]\)], ss], ss] /. {ss -> 
                0}]; 
Print[""]; 
\(i1++\)];\)\)
tk2 = SumKappa[k2];
Print["Summing these gives: kappa_0(", Ds2, ") = ", 
      tk2]
Print[""]; Print[""];
kp = k3;
Print["TABLE 4:"];
\(\(For[i1 = 1, 
        i1 ≤ Length[kp], 
pr = 
          NeededPrimes[
            kp[\([i1, 
                1]\)]]; 
Print["For u = ", \(LmBasis[
              kp[\([i1, 2]\)]]\)[\([1]\)], "*Lm1+", \(LmBasis[
              kp[\([i1, 
                  2]\)]]\)[\([2]\)], "*Lm2 = u1*Lm1_p+u2*Lm2_p  and m = ", kp[\([i1, 1]\)]]; 
Tb = 
          Table[{pr[\([j1]\)], 
              CreateOrtho[
                pr[\([j1]\)]]; \(OLBasis[
                  kp[\([i1, 2]\)]]\)[\([1]\)], \(OLBasis[
                  kp[\([i1, 2]\)]]\)[\([2]\)], 
              Simplify[
                Wmp[kp[\([i1, 1]\)], pr[\([j1]\)], kp[\([i1, 2]\)], X]], 
              Simplify[
                Wmp[kp[\([i1, 1]\)], pr[\([j1]\)], kp[\([i1, 2]\)], 
                    X] /. {X -> 1}], 
              If[Simplify[
                    Wmp[kp[\([i1, 1]\)], pr[\([j1]\)], kp[\([i1, 2]\)], 
                        X] /. {X -> 1}] == 0, 
                Simplify[
                  D[Wmp[kp[\([i1, 1]\)], pr[\([j1]\)], kp[\([i1, 2]\)], 
                        pr[\([j1]\)]^\((\(-X\))\)], 
                      X] /. {X -> 0}], ""]}, {j1, 1, 
              Length[pr]}]; 
Print[
          Prepend[Tb, Leg] // 
            TableForm]; 
Print["Thus k_u(", 
          kp[\([i1, 
              1]\)], ") = ", \ \ \((\(-2\)*DK^\((\((ss + 1)\)/2)\)/hk*
                D[W[kp[\([i1, 1]\)], kp[\([i1, 2]\)], ss], 
                  ss])\) /. {ss -> 
                0}]; 
Print[""]; 
\(i1++\)];\)\)
tk3 = SumKappa[k3];
Print["Summing these gives: kappa_0(", Ds3, ") = ", 
      tk3]
Print[""]; Print[""];
Print["To Recap:"];
Print["kappa_0(", Ds2, ") = ", tk2];
Print["kappa_0(", Ds3, ") = ", tk3];
Print["Ideal Class Number: ", hk];
Print["Field Discriminant: ", DK];
Print["Normalizing Constant: C2 = ", c2];
c1 = -SizeOfZU/2^(Length[PrimeDivisors[DB]]);
Print["Final Result:"];
FinRes = 
        Expand[c1*(tk3c*tk3 + tk2c*tk2) - 
            SizeOfZU*ExpandLog[c2]];
Print[c1, "*(", tk2c, "*kappa_0(", Ds2, ")+", 
      tk3c, "*kappa_0(", Ds3, "))-", SizeOfZU, "*Log[C2]= ", 
      FinRes]; Print[];
OutDen = "";
OutNum = "";
For[i = 1, 
        i ≤ Length[FinRes], 
If[
          FinRes[[i, 1]] < 0, 
OutDen = 
            OutDen <> ToString[FinRes[[i, 2, 1]]] <> "^{" <> 
              ToString[-FinRes[[i, 
                      1]]] <> "}", 
If[
            Length[FinRes[[i]]] == 1, 
            OutNum = OutNum <> ToString[FinRes[[i, 1]]] <> "^{1}", 
            OutNum = 
              OutNum <> ToString[FinRes[[i, 2, 1]]] <> "^{" <> 
                ToString[
                  FinRes[[i, 
                      1]]] <> "}"]]; 
i++];
If[OutNum == "", OutNum = "1"];
If[OutDen == "", OutDen = "1"];
Print[OutNum];
\(\(FracBar := 
        Module[{i, outbar = ""}, 
For[i = 1, 
            i < Max[StringLength[OutNum], StringLength[OutDen]] + 
                2, 
outbar = 
              outbar <> "-"; 
\(i++\)]; 
outbar];\)\)
Print[FracBar];
Print[OutDen];
Print["-------------------------------------------------------------\
------------------------------------"];