Skip to content

Commit

Permalink
Merge pull request #1 from assaferan/newforms
Browse files Browse the repository at this point in the history 
Newforms
  • Loading branch information
@assaferan
assaferan committed Jun 23, 2023
2 parents 387d2e6 + dd6cb6c commit 6e284ae
Show file tree
Hide file tree
Showing 29 changed files with 1,081 additions and 1,987 deletions.
1 change: 1 addition & 0 deletions .gitattributes
View file
Original file line number Diff line number Diff line change
@@ -1,3 +1,4 @@
# Auto detect text files and perform LF normalization
* text=auto
data/* filter=lfs diff=lfs merge=lfs -text
up_to_600000.tar.gz filter=lfs diff=lfs merge=lfs -text
2 changes: 1 addition & 1 deletion .github/workflows/test.yml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@ name: Tests

on:
push:
branches: [ 'composite_level', 'main' ]
branches: [ 'newforms', 'composite_level', 'main' ]

pull_request_target:

Expand Down
2 changes: 2 additions & 0 deletions .gitignore
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,3 +21,5 @@ m4
Makefile
configure
ltmain.sh
*.sig
Tests/*.sig
285 changes: 36 additions & 249 deletions Tests/MagmaFormula.m
View file
Original file line number Diff line number Diff line change
@@ -1,263 +1,50 @@
function Sfast(N, u, t, n)
fac := Factorization(N*u);
num_sols := 1;
for f in fac do
p,e := Explode(f);
if p eq 2 then
e2 := Valuation(1-t+n, 2);
if (e2 eq 0) or (e gt e2) then
return 0;
end if;
if IsEven(t) then
num_sols *:= 2^(e-1);
end if;
else
y := t^2-4*n;
e_y := Valuation(y, p);
if e_y lt e then
if IsOdd(e_y) then
return 0;
end if;
y_0 := y div p^e_y;
// compare timings
// is_sq := KroneckerSymbol(y_0,p);
is_sq := IsSquare(Integers(p)!y_0);
if not is_sq then
return 0;
end if;
num_sols *:= p^(e_y div 2) * 2;
else
num_sols *:= p^(e div 2);
end if;
end if;
end for;
return num_sols div u;
end function;

function S(N, u, t, n)
assert N mod u eq 0;
assert (t^2 - 4 * n) mod u^2 eq 0;
return [x : x in [0..N-1] | GCD(x,N) eq 1 and (x^2 - t*x + n) mod N*u eq 0];
end function;

function phi1(N)
primes := [f[1] : f in Factorization(N)];
if IsEmpty(primes) then return N; end if;
return Integers()!(N * &*[1 + 1/p : p in primes]);
end function;

function B(N, u, t, n)
// assert Sfast(N,u,t,n) eq #S(N,u,t,n);
//return #S(N,u,t,n) * phi1(N) div phi1(N div u);
return Sfast(N,u,t,n) * phi1(N) div phi1(N div u);
end function;

function C(N, M)
a, b, c, d := Explode(M);
G_M := GCD(c, d-a, b);
return B(N, GCD(G_M, N), Trace(M), Determinant(M));
end function;

function C(N, u, t, n)
return &+[B(N, u div d, t, n) * MoebiusMu(d) : d in Divisors(u)];
end function;

function H(n)
if n lt 0 then
is_sq, u := IsSquare(-n);
return (is_sq select -u/2 else 0);
end if;
if n eq 0 then
return -1/12;
end if;
if n mod 4 in [1,2] then
return 0;
end if;

//ret := &+[ClassNumber(BinaryQuadraticForms(-n div d)) : d in Divisors(n) | IsSquare(d) and (n div d) mod 4 in [0,3] ];
ret := &+[ClassNumber(-n div d) : d in Divisors(n)
| IsSquare(d) and (n div d) mod 4 in [0,3] ];
if IsSquare(n) and IsEven(n) then
ret -:= 1/2;
end if;
if n mod 3 eq 0 and IsSquare(n div 3) then
ret -:= 2/3;
end if;
return ret;
end function;

function Phi(N, a, d)
ret := 0;
for r in Divisors(N) do
s := N div r;
// scalar := r eq s select 1/2 else 1;
g := GCD(r,s);
if GCD(N, a-d) mod g eq 0 then
alpha := CRT([a,d],[r,s]);
if (alpha ne 0) or (N eq 1) then
ret +:= EulerPhi(g);
end if;
end if;
end for;
return ret;
end function;

// Gegenbauer polynomials
function P(k, t, m)
R<x> := PowerSeriesRing(Rationals());
return Coefficient((1 - t*x+m*x^2)^(-1), k-2);
end function;

// At the moment, this seems to work when N is prime
function TraceFormulaGamma0(n, N, k)
S1 := 0;
max_abst := Floor(SquareRoot(4*n));
for t in [-max_abst..max_abst] do
for u in Divisors(N) do
if ((4*n-t^2) mod u^2 eq 0) then
S1 +:= P(k,t,n)*H((4*n-t^2) div u^2)*C(N,u,t,n);
end if;
end for;
end for;
S2 := 0;
for d in Divisors(n) do
a := n div d;
S2 +:= Minimum(a,d)^(k-1)*Phi(N,a,d);
end for;
ret := -S1 / 2 - S2 / 2;
if k eq 2 then
ret +:= &+[n div d : d in Divisors(n) | GCD(d,N) eq 1];
end if;
return ret;
end function;
import "magma_functions.m" : TraceFormulaGamma0AL, TraceFormulaGamma0ALNew;

function PhiAL(N, a, d)
return EulerPhi(N) / N;
end function;

// Now this seems to work - tested for N <= 100, even k, 2<=k<=12 and 1 <= n <= 10, n = N
// This formula follows Popa - On the Trace Formula for Hecke Operators on Congruence Subgroups, II
// Theorem 4.
// (Also appears in Skoruppa-Zagier, but this way of stating the formula was easier to work with).
function TraceFormulaGamma0AL(n, N, k)
S1 := 0;
// max_abst := Floor(SquareRoot(4*N*n));
max_abst := Floor(SquareRoot(4*N*n)) div N;
// ts := [t : t in [-max_abst..max_abst] | t mod N eq 0];
// for t in ts do
for tN in [-max_abst..max_abst] do
t := tN*N;
for u in Divisors(N) do
if ((4*n*N-t^2) mod u^2 eq 0) then
S1 +:= P(k,t,N*n)*H((4*N*n-t^2) div u^2)*C(1,1,t,N*n)
*MoebiusMu(u) / N^(k div 2-1);
end if;
end for;
end for;
S2 := 0;
for d in Divisors(n*N) do
a := n*N div d;
if (a+d) mod N eq 0 then
S2 +:= Minimum(a,d)^(k-1)*PhiAL(N,a,d) / N^(k div 2-1);
procedure testPariVSMagma(N,k : New := false, OnlyPrimes := false, Prec := 1000)
if (New and (not OnlyPrimes)) then
printf("Error! Only supports newsubspace for prime Hecke operators!\n");
assert false;
end if;
end for;
ret := -S1 / 2 - S2 / 2;
if k eq 2 then
ret +:= &+[n div d : d in Divisors(n) | GCD(d,N) eq 1];
end if;
return ret;
end function;

function A1(n,N,k)
if not IsSquare(n) then
return 0;
end if;
return n^(k div 2 - 1)*phi1(N)*(k-1)/12;
end function;

function mu(N,t,f,n)
N_f := GCD(N,f);
primes := [x[1] : x in Factorization(N) | (N div N_f) mod x[1] ne 0];
s := #[x : x in [0..N-1] | (x^2 - t*x+n) mod (GCD(f*N, N^2)) eq 0];
prod := IsEmpty(primes) select 1 else &*[ 1 + 1/p : p in primes];
return N_f * prod * s;
end function;

// This is not working yet. Probably due to class number issues.
function A2(n,N,k)
R<x> := PolynomialRing(Rationals());
max_abst := Floor(SquareRoot(4*n));
if IsSquare(n) then max_abst -:= 1; end if;
ret := 0;
for t in [-max_abst..max_abst] do
F<rho> := NumberField(x^2 - t*x+n);
rho_bar := t - rho;
p := Rationals()!((rho^(k-1) - rho_bar^(k-1)) / (rho - rho_bar));
assert p eq P(k,t,n);
t_sum := 0;
for d in Divisors(4*n - t^2) do
is_sq, f := IsSquare(d);
if is_sq then
D := (t^2-4*n) div f^2;
h := (D mod 4 in [0,1]) select ClassNumber(D) else 0;
w := #UnitGroup(Integers(QuadraticField(D)));
t_sum +:= (h/w) * mu(N,t,f,n);
end if;
end for;
ret -:= p*t_sum;
end for;
return ret;
end function;

function A3(n,N,k)
g := 1;
ds := [d : d in Divisors(N) | d^2 lt n];
ret := 0;
for d in ds do
cs := [c : c in Divisors(N) |
GCD(N div g, n div d - d) mod GCD(c, N div c) eq 0];
ret -:= d^(k-1) * &+[EulerPhi(GCD(c, N div c)) : c in cs];
end for;
is_sq, d := IsSquare(n);
if is_sq then
cs := [c : c in Divisors(N) |
GCD(N div g, n div d - d) mod GCD(c, N div c) eq 0];
ret -:= 1/2 * d^(k-1) * &+[EulerPhi(GCD(c, N div c)) : c in cs];
end if;
return ret;
end function;

function A4(n,N,k)
if k eq 2 then
return &+[t : t in Divisors(n) | GCD(N, n div t) eq 1];
end if;
return 0;
end function;

function TraceCohen(n,N,k)
return A1(n,N,k) + A2(n,N,k) + A3(n,N,k) + A4(n,N,k);
end function;

procedure testPariVSMagma(N)
cmd := Sprintf("./src/traceALbatch_sta %o %o", N, N+1);
cmd := Sprintf("./src/traceALbatch_sta %o %o %o %o %o %o", N, N+1, k, Prec,
OnlyPrimes select 1 else 0, New select 1 else 0);
System(cmd);
fname := Sprintf("./data/traces_%o.m", N);
fname := Sprintf("./data/traces_%o_%o.m", k, N);
r := Read(fname);
return_cmd := Sprintf("return traces_%o, tracesAL_%o;", N, N);
traces_full, traces_AL := eval(r cat return_cmd);
traces_magma := [TraceFormulaGamma0AL(n, N, 2) : n in [1..1000]];
assert traces_magma eq traces_AL[2..1001];
al_start := Index(r, "tracesAL");
r := r[al_start..#r];
start := Index(r, "[") + 1;
fin := Index(r, "]") - 1;
r := r[start..fin];
traces_AL := [StringToInteger(x) : x in Split(r, ",")];
trace_func := New select TraceFormulaGamma0ALNew else TraceFormulaGamma0AL;
val_list := OnlyPrimes select PrimesUpTo(Prec) else [0..Prec];
traces_magma := [trace_func(n, N, k) : n in val_list];
assert traces_magma eq traces_AL[1..#val_list];
return;
end procedure;

num_tests := 10;
max_Nk2 := 10^5;
max_level := 10000;
Ns := PrimesUpTo(max_level);
printf "testing pari vs magma implementation... N = ";
Ns := [1..max_level];
printf "testing pari vs magma implementation... (N;k) = ";
for i in [1..num_tests] do
N := Random(Ns);
max_weight := Floor(Sqrt(max_Nk2/N));
ks := [2..max_weight by 2];
k := Random(ks);
printf "(%o;%o),", N, k;
testPariVSMagma(N,k);
end for;

printf "testing pari vs magma implementation on new subspaces... (N;k) = ";
for i in [1..num_tests] do
N := Random(Ns);
printf "%o,", N;
testPariVSMagma(N);
max_weight := Floor(Sqrt(max_Nk2/N));
ks := [2..max_weight by 2];
k := Random(ks);
printf "(%o;%o),", N, k;
testPariVSMagma(N,k : New, OnlyPrimes);
end for;
printf "\n";

Loading

0 comments on commit 6e284ae

Please sign in to comment.