//First compute the set of primes of bad reduction. We take advantage of the factorization of the polynomial to make the calculation easier

P<x>:=PolynomialRing(Integers());

h:=3*(2156368971259291788768*x^2 + 17631006366692164075730*x + 29699355872439635532575)*

(11059757352967285969512842208*x^2 - 7091311407633771522260131534*x + 890759660375643943081974383)*

(72471395366770846855160511270238311434544819644269483224510638908380448975178286675*x^4 + 220412834417638726170369528388500737059808808950881801569083924119940345487077117130*x^3 + 349028896642164971206777347398916786554047452255452799240009025234448247634394895917*x^2 - 226266297687093053128869990785216720827777361022207859381240801920732109106262460538*x + 92436361129062375107961184178560590459501276670160242498848213720499182945004486643);

h0:=Factorization(Factorization(h)[1][1]); // constant term

h1:=Factorization(h)[2][1];

h2:=Factorization(h)[3][1];

h3:=Factorization(h)[4][1];

F1:=Factorization(Discriminant(h1));

F2:=Factorization(Discriminant(h2));

F3:=Factorization(Discriminant(h3));

F4:=Factorization(Resultant(h1,h2));

F5:=Factorization(Resultant(h1,h3));

F6:=Factorization(Resultant(h2,h3));

FF:=F1 cat F2 cat F3 cat F4 cat F5 cat F6 cat h0;

BadPrimes:={FF[i][1]:i in [1..#FF]};

BadPrimes:=Sort([Integers()!p:p in BadPrimes]);

//Compute the primes of bad reduction for which the polynomial does admit a p-adic root

BadPrimes2:=[];

for p in BadPrimes do

R := pAdicField(p : Precision := 550);

Q<t>:=PolynomialRing(R);

f:=3*(2156368971259291788768*t^2 + 17631006366692164075730*t + 29699355872439635532575)*

(11059757352967285969512842208*t^2 - 7091311407633771522260131534*t + 890759660375643943081974383)*

(72471395366770846855160511270238311434544819644269483224510638908380448975178286675*t^4 + 220412834417638726170369528388500737059808808950881801569083924119940345487077117130*t^3 + 349028896642164971206777347398916786554047452255452799240009025234448247634394895917*t^2 - 226266297687093053128869990785216720827777361022207859381240801920732109106262460538*t + 92436361129062375107961184178560590459501276670160242498848213720499182945004486643);

if #Roots(f) eq 0 then

Append(~BadPrimes2, p);

end if;

end for;

//The remaining primes of bad reduction: [ 7, 11, 17, 8221, 63596671, 5817238919 ]

BadPrimes2;

//For each remaining prime p of bad reduction, verify that h(2) or h(4) is a quadratic nonresidue mod p

a2:=Evaluate(h,2);

a4:=Evaluate(h,4);

// for each p, at least one of the Legendre symbols is -1

for p in BadPrimes2 do

p, LegendreSymbol(a2,p), LegendreSymbol(a4,p);

end for;