-- | 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@).

{-# LANGUAGE DataKinds #-}
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.Algebra.Polynomial.FreeModule as ZMod

import Math.Algebra.Polynomial.Univariate
import Math.Algebra.Polynomial.Univariate.Lagrange

import Math.RootLoci.Algebra
import Math.RootLoci.Classic

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

type X = U "x"

mkX :: Int -> X
mkX = U

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

-- | 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 (U 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 = unUni $ 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 ++ " ]"
-}