DoublePermsSum.hs

-- A Haskell based SMT solution attempt to:
-- https://stackoverflow.com/questions/69192991/how-many-different-sums-can-we-get-from-very-few-floats

import Data.List
import Control.Monad

import Data.SBV

-- Make a permutation from 0 to k-1
mkPerm :: Word8 -> Word8 -> Symbolic [SWord8]
mkPerm k i = do
   ps <- mapM (\ind -> free ("perm_" ++ show i ++ "_" ++ show ind)) [1 .. k]
   constrain $ distinct ps
   mapM_ (\p -> constrain (p .< literal k)) ps
   pure ps

add :: Word8 -> [SDouble] -> [SWord8] -> SDouble
add n ds ps = walk ps 0
 where walk []     sofar = sofar
       walk (p:ps) sofar = walk ps (pick ds p + sofar)

       pick []     _ = error "impossible"
       pick [d]    _ = d
       pick (d:ds) p = ite (p .== 0) d (pick ds (p-1))

-- Given n doubles, are there k permutations of them
-- such that each sum to a different value
gen :: Word8 -> Word8 -> Symbolic ()
gen n k = do
   ds <- mapM (\i -> free ("v" ++ show (i-1))) [1..n]
   mapM_ (constrain . fpIsPoint) ds

   pps  <- mapM (mkPerm n) [1..k]

   let sums = map (add n ds) pps
   constrain $ distinct sums

disp :: Maybe ([Double], [Word8]) -> IO ()
disp Nothing          = putStrLn "No satisfying model found"
disp (Just (xs, pps)) = do
   zipWithM_ (\i x -> putStrLn ("v" ++ show i ++ " = " ++ show x)) [0..] xs
   let sp ps = intercalate " + " ["v" ++ show p | p <- ps] ++ " = " ++ show (sum [xs !! fromIntegral p | p <- ps])
       by _ [] = []
       by n xs = take n xs : by n (drop n xs)
   mapM_ (putStrLn . sp) (by (length xs) pps)

-- Find 4 doubles with 4 combinations summing to different values
-- This takes about 25 seconds on my machine to find the following:
--
-- v0 = 1.2499704954752813
-- v1 = -0.24998477360168156
-- v2 = -2.1151680754116775e-10
-- v3 = -0.9999857218953875
-- v2 + v1 + v3 + v0 = -2.3330448684077965e-10
-- v3 + v0 + v1 + v2 = -2.3330457357695344e-10
-- v1 + v3 + v0 + v2 = -2.333044903102266e-10
-- v2 + v0 + v1 + v3 = -2.333045978630821e-10
easy :: IO ()
easy = (disp . extractModel) =<< sat (gen 4 4)

-- Find 4 doubles with 10 combinations summing to different values
-- I ran this for more than 2 hours, but it did not end.
hard :: IO ()
hard = (disp . extractModel) =<< sat (gen 4 10)