-- | 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 :: Int -> X
mkX = Int -> X
forall (var :: Symbol). Int -> U var
U

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

-- | The localization formula as a string which Mathematica can parse
localizeMathematica :: Partition -> String
localizeMathematica :: Partition -> String
localizeMathematica (Partition [Int]
xs) = String
formula where
  n :: Int
n   = [Int] -> Int
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Int]
xs
  d :: Int
d   = [Int] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Int]
xs
  ies :: [(Int, Int)]
ies = [ ([Int] -> Int
forall a. [a] -> a
head [Int]
ys, [Int] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Int]
ys) | [Int]
ys <- [Int] -> [[Int]]
forall a. Eq a => [a] -> [[a]]
group ([Int] -> [Int]
forall a. Ord a => [a] -> [a]
sort [Int]
xs) ]
  es :: [Int]
es  = ((Int, Int) -> Int) -> [(Int, Int)] -> [Int]
forall a b. (a -> b) -> [a] -> [b]
map (Int, Int) -> Int
forall a b. (a, b) -> b
snd [(Int, Int)]
ies

  paren :: String -> String
paren String
str = Char
'(' Char -> String -> String
forall a. a -> [a] -> [a]
: String
str String -> String -> String
forall a. [a] -> [a] -> [a]
++ String
")"
  wt :: Int -> String
wt Int
j = String -> String
paren (String -> String) -> String -> String
forall a b. (a -> b) -> a -> b
$ Int -> String
forall a. Show a => a -> String
show Int
j String -> String -> String
forall a. [a] -> [a] -> [a]
++ String
"a+" String -> String -> String
forall a. [a] -> [a] -> [a]
++ Int -> String
forall a. Show a => a -> String
show (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
j) String -> String -> String
forall a. [a] -> [a] -> [a]
++ String
"b"

  sumOver :: [[Int]]
sumOver = [[Int]] -> [[Int]]
forall a. [[a]] -> [[a]]
listTensor [ [Int
0..Int
e] | Int
e<-[Int]
es ] 
  formula :: String
formula = String
global String -> String -> String
forall a. [a] -> [a] -> [a]
++ String
" * " String -> String -> String
forall a. [a] -> [a] -> [a]
++ String -> String
paren (String -> [String] -> String
forall a. [a] -> [[a]] -> [a]
intercalate String
" + " (([Int] -> String) -> [[Int]] -> [String]
forall a b. (a -> b) -> [a] -> [b]
map [Int] -> String
term [[Int]]
sumOver)) 

  global :: String
global = String -> [String] -> String
forall a. [a] -> [[a]] -> [a]
intercalate String
"*" [ Int -> String
wt Int
j | Int
j<-[Int
0..Int
n] ] String -> String -> String
forall a. [a] -> [a] -> [a]
++ String
" / (b-a)^" String -> String -> String
forall a. [a] -> [a] -> [a]
++ Int -> String
forall a. Show a => a -> String
show Int
d

  rkonst :: [Int] -> Integer
rkonst [Int]
ss = [Integer] -> Integer
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [ Int -> Integer
forall a. Integral a => a -> Integer
factorial Int
s Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Int -> Integer
forall a. Integral a => a -> Integer
factorial (Int
eInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
s) | (Int
s,Int
e) <- [Int] -> [Int] -> [(Int, Int)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int]
ss [Int]
es ]
  konst :: [Int] -> String
konst  [Int]
ss = Integer -> String
forall a. Show a => a -> String
show (Int -> Integer
forall a. Integral a => a -> Integer
paritySignValue ([Int] -> Int
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Int]
ss)) 
            String -> String -> String
forall a. [a] -> [a] -> [a]
++ String
"/" String -> String -> String
forall a. [a] -> [a] -> [a]
++ Integer -> String
forall a. Show a => a -> String
show ([Int] -> Integer
rkonst [Int]
ss)              
  denom :: [Int] -> String
denom  [Int]
ss = Int -> String
forall a. Show a => a -> String
show Int
n String -> String -> String
forall a. [a] -> [a] -> [a]
++ String
"*a - " String -> String -> String
forall a. [a] -> [a] -> [a]
++ Int -> String
forall a. Show a => a -> String
show ([Int] -> Int
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [ Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int
s | ((Int
i,Int
e),Int
s) <- [(Int, Int)] -> [Int] -> [((Int, Int), Int)]
forall a b. [a] -> [b] -> [(a, b)]
zip [(Int, Int)]
ies [Int]
ss ]) String -> String -> String
forall a. [a] -> [a] -> [a]
++ String
"*(a-b)"
  term :: [Int] -> String
term   [Int]
ss = [Int] -> String
konst [Int]
ss String -> String -> String
forall a. [a] -> [a] -> [a]
++ String
" / " String -> String -> String
forall a. [a] -> [a] -> [a]
++ String -> String
paren ([Int] -> String
denom [Int]
ss)

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

-- | The localization formula evaluated at given values of @a@ and @b@
localizeEval :: Fractional q => Partition -> q -> q -> q
localizeEval :: Partition -> q -> q -> q
localizeEval (Partition [Int]
xs) q
a q
b = q
formula where
  n :: Int
n   = [Int] -> Int
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Int]
xs
  d :: Int
d   = [Int] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Int]
xs
  ies :: [(Int, Int)]
ies = [ ([Int] -> Int
forall a. [a] -> a
head [Int]
ys, [Int] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Int]
ys) | [Int]
ys <- [Int] -> [[Int]]
forall a. Eq a => [a] -> [[a]]
group ([Int] -> [Int]
forall a. Ord a => [a] -> [a]
sort [Int]
xs) ]
  es :: [Int]
es  = ((Int, Int) -> Int) -> [(Int, Int)] -> [Int]
forall a b. (a -> b) -> [a] -> [b]
map (Int, Int) -> Int
forall a b. (a, b) -> b
snd [(Int, Int)]
ies

  wt :: Int -> q
wt Int
j = Int -> q
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
j q -> q -> q
forall a. Num a => a -> a -> a
* q
a q -> q -> q
forall a. Num a => a -> a -> a
+ Int -> q
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
j) q -> q -> q
forall a. Num a => a -> a -> a
* q
b

  sumOver :: [[Int]]
sumOver = [[Int]] -> [[Int]]
forall a. [[a]] -> [[a]]
listTensor [ [Int
0..Int
e] | Int
e<-[Int]
es ] 
  formula :: q
formula = q
global q -> q -> q
forall a. Num a => a -> a -> a
* [q] -> q
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum (([Int] -> q) -> [[Int]] -> [q]
forall a b. (a -> b) -> [a] -> [b]
map [Int] -> q
term [[Int]]
sumOver)

  global :: q
global = [q] -> q
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [ Int -> q
wt Int
j | Int
j<-[Int
0..Int
n] ] q -> q -> q
forall a. Fractional a => a -> a -> a
/ (q
bq -> q -> q
forall a. Num a => a -> a -> a
-q
a)q -> Int -> q
forall a b. (Num a, Integral b) => a -> b -> a
^Int
d

  rkonst :: [Int] -> Integer
rkonst [Int]
ss = [Integer] -> Integer
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [ Int -> Integer
forall a. Integral a => a -> Integer
factorial Int
s Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Int -> Integer
forall a. Integral a => a -> Integer
factorial (Int
eInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
s) | (Int
s,Int
e) <- [Int] -> [Int] -> [(Int, Int)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int]
ss [Int]
es ]
  konst :: [Int] -> a
konst  [Int]
ss = Integer -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int -> Integer
forall a. Integral a => a -> Integer
paritySignValue ([Int] -> Int
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Int]
ss)) 
            a -> a -> a
forall a. Fractional a => a -> a -> a
/ Integer -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral ([Int] -> Integer
rkonst [Int]
ss)              
  denom :: [Int] -> q
denom  [Int]
ss = Int -> q
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n q -> q -> q
forall a. Num a => a -> a -> a
* q
a q -> q -> q
forall a. Num a => a -> a -> a
- Int -> q
forall a b. (Integral a, Num b) => a -> b
fromIntegral ([Int] -> Int
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [ Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int
s | ((Int
i,Int
e),Int
s) <- [(Int, Int)] -> [Int] -> [((Int, Int), Int)]
forall a b. [a] -> [b] -> [(a, b)]
zip [(Int, Int)]
ies [Int]
ss ]) q -> q -> q
forall a. Num a => a -> a -> a
* (q
aq -> q -> q
forall a. Num a => a -> a -> a
-q
b)
  term :: [Int] -> q
term   [Int]
ss = [Int] -> q
forall a. Fractional a => [Int] -> a
konst [Int]
ss q -> q -> q
forall a. Fractional a => a -> a -> a
/ [Int] -> q
denom [Int]
ss

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

-- | The dual class, recovered via from the localization formula via Lagrange
-- interpolation
localizeDual :: Partition -> ZMod AB
localizeDual :: Partition -> ZMod AB
localizeDual Partition
part = (X -> AB) -> FreeMod Integer X -> ZMod AB
forall a b c.
(Ord a, Ord b, Eq c, Num c) =>
(a -> b) -> FreeMod c a -> FreeMod c b
ZMod.mapBase X -> AB
forall (var :: Symbol). U var -> AB
worker (FreeMod Integer X -> ZMod AB) -> FreeMod Integer X -> ZMod AB
forall a b. (a -> b) -> a -> b
$ Partition -> FreeMod Integer X
localizeInterpolatedZ Partition
part where
  c :: Int
c = Partition -> Int
codim Partition
part
  worker :: U var -> AB
worker (U Int
i) = Int -> Int -> AB
AB (Int
cInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
i) Int
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 :: Partition -> QMod X
localizeInterpolatedQ part :: Partition
part@(Partition [Int]
xs) = QMod X
final where
  codim :: Int
codim = [Int] -> Int
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Int]
xs Int -> Int -> Int
forall a. Num a => a -> a -> a
- [Int] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Int]
xs
  bs :: [Rational]
bs = (Int -> Rational) -> [Int] -> [Rational]
forall a b. (a -> b) -> [a] -> [b]
map Int -> Rational
forall a b. (Integral a, Num b) => a -> b
fromIntegral [ Int
2..Int
codimInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
2 ]    :: [Rational]
  qs :: [Rational]
qs = [ Partition -> Rational -> Rational -> Rational
forall q. Fractional q => Partition -> q -> q -> q
localizeEval Partition
part Rational
1 Rational
b | Rational
b<-[Rational]
bs ] :: [Rational]
  final :: QMod X
final = Univariate Rational "x" -> QMod X
forall c (v :: Symbol). Univariate c v -> FreeMod c (U v)
unUni (Univariate Rational "x" -> QMod X)
-> Univariate Rational "x" -> QMod X
forall a b. (a -> b) -> a -> b
$ [(Rational, Rational)] -> Univariate Rational "x"
forall (var :: Symbol).
KnownSymbol var =>
[(Rational, Rational)] -> Univariate Rational var
lagrangeInterp ([Rational] -> [Rational] -> [(Rational, Rational)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Rational]
bs [Rational]
qs)

localizeInterpolatedZ :: Partition -> ZMod X
localizeInterpolatedZ :: Partition -> FreeMod Integer X
localizeInterpolatedZ = ((Rational -> Integer) -> QMod X -> FreeMod Integer X
forall b c2 c1.
(Ord b, Eq c2, Num c2) =>
(c1 -> c2) -> FreeMod c1 b -> FreeMod c2 b
ZMod.mapCoeff Rational -> Integer
f (QMod X -> FreeMod Integer X)
-> (Partition -> QMod X) -> Partition -> FreeMod Integer X
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Partition -> QMod X
localizeInterpolatedQ) where
  f :: Rational -> Integer
  f :: Rational -> Integer
f Rational
q = case Rational -> Integer
forall a. Ratio a -> a
denominator Rational
q of
          Integer
1 -> Rational -> Integer
forall a. Ratio a -> a
numerator Rational
q
          Integer
_ -> String -> Integer
forall a. HasCallStack => String -> a
error String
"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 ++ " ]"
-}