numeric basis issues

[barak]

This issue is sort of a placeholder, but it is hard for a plugin like this to operate in an environment where the numeric basis itself is rife with unnecessary loss of precision. Below are two examples of extremely poor and unnecessary truncation, of small values to zero, in atanh on Double and in asin and atan on Complex Double. It would be great if someone would let this tool run wild on the prelude!

$ ghci
GHCi, version 7.8.4: http://www.haskell.org/ghc/  :? for help
Loading package ghc-prim ... linking ... done.
Loading package integer-gmp ... linking ... done.
Loading package base ... linking ... done.
Prelude> tan 5e-17
5.0e-17
Prelude> atan 5e-17
5.0e-17
Prelude> tanh 5e-17
5.0e-17
Prelude> atanh 5e-17
0.0
Prelude> :m + Data.Complex
Prelude Data.Complex> sin (0 :+ 5e-17)
0.0 :+ 5.0e-17
Prelude Data.Complex> asin (0 :+ 5e-17)
0.0 :+ (-0.0)
Prelude Data.Complex> tan (0 :+ 5e-17)
0.0 :+ 5.0e-17
Prelude Data.Complex> atan (0 :+ 5e-17)
0.0 :+ (-0.0)

[mikeizbicki]

FWIW, the standard prelude calculates

atanh x = 0.5 * (log ((1 + x) / (1 - x)))

and herbie uses a taylor expansion to improve it to

atanh x =  0.5 * (((2 / 5) * x * x * x * x * x) + ((2 / 3) * x * x * x) + (2 * x))

Herbie reports the original expression has 58 bits of loss, and the improved has only 8e-2 bits of loss. This is probably worth reporting as a GHC bug.

I imagine these functions are used often enough there are pretty standard implementations for them. I'm actually surprised they aren't implemented as primops within GHC.

[barak]

Agreed. These should all use standard numerically safe expressions. And Double etc should call out to the math library.

Even functional language implementations that allegedly cater to numeric accuracy, like MIT Scheme, often get dramatic loss of precision on atan of a small imaginary number. I guess people are getting a bit tired of filling in all these formulas by the time they get to atan of a complex number.

PS Shouldn't

atanh x =  0.5 * (((2 / 5) * x * x * x * x * x) + ((2 / 3) * x * x * x) + (2 * x))

be

atanh x =  ((1 / 5) * x * x * x * x * x) + ((1 / 3) * x * x * x) + x

or am I missing something?

[mikeizbicki]

I think you're version is better. It's just not what Herbie reports.

[barak]

Not to belabor this too much, but there is a reason numeric libraries include log1p and friends, and the definition of atanh is a perfect example. If we have that primitive, more precise definitions practically write themselves:

atanh x = 0.5 * log ((1 + x) / (1 - x))
        = 0.5 * (log(1+x) - log(1-x))
        = 0.5 * (log1p(x) - log1p(-x))

or if we'd like to avoid two invocations of a special function,

atanh x = 0.5 * log ((1 + x) / (1 - x))
        = 0.5 * log ((1-x)/(1-x) + 2*x/(1-x))
        = 0.5 * log1p(2*x/(1-x))

or just for fun,

        = - log1p(-x) + 0.5*log1p(-x^2)

[mikeizbicki]

Of course! I would very much like to get these functions supported because
they are very important.

On Wed, Sep 30, 2015 at 3:02 PM, Barak A. Pearlmutter <
notifications@github.com> wrote:

Not to belabor this too much, but there is a reason numeric libraries
include log1p and friends, and the definition of atanh is a perfect
example. If we have that primitive, more precise definitions practically
writes themselves:

atanh x = 0.5 * log ((1 + x) / (1 - x))
= 0.5 * (log(1+x) - log(1-x))
= 0.5 * (log1p(x) - log1p(-x))

or if we'd like to avoid two invocations of a special function,

atanh x = 0.5 * log ((1 + x) / (1 - x))
= 0.5 * log ((1-x)/(1-x) + 2_x/(1-x))
= 0.5 * log1p(2_x/(1-x))

or just for fun,

    = - log1p(-x) + 0.5*log1p(-x^2)

—
Reply to this email directly or view it on GitHub
#15 (comment)
.