I have just spent a month applying polredabs to about 7 million polynomials of degree up to 21, with 100% success even though some took a long time (and I could be patient).
Each of these polynomials (monic and integral) defines an extension of Q which I know to be unramified at primes > 500000 (these are the fields over which the torsion of an elliptic curve in my database acquires extra torsion, and were supplied by Enrique Gonzalez and Filip Najman), so I was able to help the process by dividing out from the polynomial discriminant all small primes, seeing that the quotient was a square, and taking its square root before factoring it (to give the factors to polredabs).
But now although I have all the polredabs-ed polynomials (in about 92.5% of the cases the original polynomial was returned, in 7.5% of cases a new polynomial was returned), we want to add some or all of the fields to the LMFDB which will require doing many more computations with these number fields, starting with finding the maximal order. It seems a pity that polredabs() has laready computed the maximal order and the field discriminant, but that has been discarded -- as far as I can tell there is no option to polredabs() to return anything other than the new polynomial and an _expression_ giving its root in terms of the original (or vice versa). And the integral basis for the field -- but not the factored discriminant.
Please can we have an option to polredabs to return more? At least: the output polynomial's factored discriminant.
I do see from the documentation that it can return the integral basis, from which one can recover the field discriminant, and in may case this is easy to factor. But I did not think of that until after running it 7 million times...
John
PS degree 16 was by far the worst
--