[issue31630] math.tan has poor accuracy near pi/2 on OpenBSD

[Tim Peters]

Tim Peters <tim at python.org> added the comment:

When Sun was developing fdlibm, I was (among other things) working on a proprietary libm for Kendall Square Research.  I corresponded with fdlibm's primary author (KC Ng) often at the time.  There's no way he would have left errors this egregious slip by, or neglect careful testing near singularities.

But writing portable code at this level is very delicate, and fdlibm failed in many ways when I first compiled it for the KSR machine.  That was the first 64-bit box fdlibm was tried on, & it exposed many unwarranted assumptions in the code.

Browsing the OpenBSD port, they bumped into more, seemingly driven by gcc violating assumptions about how tricky casting would work - there are many places in the code that, e.g., want to view a single double as a pair of int32s (or vice versa), at the bit level.  An aggressive compiler can easily be provoked into reordering operations in fatally damaging ways, if it doesn't realize that, under the maze of cast and union tricks, successive operations are actually operating on the _same_ bits.

So, if possible, it would be great to try compiling OpenBSD's libm with every conceivable optimization disabled, to see if that changes results.

Why I suspect that may be promising:  as Stefan pointed out, for arguments near pi/2 fdlibm's tan() subtracts pi/2 and evaluates the result as tan(x) = -1/tan(x - pi/2).  But these arguments are _so_ close to pi/2 that |x - pi/2| is "tiny".  But for tiny arguments z, tan(z) == z to machine precision (tan(z) = sin(z) / cos(z) ~= z/1 = z for tiny |z|).  So it's effectively using tan(x) = -1/(x - pi/2) in this range.

In the original failing, test, behold:

>>> (1.0 / 1978937966095219.0).hex() # correct
'0x1.234c4c6628b81p-51'

>>> (1.0 / 1978945885716843.0).hex() # OpenBSD
'0x1.234c000000000p-51'

It's very hard to believe all those trailing zero bits in the OpenBSD reciprocal are a coincidence.  It's extremely hard to believe that fdlibm _intended_ them to be all zeroes.  But it's pretty easy to believe an aggressive optimizer managed to lose the intent of the code.

----------

_______________________________________
Python tracker <report at bugs.python.org>
<https://bugs.python.org/issue31630>
_______________________________________