overcompensated is incorrect

[claudeha]

I'm not sure where the problem lies, but I initially puzzled over incorrect rendering in a port to another language, before testing the upstream (this package).

Tested with a small Mandelbrot set renderer. Until zoom depth 1e30 or so images look fine, but then they drift to feasible-looking images centred at the wrong point (eventually becoming completely wrong):

Incorrect image using Compensated (Compensated Double):

Correct image rendered using Numeric.QD.QuadDouble:

source code for test program:

-- A simple Mandelbrot set rendering test
-- shows Compensated (Compensated Double) is incorrect.
-- Swap the three commented lines out for the following ones
-- to compare with the correct `qd` package (bindings to C libqd).
-- The coordinates of the center at zoom depth 1e100 (or so)
-- is the artwork "Polefcra" by Jonathan Leavitt.
-- Usage: ./mandelbrot 1e32 > 1e32.pgm && display 1e32.pgm
-- The emblem should be centered in the view.

import Data.ByteString.Char8
import Data.Char
import System.Environment (getArgs)

-- import Numeric.QD.QuadDouble
import Numeric.Compensated

-- type R = QuadDouble
type R = Compensated (Compensated Double)

data Complex a = Complex !a !a
  deriving (Show)

instance Num a => Num (Complex a) where
  Complex a b + Complex c d = Complex (a + c) (b + d)
  Complex a b - Complex c d = Complex (a - c) (b - d)
  Complex a b * Complex c d = Complex (a * c - b * d) (a * d + b * c)
  fromInteger a = Complex (fromInteger a) 0

norm (Complex a b) = a * a + b * b

data Dual a = Dual !a !a
  deriving (Show)

instance Num a => Num (Dual a) where
  Dual a b + Dual c d = Dual (a + c) (b + d)
  Dual a b - Dual c d = Dual (a - c) (b - d)
  Dual a b * Dual c d = Dual (a * c) (a * d + b * c)
  fromInteger a = Dual (fromInteger a) 0

r2 = 100
n = 10100

mandelbrot c i z@(Dual zz zdz)
  | i < n && norm zz < r2 = mandelbrot c (i + 1) (z * z + c)
  | i < n = let z2 = norm zz in sqrt (z2 / norm zdz) * log z2
  | otherwise = 0

c0 :: Dual (Complex R)
c0 = Dual (Complex (-1.74969984376725522654321799095190427951610093088539085379101021670246103474736800872666844130359169734068) 9.25412893997304893805990893994351884853682218769826224511888835093933795525097546683325515581321148e-7) 0
w = 64
h = 36

image :: R -> ByteString
image zoom = pack ("P5\n" ++ show w ++ " " ++ show h ++ "\n255\n" ++ [ m (c0 + Dual (Complex (fromInteger i + 0.5 - fromInteger w / 2) (fromInteger j + 0.5 - fromInteger h / 2)) 1 * Dual (Complex (1 / (fromInteger h * zoom)) 0) 0) | j <- [0..h-1], i <- [0..w-1] ])
 where
  m c = case mandelbrot c 0 0 of
    d | d > 0 -> let k = 255 + 16 * log d in chr (round k `mod` 256)
      | otherwise -> chr 0

main = do
  [sZoom] <- getArgs
  -- let zoom = fromDouble (read sZoom)
  let zoom = compensated (compensated (read sZoom) 0) 0
  Data.ByteString.Char8.putStr $ image zoom

[claudeha]

Problem is in (*) or one of its dependencies

let a = compensated (compensated (-9.231236890307623e12::Double) 6.897617722007686e-4) (compensated (-2.6232210439034232e-20) (-9.294634128446337e-37)) ; b = compensated (compensated (-9.060527246314555e-31) 4.540998182136166e-47) (compensated (-3.1624962852246887e-63) (3.469844214603107e-80 :: Double)) ; rnd = toRational (a * b) ; ext = toRational  a* toRational b in logBase 2 (fromRational $ abs (rnd - ext) /  abs ext)

Prints -108.26948516958842 instead of about -215. Test case found by QuickCheck. (+) passed tests.

[claudeha]

Problem is in times or one of its dependencies (but split passes rudimentary tests).

main = do
  let a,b :: Compensated Double
      a = compensated (-211.11871471905192) 8.926709967512482e-15
      b = compensated 3.9220612083458963e-19 (-1.80497973320049e-35)
      rounded = toRational (times a b compensated)
      exact = toRational a * toRational b
      err = logBase 2 . fromRational $ abs (rounded - exact) / abs exact
  print err

Prints -109 instead of ~-215

[claudeha]

Making split be based on with fixes the last failure case (those numbers now print -Infinity), but there are still failures in times:

=== prop_times from tests/quickcheck.hs:80 ===
*** Failed! Falsified (after 158 tests):
compensated 0.13127052166409675 (-5.807852802487686e-20)                     
compensated (-2.730255721631754e-2) (-9.105471357507769e-19)

I think stacking compensation is ill-founded because add can give results that differ greatly in magnitude, and then times of such compensateds gives results with too many significant bits to fit into just two parts: truncating the least significant bits of the most significant half matches the behavior I saw initially..

[cartazio]

Roughly speaking: do we need to have a sortah vector rep then? I actually started some experiments Towards that for some of my own stuff?

[claudeha]

I don't know. The result of times needs to be exact afaict, otherwise stacked (*) breaks...

In the worst case I think that means times (compensated a b) (compensated c d) needs 8 limbs (instead of just four with compensated (compensated ...)) to store the results (if a,b,c,d are all widely spaced in magnitude, then a*c,a*d,b*c,b*d might not overlap, and they need to be two limbs each to be exactly representable...).

[claudeha]

An alternative idea I haven't tried is to truncate or round the least significant half, so that the total precision is fixed - the cost may be prohibitive though. Something like this perhaps?

normalizeAfterAdd hi lo f = f hi lo'
  where
    p = copysign ((nextafter hi (2*hi) - hi) * 0.5) lo
    lo' = (lo + p) - p

This would bring it in line with RealFloat, too.