-- | Localization formula for the dual class from:
--
-- L. M. Feher, A. Nemethi, R. Rimanyi: Coincident root loci of binary forms;
-- Michigan Math. J. Volume 54, Issue 2 (2006), 375--392.
--
-- Note: This formula is in the form of /rational function/ (as opposed to 
-- a polynomial). Since we don't have polynomial factorization implemented here,
-- instead we /evaluate/ it substituting different rational numbers
-- into @alpha@ and @beta@, and then use Lagrange interpolation to figure
-- out the result (we know a priori that it is a homogenenous polynomial
-- in @alpha@ and @beta@).

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

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

-- | The localization formula as a string which Mathematica can parse
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 (n-j) ++ "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 (e-s) | (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)

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

-- | The localization formula evaluated at given values of @a@ and @b@
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 (n-j) * b

  sumOver = listTensor [ [0..e] | e<-es ] 
  formula = global * sum (map term sumOver)

  global = product [ wt j | j<-[0..n] ] / (b-a)^d

  rkonst ss = product [ factorial s * factorial (e-s) | (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 ]) * (a-b)
  term   ss = konst ss / denom ss

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

-- | The dual class, recovered via from the localization formula via Lagrange
-- interpolation
localizeDual :: Partition -> ZMod AB
localizeDual part = ZMod.mapBase worker $ localizeInterpolatedZ part where
  c = codim part
  worker (X i) = AB (c-i) i 

-- | We can use Lagrange interpolation to express the dual class from the
-- localization formula (since we know a priori that the result is a homogeneous
-- polynomial in @a@ and @b@)
--
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"

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

{-
main = do
  forM_ (partitions 9) $ \part@(Partition xs) -> do
    putStrLn $ "X" ++ concatMap show xs ++ " = Factor[ " ++ tp_local_mathematica part ++ " ]"
-}