{-# 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
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)
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
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
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"