cartazio / numbers

Other
29 stars 13 forks source link

Data.Number.Fixed: loss of accuracy in exp() function #5

Closed selinger closed 10 years ago

selinger commented 11 years ago

The exp(x) function in Data.Number.Fixed does not have anywhere near the advertised accuracy. This is due to the fact that the implementation multiplies very large numbers by limited-precision numbers, therefore amplifying the round-off errors.

To reproduce:

Prelude Data.Number.Fixed> exp 100 :: Fixed Prec50 26881171418161354484126255515800135873611118.77373326422066361293765446985616509052697695397433

Compare:

Prelude Data.Number.Fixed> exp 100 :: Fixed (PrecPlus20 (PrecPlus20 (PrecPlus20 Prec50))) 26881171418161354484126255515800135873611118.77374192241519160861528028703490956491415887109721984571081167088634206600370770528301079734543724490162990550 Prelude Data.Number.Fixed> :m Data.Number.CReal Prelude Data.Number.CReal> putStrLn $ showCReal 50 $ exp 100 26881171418161354484126255515800135873611118.77374192241519160861528028703490956491415887109722

As the example shows, the "Fixed" implementation loses 47 digits of accuracy compared to the true value. This is approximately equal to the number of digits before the decimal point.

To fix: temporarily increase the accuracy during the computation.

Add the auxiliary function:

-- | The call @with_added_precision r f v@ evaluates @f v@, while
-- temporarily multiplying the precision of /v/ by /r/.
with_added_precision :: forall a f.(Epsilon f) => Rational -> (forall e.(Epsilon e) => Fixed e -> a) -> Fixed f -> a
with_added_precision r f v = dynamicEps (p*r) f (toRational v) where
  p = precision v

And replace

    exp = toFixed1 F.exp

by:

    exp x = with_added_precision r (convertFixed . (toFixed1 F.exp)) x where
      r = if x < 0 then 1 else 0.1 ^ (ceiling (x * 0.45))

Note: the issue also affects cosh(x), sinh(x), and potentially tanh(x), but only because the latter are defined in terms of exp(x).