
(* lambda_16 *)
<<cox.m
Print["\n\n             SALEM 16\n"]

DiStrat[salem16]

(*---------------------------------Salem factor---------------------------*)

prime=3
s = salem16
qs = QEuler[s]
fs = Mat[s]
ps=PrimeTwists[qs,fs,prime]
ps=Select[ps,DimNeeded[qs,fs,#] <= 22&]
Print["Number of primes: ",Length[ps]]

polys = ps[[1]]
{ob,certificate} = AllOb[qs,fs,polys]
Print["All twists at this prime are obstructed: ",ob]
