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 (\n -> [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 d :: [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 (\u -> 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 h :: [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