Subs.hs

-- | We need to find (small) integer substitution such
-- that the denominator in our formula never vanishes.
--
-- That is, for @n@ we need to find @n@ integers @a[i]@ such that:
--
-- * @0  /=  a[i] -  a[j]@
--   
-- * @0  /=  a[i] - (a[j] + a[k])@
--
-- * @0  /=  (a[i] + a[j]) - (a[k] + a[l])@
--

{-# LANGUAGE ScopedTypeVariables, BangPatterns #-}
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

--------------------------------------------------------------------------------
-- * Lazily cached tables of substitutions

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)

-- | We cache a substitution table
theSubsTable :: [Array Int Integer]
theSubsTable = [ listArray (1,n) (Unsafe.unsafePerformIO (findSubs n)) | n <- [0..] ]

--------------------------------------------------------------------------------
-- * Find substitutions

-- | Select @k@ elements from a list in all possible orders
choosePerm :: Int -> [a] -> [[a]]
choosePerm n xs = concatMap Data.List.permutations (choose n xs)

-- | Checks if a substitution satisfies the constraints
{-# SPECIALIZE checkSubs :: [Int]     -> Bool #-}
{-# SPECIALIZE checkSubs :: [Integer] -> Bool #-}
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

-- | Find random substitution which satisfies the constraints
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 (cnt-1) bound
      True  -> do
        putStrLn $ "good substitution found! " ++ show subs
        return subs

--------------------------------------------------------------------------------