Ah!  So it looks like OpenBSD took its math functions from Sun's fdlibm.  I believe wxMaxima does too.  That would certainly explain why they give the same answer ;-)

So who's right?  Using "bigfloats" in wxMaxima and feeding its bigfloat tan() the exact decimal value corresponding to the binary double denoted by 1.5707963267948961, it reproduces the result the test is expecting.  Ditto by writing my own tan function, using `decimal` with 50 digits, using the straightforward series expansions for sin and cos, and then dividing to get tan (slow, but without any attempt at "cleverness" it's easy to trust it).

So if fdlibm's author were still active, I'm sure he'd consider this to be a bug - it _intended_ to have < 1 ulp error in almost all cases.

It's not argument reduction that's hurting it.  fdlibm stores pi/2 to waaaay more bits than are needed to get the right answer via -1/tan(x - pi/2).  

Regardless, there's nothing Python can do about it short of changing the test.  As the Python results I pasted before show, the OpenBSD answer _is_ correct for some (non-representable in double precision) input value within 1 ulp of the actual input value.  So it's a very demanding test to get right.