module Math.ThomPoly.Subs where
import Data.Array
import Data.List
import Control.Monad
import System.Random
import System.IO.Unsafe as Unsafe
import Math.Combinat.Sets
import Math.Algebra.ModP
getSubs :: Int -> Array Int Integer
getSubs n = theSubsTable!!n
getSubsNum :: Num a => Int -> Array Int a
getSubsNum n = fmap fromInteger (theSubsTable!!n)
getSubsZp :: Int -> Array Int Zp
getSubsZp n = fmap mkZp (theSubsTable!!n)
theSubsTable :: [Array Int Integer]
theSubsTable = [ listArray (1,n) (Unsafe.unsafePerformIO (findSubs n)) | n <- [0..] ]
choosePerm :: Int -> [a] -> [[a]]
choosePerm n xs = concatMap Data.List.permutations (choose n xs)
checkSubs :: forall a. (Eq a, Num a) => [a] -> Bool
checkSubs input = ok2 && ok3 && ok4 where
n = length input
nn = [1..n]
arr = listArray (1,n) input :: Array Int a
ok2 = and [ 0 /= (arr!i arr!j) | [i,j] <- choose 2 nn ]
ok3 = and [ 0 /= (arr!i arr!j arr!k) | [i,j,k] <- choosePerm 3 nn ] &&
and [ 0 /= (arr!i 2 * arr!j ) | [i,j] <- choosePerm 2 nn ]
ok4 = and [ 0 /= (arr!i + arr!j arr!k arr!l) | [i,j,k,l] <- choosePerm 4 nn ] &&
and [ 0 /= (arr!i + arr!j 2 * arr!k ) | [i,j,k] <- choosePerm 3 nn ]
findSubsZp :: Int -> IO [Zp]
findSubsZp = liftM (map mkZp) . findSubs
findSubs :: Int -> IO [Integer]
findSubs n = go 25 where
go !bound = tryN 100 bound
tryN 0 !bound = go (div (bound*3) 2)
tryN !cnt !bound = do
subs <- replicateM n $ randomRIO (bound,bound)
case checkSubs subs of
False -> tryN (cnt1) bound
True -> do
putStrLn $ "good substitution found! " ++ show subs
return subs