1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
|
{-# OPTIONS_GHC -fno-warn-orphans #-}
module Tests.Sum (tests) where
import Control.Applicative ((<$>))
import Numeric.Sum as Sum
import Numeric.MathFunctions.Comparison
import Prelude hiding (sum)
import Test.Tasty (TestTree, testGroup)
import Test.Tasty.QuickCheck
import Test.QuickCheck (Arbitrary(..))
import qualified Prelude
-- Test that summation result is same as exact sum. That should pass
-- if we're effectively working with quad precision
t_sum :: ([Double] -> Double) -> [Double] -> Property
t_sum f xs
= counterexample ("APPROX = " ++ show approx)
$ counterexample ("EXACT = " ++ show exact)
$ counterexample ("DELTA = " ++ show (approx - exact))
$ counterexample ("ULPS = " ++ show (ulpDistance approx exact))
$ approx == exact
where
approx = f xs
exact = trueSum xs
-- Test that summation has smaller error than naive summation or no
-- worse than given number of ulps. If we're close enough to exact
-- answer naive may get ahead
t_sum_error :: ([Double] -> Double) -> [Double] -> Property
t_sum_error f xs
= counterexample ("APPROX = " ++ show approx)
$ counterexample ("NAIVE = " ++ show naive)
$ counterexample ("EXACT = " ++ show exact)
$ counterexample ("A-EXACT = " ++ show (approx - exact))
$ counterexample ("N-EXACT = " ++ show (naive - exact))
$ counterexample ("ULPS[A] = " ++ show (ulpDistance approx exact))
$ counterexample ("ULPS[N] = " ++ show (ulpDistance naive exact))
$ abs (exact - approx) <= abs (exact - naive)
where
naive = Prelude.sum xs
approx = f xs
exact = trueSum xs
t_sum_shifted :: ([Double] -> Double) -> [Double] -> Property
t_sum_shifted f = t_sum_error f . zipWith (+) badvec
trueSum :: (Fractional b, Real a) => [a] -> b
trueSum xs = fromRational . Prelude.sum . map toRational $ xs
badvec :: [Double]
badvec = cycle [1, 1e14, -1e14]
tests :: TestTree
tests = testGroup "Summation"
[ testGroup "Kahan" [
-- Kahan summation only beats naive summation when truly
-- catastrophic cancellation occurs
testProperty "t_sum_shifted" $ t_sum_shifted (sum kahan)
]
, testGroup "KBN" [
testProperty "t_sum" $ t_sum (sum kbn)
, testProperty "t_sum_error" $ t_sum_error (sum kbn)
, testProperty "t_sum_shifted" $ t_sum_shifted (sum kbn)
]
, testGroup "KB2" [
testProperty "t_sum" $ t_sum (sum kb2)
, testProperty "t_sum_error" $ t_sum_error (sum kb2)
, testProperty "t_sum_shifted" $ t_sum_shifted (sum kb2)
]
]
instance Arbitrary KahanSum where
arbitrary = toKahan <$> arbitrary
shrink = map toKahan . shrink . fromKahan
toKahan :: (Double, Double) -> KahanSum
toKahan (a,b) = KahanSum a b
fromKahan :: KahanSum -> (Double, Double)
fromKahan (KahanSum a b) = (a,b)
instance Arbitrary KBNSum where
arbitrary = toKBN <$> arbitrary
shrink = map toKBN . shrink . fromKBN
toKBN :: (Double, Double) -> KBNSum
toKBN (a,b) = KBNSum a b
fromKBN :: KBNSum -> (Double, Double)
fromKBN (KBNSum a b) = (a,b)
instance Arbitrary KB2Sum where
arbitrary = toKB2 <$> arbitrary
shrink = map toKB2 . shrink . fromKB2
toKB2 :: (Double, Double, Double) -> KB2Sum
toKB2 (a,b,c) = KB2Sum a b c
fromKB2 :: KB2Sum -> (Double, Double, Double)
fromKB2 (KB2Sum a b c) = (a,b,c)
|
