The Magma code used for the lower bound in Example 7.1:

%% This script computes the lower bound on a_Y^n for fixed n.

p:=3;
d:=100;
Genus:=(p-1)/2*(d-1);
prank:=0;
print "Genus", Genus, "p-rank", prank, "\n";

%% the a_Y^n should approach g_Y, since the Cartier-Manin matrix is nilpotent.

print "n Lower bound \n";
n:=1;
%% we compute the lower bound first.
maxlowerbound:=0;
for j in [0.. p-1] do
lowerboundsumj:=0;
for i in [j.. p-1] do
termij:=Floor(i*d/p) - Floor( i*d/p - (1-p^(-n))*j*d/p);
lowerboundsumj:=lowerboundsumj+ termij;
end for;
if lowerboundsumj gt maxlowerbound then
maxlowerbound:=lowerboundsumj;
end if; end for;

print n, maxlowerbound;

---

The Magma code used for the upper bound in Example 7.2:

%% This script computes the upper bound for a_Y^n. Annoyingly i cannot be
%% put in a for loop, so you have to manually paste the code p-1 times.

p:=3;
d:=100;
Genus:=(p-1)/2*(d-1);
prank:=0;
print "Genus", Genus, "p-rank", prank, "\n";
print "n", "upper bound";

n:=2;
upperbound:=0;

i:=1;
%% we compute sigma_p(d,i,n).
sigmapdin:=0;
for l in [2.. Ceiling(i*d/p)] do
if l mod p^n ne 1 then
QljinT:=1;
newj:=p-1-i;
quantity:=l-1;
for k:=0 to n-1 do
oldj:=newj;
newjbutmightexceedp:=Integers()!(quantity/p^k)*Modinv(d,p)+oldj;
newj:=newjbutmightexceedp-p*Floor(newjbutmightexceedp/p);
if newj gt oldj then
QljinT:=0;
break k;
else
quantity:=quantity+(oldj-newj)*d*p^k;
end if;
end for;
if QljinT eq 1 then
sigmapdin:=sigmapdin+1;
end if;
end if;
end for;
upperbound:=upperbound+Floor(i*d/p)-Floor(i*d/p^(n+1)) - sigmapdin;

i:=2;
%% we compute sigma_p(d,i,n).
sigmapdin:=0;
for l in [2.. Ceiling(i*d/p)] do
if l mod p^n ne 1 then
QljinT:=1;
newj:=p-1-i;
quantity:=l-1;
for k:=0 to n-1 do
oldj:=newj;
newjbutmightexceedp:=Integers()!(quantity/p^k)*Modinv(d,p)+oldj;
newj:=newjbutmightexceedp-p*Floor(newjbutmightexceedp/p);
if newj gt oldj then
QljinT:=0;
break k;
else
quantity:=quantity+(oldj-newj)*d*p^k;
end if;
end for;
if QljinT eq 1 then
sigmapdin:=sigmapdin+1;
end if;
end if;
end for;
upperbound:=upperbound+Floor(i*d/p)-Floor(i*d/p^(n+1)) - sigmapdin;

print n, upperbound;