Special Hilbert Class Field computation

[Kurt Foster <drsardonicus@earthlink.net>]

I've got a bnf structure, field = bnfinit(poly) with field,clgp,no =  
2.  I want the Hilbert Class Field.  A classical solution (due to  
Hilbert himself) is to find an integral element which is not a square,  
but is a square mod 4, i.e. mod 4zk.  The bnf structure gives two  
lists of suspects: one is units mod squares of units.  The other is a  
list of principal generators of P^2, where P is the ideal given by  
field.clgp.gen -- such a list can be produced using bnfisintnorm().

The thing that's got me hung up is checking for squareness mod 4.  I  
can produce a list of squares mod 4 (the square roots are defined mod  
2) and use nfbasistoalg() and nfalgtobasis() to get everything in  
terms of t_COL's of integer coefficients of basis elements to compare,  
but I'm probably being very clumsy about it.  Besides, I'd probably be  
reinventing the wheel.  Does anyone have a slick and quick way of  
doing this special computation?