module Math.RootLoci.Dual.Localization where
import Control.Monad
import Data.List
import Data.Ratio
import Math.Combinat.Numbers
import Math.Combinat.Sets
import Math.Combinat.Sign
import Math.Combinat.Partitions
import qualified Data.Map as Map
import qualified Math.RootLoci.Algebra.FreeMod as ZMod
import Math.RootLoci.Algebra
import Math.RootLoci.Classic
localizeMathematica :: Partition -> String
localizeMathematica (Partition xs) = formula where
n = sum xs
d = length xs
ies = [ (head ys, length ys) | ys <- group (sort xs) ]
es = map snd ies
paren str = '(' : str ++ ")"
wt j = paren $ show j ++ "a+" ++ show (nj) ++ "b"
sumOver = listTensor [ [0..e] | e<-es ]
formula = global ++ " * " ++ paren (intercalate " + " (map term sumOver))
global = intercalate "*" [ wt j | j<-[0..n] ] ++ " / (b-a)^" ++ show d
rkonst ss = product [ factorial s * factorial (es) | (s,e) <- zip ss es ]
konst ss = show (paritySignValue (sum ss))
++ "/" ++ show (rkonst ss)
denom ss = show n ++ "*a - " ++ show (sum [ i*s | ((i,e),s) <- zip ies ss ]) ++ "*(a-b)"
term ss = konst ss ++ " / " ++ paren (denom ss)
localizeEval :: Fractional q => Partition -> q -> q -> q
localizeEval (Partition xs) a b = formula where
n = sum xs
d = length xs
ies = [ (head ys, length ys) | ys <- group (sort xs) ]
es = map snd ies
wt j = fromIntegral j * a + fromIntegral (nj) * b
sumOver = listTensor [ [0..e] | e<-es ]
formula = global * sum (map term sumOver)
global = product [ wt j | j<-[0..n] ] / (ba)^d
rkonst ss = product [ factorial s * factorial (es) | (s,e) <- zip ss es ]
konst ss = fromIntegral (paritySignValue (sum ss))
/ fromIntegral (rkonst ss)
denom ss = fromIntegral n * a fromIntegral (sum [ i*s | ((i,e),s) <- zip ies ss ]) * (ab)
term ss = konst ss / denom ss
localizeDual :: Partition -> ZMod AB
localizeDual part = ZMod.mapBase worker $ localizeInterpolatedZ part where
c = codim part
worker (X i) = AB (ci) i
localizeInterpolatedQ :: Partition -> QMod X
localizeInterpolatedQ part@(Partition xs) = final where
codim = sum xs length xs
bs = map fromIntegral [ 2..codim+2 ] :: [Rational]
qs = [ localizeEval part 1 b | b<-bs ] :: [Rational]
final = lagrangeInterp' (zip bs qs)
localizeInterpolatedZ :: Partition -> ZMod X
localizeInterpolatedZ = (ZMod.mapCoeff f . localizeInterpolatedQ) where
f :: Rational -> Integer
f q = case denominator q of
1 -> numerator q
_ -> error "non-integral coefficient in the result"