Skip to content

Instantly share code, notes, and snippets.

@LeventErkok
Created September 19, 2021 04:35
Show Gist options
  • Save LeventErkok/75d0c0bad4a4415ca26d0a5a7eda95af to your computer and use it in GitHub Desktop.
Save LeventErkok/75d0c0bad4a4415ca26d0a5a7eda95af to your computer and use it in GitHub Desktop.
-- 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)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment