As with any code, buyer beware.

function LocalDualGram(L,q)

Mat := MatrixAlgebra(ResidueClassRing(q), Rank(L));

U := Kernel(Mat!GramMatrix(L));

B := [ L!v : v in Basis(U) ] cat [ q*v : v in Basis(L) ];

return MinkowskiGramReduction((1/q)*GramMatrix(sub< L | B >):

Canonical := true);

end function;

function fiber(D,N,p)

Q := QuaternionAlgebra(D);

O := QuaternionOrder(Q,N);

L := LeftIdealClasses(O);

OO := QuaternionOrder(Q,N*p);

LL := LeftIdealClasses(OO);

OriginList := [ Explode([Index(L,I): I in L | IsIsomorphic(I,lideal

ALList := [0 : II in LL];

OrderGrams := [ ReducedGramMatrix(RightOrder(II)) : II in LL];

for A in SequenceToSet(OrderGrams) do

Lat := LatticeWithGram(A);

K := [ j : j in [1..#OrderGrams] | OrderGrams[j] eq A ];

B := LocalDualGram(Lat,p);

for i in K do

ALList[i] := Explode([ j : j in K | IsIsometric(LatticeWithGram(1/(Norm(I)*Norm(J)) * GramMatrix(Conjugate(J)*I)),LatticeWithGram(B)) where I := LL[i] where J := LL[j] ]);

end for;

end for;

W := ZeroMatrix(Integers(),#LL,#LL);

for i in [1..#LL] do

W[i,ALList[i]] := 1;

end for;

TerminusList := RowSequence(Transpose(W*Matrix(Integers(), #LL,1,OriginList)))[1];

ALOperators := AssociativeArray();

for l in PrimeDivisors(D*N) do

ALList := [0 : II in LL];

for A in SequenceToSet(OrderGrams) do

Lat := LatticeWithGram(A);

K := [ j : j in [1..#OrderGrams] | OrderGrams[j] eq A ];

B := LocalDualGram(Lat,l);

for i in K do

ALList[i] := Explode([ j : j in K | IsIsometric(LatticeWithGram(1/(Norm(I)*Norm(J)) * GramMatrix(Conjugate(J)*I)),LatticeWithGram(B)) where I := LL[i] where J := LL[j] ]);

end for;

end for;

ALOperators[l] := ALList;

end for;

return OriginList , TerminusList, [ExactQuotient(#Units(RightOrder(II)),2) : II in LL], ALOperators;

end function;