1 -- References are to
2 --     http://arxiv.org/abs/0912.1436
3 -- or
4 --     http://zeros.spag.dk
5 import List
6
7 -- Note the ordering of the lists:
8 -- [i_m, ..., i_1], [s_m, ..., s_1], [u_1, ..., u_r]
9
10 -- sumMult [a_1, ..., a_n] = \sum_{i=1}^n i a_i
11 sumMult = sum . zipWith (*) [1..]
12
13 -- maximize f [a_1, ..., a_n] = max{f(a_1), ..., f(a_n)}
14 maximize f = maximum . map f
15
16 -- Given a list [s1, ..., sn] produce
17 -- [ [x1, ..., xn] | x1 <- [0..s1], ..., xn <- [0..sn] ].
18 cartesian :: [Integer] -> [[Integer]]
19 cartesian = mapM (\-> [0..n])
20
21 cart dim size = cartesian (take (fromInteger dim) (repeat (fromInteger size)))
22
23 -- Check if the i vector is outside the stairs form.
24 -- See Remark 12.
25 outside :: [Integer] -> Integer -> [Integer] -> Bool
26 outside is r ss = r <= sum (zipWith div is ss)
27
28 -- See Definition 10 and Theorem 11.
29 :: [Integer] -> Integer -> [Integer] -> Integer
30 d is r ss | outside is r ss = product ss
31           | otherwise       = dr (reverse is) r (reverse ss)
32     where
33     dr [i] r [s] = min (i div r) s
34     dr (i:is) r (s:ss) = maximize f (setA i r s)
35         where
36             f us = (s - sum us) * (dr is r ss)
37                    + summing us is r ss
38             -- Note that [u_1, ..., u_r].
39             summing [u] is 1 ss = u * product ss
40             summing (u:us) is r ss = u * (dr is (r-1) ss)
41                                      + summing us is (r-1) ss
42
43 -- See Definition 10.
44 setA :: Integer -> Integer -> Integer -> [[Integer]]
45 setA i r s = filter (\-> sum u <= s && sumMult u <= i)
46                     (cart r s)
47
48 -- d' [i_1, ..., i_m] r q = d [i_1, ..., i_m] r [q, ..., q]
49 d' :: [Integer] -> Integer -> Integer -> Integer
50 d' is r q = d is r (replicate (length is) q)
51
52 -- See Proposition 39.
53 htilde :: [[Integer]] -> Integer -> [Integer] -> Integer
54 htilde ws k ss = hr (reverse ws) k (reverse ss)
55     where
56     -- [s_m, ..., s_2, s_1]
57     -- [w_m, ..., w_2, w_1] = [[v^m_1, ..., v^m_r], ..., [v^1_1, ..., v^1_r]]
58     hr [vs] k [s] = sum $drop (fromIntegral (k-1)) vs 59 hr (w:ws) k (s:ss) = (s - (sum w)) * (hr ws k ss) 60 + summing (take (fromIntegral (k-1)) w) k 61 + (sum (drop (fromIntegral (k-1)) w)) * (product ss) 62 where 63 summing [] 1 = 0 64 summing (v:vs) k' = v * (hr ws (k'-1) ss) 65 + summing vs (k'-1) 66 hr [] k s = fromIntegral (40 + length s) 67 68 -- Lower Bound, H. See Proposition 40. 69 :: [Integer] -> Integer -> [Integer] -> Integer 70 h is r ss = maximize (\ws -> htilde ws r ss) (sequence (set is ss)) 71 where 72 set [] [] = [] 73 set (i:is) (s:ss) = 74 (filter (\vs -> sum vs <= s && sumMult vs == i) (cart r s)) 75 : set is ss 76 77 -- h' [i_1, ..., i_m] r q = h [i_1, ..., i_m] r [q, ..., q] 78 h' :: [Integer] -> Integer -> Integer -> Integer 79 h' is r q = h is r (replicate (length is) q) 80 81 -- See Corollary 9. 82 -- min { (i_1·s_1···s_m + s_1·i_2·s_3···s_m + ··· + s_1···s_{m-1}·i_m) / r 83 -- , s_1···s_m } 84 sz :: [Integer] -> Integer -> [Integer] -> Integer 85 sz is r ss = min cor9 (product ss) 86 where 87 cor9 = (sum (zipWith (*) is (map product$ f ss))) div r
88         f [] = []
89         f (x:xs) = xs : map (x:) (f xs)
90
91 -- sz' [i_1, ..., i_m] r q
92 -- = min{floor((i_1 + ... + i_m) * q^(m-1) / r), q^m}
93 -- See (2.1).
94 sz' :: [Integer] -> Integer -> Integer -> Integer
95 sz' is r q = min (((sum is) * q^(m-1)) div r) (q^m)
96     where m = length is
97
98 -- Closed Formula; two variables.
99 cf :: [Integer] -> Integer -> [Integer] -> Integer
100 cf is r ss = truncate (cfrat is r ss)
101
102 cf' :: [Integer] -> Integer -> Integer -> Integer
103 cf' is r q = truncate (cfrat is r [q,q])
104
105 -- See Proposition 19.
106 cfrat :: [Integer] -> Integer -> [Integer] -> Rational
107 cfrat [i_i1,i_i2] i_r [i_s1,i_s2] =
108     minimum (l c1 ++ l c2 ++ l c3 ++ [c4])
109     where
110         i1 = fromIntegral i_i1; i2 = fromIntegral i_i2
111         r = fromIntegral i_r
112         s1 = fromIntegral i_s1; s2 = fromIntegral i_s2
113         l c = map c [1..r-1]
114         c1 k = if (r-k)*(r/(r+1))*s1 <= i1 && i1 < (r-k)*s1
115                   && 0 <= i2 && i2 < k*s2
116                    then s2*(i1/r) + (i2/r)*(i1/(r-k))
117                    else s1*s2
118         c2 k = if (r-k)*(r/(r+1))*s1 <= i1 && i1 < (r-k)*s1
119                   && k*s2 <= i2 && i2 < (k+1)*s2
120                    then s2*(i1/r)
121                         + ((k+1)*s2 - i2) * ((i1/(r-k))-(i1/r))
122                         + (i2-k*s2)*(s1-(i1/r))
123                    else s1*s2
124         c3 k = if (r-k-1)*s1 <= i1 && i1 < (r-k)*(r/(r+1))*s1
125                   && 0 <= i2 && i2 < (k+1)*s2
126                    then s2*(i1/r) + (i2/(k+1))*(s1-(i1/r))
127                    else s1*s2
128         c4 = if s1*(r-1<= i1 && i1 < s1*r
129                 && 0 <= i2 && i2 < s2
130                  then s2*(fromIntegral (floor (i1/r)))
131                       + i2*(s1 - fromIntegral ( floor (i1/r)))
132                  else s1*s2