[Issue 15365] New: std.math: 80-bit exp() tests are utterly wrong when returning subnormals

[via Digitalmars-d-bugs]

https://issues.dlang.org/show_bug.cgi?id=15365

          Issue ID: 15365
           Summary: std.math: 80-bit exp() tests are utterly wrong when
                    returning subnormals
           Product: D
           Version: D2
          Hardware: x86_64
                OS: Linux
            Status: NEW
          Severity: enhancement
          Priority: P1
         Component: phobos
          Assignee: nobody at puremagic.com
          Reporter: ibuclaw at gdcproject.org

When comparing the DMD-style inline assembler path vs. the generic cephes
implementation, observed the following.

These test points:

(1):
[ 0x1.1p+13L,     0x1.29aeffefc8ec645p+12557L  ], // near overflow
http://www.wolframalpha.com/input/?i=exp%288704%29
1.256523125474349284688509465050423995871881189122581533e3780L

feqrel(cephes, iasm) = 52
feqrel(iasm, wolfram) = 52
feqrel(cephes, wolfram) = 64

=========

(2):
[-0x1.18p+13L,    0x1.5e4bf54b4806db9p-12927L  ], // near underflow
http://www.wolframalpha.com/input/?i=exp%28-8960%29
5.26553067155989234457558352116540147577224389983906948e-3892L

feqrel(cephes, iasm) = 51
feqrel(iasm, wolfram) = 51
feqrel(cephes, wolfram) = 64

=========

(3):
[-0x1.625p+13L,   0x1.a6bd68a39d11f35cp-16358L ], // ditto
http://www.wolframalpha.com/input/?i=exp%28-11338%29
9.31459938556868263111811562594117020130237890448547800e-4925L

feqrel(cephes, iasm) = 53
feqrel(iasm, wolfram) = 53
feqrel(cephes, wolfram) = 64

=========

(4):
[-0x1.62dafp+13L, 0x1.96c53d30277021dp-16383L  ], // near underflow - subnormal
http://www.wolframalpha.com/input/?i=exp%28-11355.4%29
2.58487886630776605706860595662918121578535663781714980e-4932L

feqrel(cephes, iasm) = 54
feqrel(iasm, wolfram) = 5
feqrel(cephes, wolfram) = 5

=========

In (1), (2), and (3) the generic version is more accurate than the IASM
version, matching bit-for-bit the result given from wolfram.

What is interesting is that in each case, they share up to double precision in
common.  This suggests that whatever the IASM branch does (I assume this
happens in L_subnormal) it truncates or rounds the value to double.


Test (4) is simply way off the mark for both IASM and generic, and should
probably be scrapped as it's clear we can't guarantee any reasonable level of
accuracy with any lower number.  If I understand the documentation correctly,
the domain that the exp() function has been tested with is +-10000.

--