{-# LANGUAGE RebindableSyntax #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE FlexibleInstances #-}
{- |
Copyright   :  (c) Henning Thielemann 2004-2006

Maintainer  :  numericprelude@henning-thielemann.de
Stability   :  provisional
Portability :  requires multi-parameter type classes

Polynomials with negative and positive exponents.
-}
module MathObj.LaurentPolynomial where

import qualified MathObj.Polynomial  as Poly
import qualified MathObj.PowerSeries as PS
import qualified MathObj.PowerSeries.Core as PSCore

import qualified Algebra.VectorSpace  as VectorSpace
import qualified Algebra.Module       as Module
import qualified Algebra.Vector       as Vector
import qualified Algebra.Field        as Field
import qualified Algebra.Ring         as Ring
import qualified Algebra.Additive     as Additive
import qualified Algebra.ZeroTestable as ZeroTestable

import qualified Number.Complex as Complex

import qualified NumericPrelude.Numeric as NP

import NumericPrelude.Base    hiding (const, reverse, )
import NumericPrelude.Numeric hiding (div, negate, )

import qualified Data.List as List
import Data.List.HT (mapAdjacent)


{- | Polynomial including negative exponents -}

data T a = Cons {T a -> Int
expon :: Int, T a -> [a]
coeffs :: [a]}


{- * Basic Operations -}

const :: a -> T a
const :: a -> T a
const a
x = [a] -> T a
forall a. [a] -> T a
fromCoeffs [a
x]

(!) :: Additive.C a => T a -> Int -> a
(!) (Cons Int
xt [a]
x) Int
n =
   if Int
nInt -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<Int
xt
     then a
forall a. C a => a
zero
     else [a] -> a
forall a. [a] -> a
head (Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
drop (Int
nInt -> Int -> Int
forall a. C a => a -> a -> a
-Int
xt) [a]
x [a] -> [a] -> [a]
forall a. [a] -> [a] -> [a]
++ [a
forall a. C a => a
zero])

fromCoeffs :: [a] -> T a
fromCoeffs :: [a] -> T a
fromCoeffs = Int -> [a] -> T a
forall a. Int -> [a] -> T a
fromShiftCoeffs Int
0

fromShiftCoeffs :: Int -> [a] -> T a
fromShiftCoeffs :: Int -> [a] -> T a
fromShiftCoeffs = Int -> [a] -> T a
forall a. Int -> [a] -> T a
Cons

fromPolynomial :: Poly.T a -> T a
fromPolynomial :: T a -> T a
fromPolynomial = [a] -> T a
forall a. [a] -> T a
fromCoeffs ([a] -> T a) -> (T a -> [a]) -> T a -> T a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. T a -> [a]
forall a. T a -> [a]
Poly.coeffs

fromPowerSeries :: PS.T a -> T a
fromPowerSeries :: T a -> T a
fromPowerSeries = [a] -> T a
forall a. [a] -> T a
fromCoeffs ([a] -> T a) -> (T a -> [a]) -> T a -> T a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. T a -> [a]
forall a. T a -> [a]
PS.coeffs

bounds :: T a -> (Int, Int)
bounds :: T a -> (Int, Int)
bounds (Cons Int
xt [a]
x) = (Int
xt, Int
xt Int -> Int -> Int
forall a. C a => a -> a -> a
+ [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [a]
x Int -> Int -> Int
forall a. C a => a -> a -> a
- Int
1)

shift :: Int -> T a -> T a
shift :: Int -> T a -> T a
shift Int
t (Cons Int
xt [a]
x) = Int -> [a] -> T a
forall a. Int -> [a] -> T a
Cons (Int
xtInt -> Int -> Int
forall a. C a => a -> a -> a
+Int
t) [a]
x

{-# DEPRECATED translate "In order to avoid confusion with Polynomial.translate, use 'shift' instead" #-}
translate :: Int -> T a -> T a
translate :: Int -> T a -> T a
translate = Int -> T a -> T a
forall a. Int -> T a -> T a
shift


instance Functor T where
  fmap :: (a -> b) -> T a -> T b
fmap a -> b
f (Cons Int
xt [a]
xs) = Int -> [b] -> T b
forall a. Int -> [a] -> T a
Cons Int
xt ((a -> b) -> [a] -> [b]
forall a b. (a -> b) -> [a] -> [b]
map a -> b
f [a]
xs)


{- * Show -}

appPrec :: Int
appPrec :: Int
appPrec  = Int
10

instance (Show a) => Show (T a) where
  showsPrec :: Int -> T a -> ShowS
showsPrec Int
p (Cons Int
xt [a]
xs) =
    Bool -> ShowS -> ShowS
showParen (Int
p Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
appPrec)
       (String -> ShowS
showString String
"LaurentPolynomial.Cons " ShowS -> ShowS -> ShowS
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> ShowS
forall a. Show a => a -> ShowS
shows Int
xt ShowS -> ShowS -> ShowS
forall b c a. (b -> c) -> (a -> b) -> a -> c
.
        String -> ShowS
showString String
" " ShowS -> ShowS -> ShowS
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [a] -> ShowS
forall a. Show a => a -> ShowS
shows [a]
xs)

{- * Additive -}

add :: Additive.C a => T a -> T a -> T a
add :: T a -> T a -> T a
add (Cons Int
_ [])  T a
y          = T a
y
add  T a
x          (Cons Int
_ []) = T a
x
add (Cons Int
xt [a]
x) (Cons Int
yt [a]
y) =
   if Int
xt Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
yt
     then Int -> [a] -> T a
forall a. Int -> [a] -> T a
Cons Int
xt (Int -> [a] -> [a] -> [a]
forall a. C a => Int -> [a] -> [a] -> [a]
addShifted (Int
ytInt -> Int -> Int
forall a. C a => a -> a -> a
-Int
xt) [a]
x [a]
y)
     else Int -> [a] -> T a
forall a. Int -> [a] -> T a
Cons Int
yt (Int -> [a] -> [a] -> [a]
forall a. C a => Int -> [a] -> [a] -> [a]
addShifted (Int
xtInt -> Int -> Int
forall a. C a => a -> a -> a
-Int
yt) [a]
y [a]
x)

{-
Compute the value of a series of Laurent polynomials.

Requires that the start exponents constitute a (weakly) rising sequence,
where each exponent occurs only finitely often.

Cf. Synthesizer.Cut.arrange
-}
series :: (Additive.C a) => [T a] -> T a
series :: [T a] -> T a
series [] = [a] -> T a
forall a. [a] -> T a
fromCoeffs []
series [T a]
ps =
   let es :: [Int]
es = (T a -> Int) -> [T a] -> [Int]
forall a b. (a -> b) -> [a] -> [b]
map T a -> Int
forall a. T a -> Int
expon  [T a]
ps
       xs :: [[a]]
xs = (T a -> [a]) -> [T a] -> [[a]]
forall a b. (a -> b) -> [a] -> [b]
map T a -> [a]
forall a. T a -> [a]
coeffs [T a]
ps
       ds :: [Int]
ds = (Int -> Int -> Int) -> [Int] -> [Int]
forall a b. (a -> a -> b) -> [a] -> [b]
mapAdjacent Int -> Int -> Int
forall a. C a => a -> a -> a
subtract [Int]
es
   in  Int -> [a] -> T a
forall a. Int -> [a] -> T a
Cons ([Int] -> Int
forall a. [a] -> a
head [Int]
es) ([Int] -> [[a]] -> [a]
forall a. C a => [Int] -> [[a]] -> [a]
addShiftedMany [Int]
ds [[a]]
xs)

{- |
Add lists of numbers respecting a relative shift between the starts of the lists.
The shifts must be non-negative.
The list of relative shifts is one element shorter
than the list of summands.
Infinitely many summands are permitted,
provided that runs of zero shifts are all finite.


We could add the lists either with 'foldl' or with 'foldr',
'foldl' would be straightforward, but more time consuming (quadratic time)
whereas foldr is not so obvious but needs only linear time.

(stars denote the coefficients,
 frames denote what is contained in the interim results)
'foldl' sums this way:

> | | | *******************************
> | | +--------------------------------
> | |          ************************
> | +----------------------------------
> |                        ************
> +------------------------------------

I.e. 'foldl' would use much time find the time differences
by successive subtraction 1.

'foldr' mixes this way:

>     +--------------------------------
>     | *******************************
>     |      +-------------------------
>     |      | ************************
>     |      |           +-------------
>     |      |           | ************

-}
addShiftedMany :: (Additive.C a) => [Int] -> [[a]] -> [a]
addShiftedMany :: [Int] -> [[a]] -> [a]
addShiftedMany [Int]
ds [[a]]
xss =
   ((Int, [a]) -> [a] -> [a]) -> [a] -> [(Int, [a])] -> [a]
forall (t :: * -> *) a b.
Foldable t =>
(a -> b -> b) -> b -> t a -> b
foldr ((Int -> [a] -> [a] -> [a]) -> (Int, [a]) -> [a] -> [a]
forall a b c. (a -> b -> c) -> (a, b) -> c
uncurry Int -> [a] -> [a] -> [a]
forall a. C a => Int -> [a] -> [a] -> [a]
addShifted) [] ([Int] -> [[a]] -> [(Int, [a])]
forall a b. [a] -> [b] -> [(a, b)]
zip ([Int]
ds[Int] -> [Int] -> [Int]
forall a. [a] -> [a] -> [a]
++[Int
0]) [[a]]
xss)



addShifted :: Additive.C a => Int -> [a] -> [a] -> [a]
addShifted :: Int -> [a] -> [a] -> [a]
addShifted Int
del [a]
px [a]
py =
   let recurse :: Int -> [a] -> [a]
recurse Int
0 [a]
x      = [a] -> [a] -> [a]
forall a. C a => [a] -> [a] -> [a]
PSCore.add [a]
x [a]
py
       recurse Int
d []     = Int -> a -> [a]
forall a. Int -> a -> [a]
replicate Int
d a
forall a. C a => a
zero [a] -> [a] -> [a]
forall a. [a] -> [a] -> [a]
++ [a]
py
       recurse Int
d (a
x:[a]
xs) = a
x a -> [a] -> [a]
forall a. a -> [a] -> [a]
: Int -> [a] -> [a]
recurse (Int
dInt -> Int -> Int
forall a. C a => a -> a -> a
-Int
1) [a]
xs
   in  if Int
del Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
0
         then Int -> [a] -> [a]
recurse Int
del [a]
px
         else String -> [a]
forall a. HasCallStack => String -> a
error String
"LaurentPolynomial.addShifted: negative shift"


negate :: Additive.C a => T a -> T a
negate :: T a -> T a
negate (Cons Int
xt [a]
x) = Int -> [a] -> T a
forall a. Int -> [a] -> T a
Cons Int
xt ((a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map a -> a
forall a. C a => a -> a
NP.negate [a]
x)

sub :: Additive.C a => T a -> T a -> T a
sub :: T a -> T a -> T a
sub T a
x T a
y = T a -> T a -> T a
forall a. C a => T a -> T a -> T a
add T a
x (T a -> T a
forall a. C a => T a -> T a
negate T a
y)

instance Additive.C a => Additive.C (T a) where
   zero :: T a
zero   = [a] -> T a
forall a. [a] -> T a
fromCoeffs []
   + :: T a -> T a -> T a
(+)    = T a -> T a -> T a
forall a. C a => T a -> T a -> T a
add
   (-)    = T a -> T a -> T a
forall a. C a => T a -> T a -> T a
sub
   negate :: T a -> T a
negate = T a -> T a
forall a. C a => T a -> T a
negate


{- * Module -}

instance Vector.C T where
   zero :: T a
zero  = T a
forall a. C a => a
zero
   <+> :: T a -> T a -> T a
(<+>) = T a -> T a -> T a
forall a. C a => a -> a -> a
(+)
   *> :: a -> T a -> T a
(*>)  = a -> T a -> T a
forall (v :: * -> *) a. (Functor v, C a) => a -> v a -> v a
Vector.functorScale

instance (Module.C a b) => Module.C a (T b) where
    a
x *> :: a -> T b -> T b
*> Cons Int
yt [b]
ys = Int -> [b] -> T b
forall a. Int -> [a] -> T a
Cons Int
yt (a
x a -> [b] -> [b]
forall a v. C a v => a -> v -> v
*> [b]
ys)

instance (Field.C a, Module.C a b) => VectorSpace.C a (T b)


{- * Ring -}

mul :: Ring.C a => T a -> T a -> T a
mul :: T a -> T a -> T a
mul (Cons Int
xt [a]
x) (Cons Int
yt [a]
y) = Int -> [a] -> T a
forall a. Int -> [a] -> T a
Cons (Int
xtInt -> Int -> Int
forall a. C a => a -> a -> a
+Int
yt) ([a] -> [a] -> [a]
forall a. C a => [a] -> [a] -> [a]
PSCore.mul [a]
x [a]
y)

instance (Ring.C a) => Ring.C (T a) where
    one :: T a
one           = a -> T a
forall a. a -> T a
const a
forall a. C a => a
one
    fromInteger :: Integer -> T a
fromInteger Integer
n = a -> T a
forall a. a -> T a
const (Integer -> a
forall a. C a => Integer -> a
fromInteger Integer
n)
    * :: T a -> T a -> T a
(*)           = T a -> T a -> T a
forall a. C a => T a -> T a -> T a
mul


{- * Field.C -}

div :: (Field.C a, ZeroTestable.C a) => T a -> T a -> T a
div :: T a -> T a -> T a
div (Cons Int
xt [a]
xs) (Cons Int
yt [a]
ys) =
   let ([a]
xzero,[a]
x) = (a -> Bool) -> [a] -> ([a], [a])
forall a. (a -> Bool) -> [a] -> ([a], [a])
span a -> Bool
forall a. C a => a -> Bool
isZero [a]
xs
       ([a]
yzero,[a]
y) = (a -> Bool) -> [a] -> ([a], [a])
forall a. (a -> Bool) -> [a] -> ([a], [a])
span a -> Bool
forall a. C a => a -> Bool
isZero [a]
ys
   in  Int -> [a] -> T a
forall a. Int -> [a] -> T a
Cons (Int
xt Int -> Int -> Int
forall a. C a => a -> a -> a
- Int
yt Int -> Int -> Int
forall a. C a => a -> a -> a
+ [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [a]
xzero Int -> Int -> Int
forall a. C a => a -> a -> a
- [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [a]
yzero)
            ([a] -> [a] -> [a]
forall a. C a => [a] -> [a] -> [a]
PSCore.divide [a]
x [a]
y)

instance (Field.C a, ZeroTestable.C a) => Field.C (T a) where
   / :: T a -> T a -> T a
(/) = T a -> T a -> T a
forall a. (C a, C a) => T a -> T a -> T a
div

divExample :: T NP.Rational
divExample :: T Rational
divExample = T Rational -> T Rational -> T Rational
forall a. (C a, C a) => T a -> T a -> T a
div (Int -> [Rational] -> T Rational
forall a. Int -> [a] -> T a
Cons Int
1 [Rational
0,Rational
0,Rational
1,Rational
2,Rational
1]) (Int -> [Rational] -> T Rational
forall a. Int -> [a] -> T a
Cons Int
1 [Rational
0,Rational
0,Rational
0,Rational
1,Rational
1])




{- * Comparisons -}

{- |
Two polynomials may be stored differently.
This function checks whether two values of type @LaurentPolynomial@
actually represent the same polynomial.
-}
equivalent :: (Eq a, ZeroTestable.C a) => T a -> T a -> Bool
equivalent :: T a -> T a -> Bool
equivalent T a
xs T a
ys =
   let (Cons Int
xt [a]
x, Cons Int
yt [a]
y) =
          if T a -> Int
forall a. T a -> Int
expon T a
xs Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= T a -> Int
forall a. T a -> Int
expon T a
ys
            then (T a
xs,T a
ys)
            else (T a
ys,T a
xs)
       ([a]
xPref, [a]
xSuf) = Int -> [a] -> ([a], [a])
forall a. Int -> [a] -> ([a], [a])
splitAt (Int
ytInt -> Int -> Int
forall a. C a => a -> a -> a
-Int
xt) [a]
x
       aux :: [a] -> [a] -> Bool
aux (a
a:[a]
as) (a
b:[a]
bs) = a
a a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
b Bool -> Bool -> Bool
&& [a] -> [a] -> Bool
aux [a]
as [a]
bs
       aux []     [a]
bs     = (a -> Bool) -> [a] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all a -> Bool
forall a. C a => a -> Bool
isZero [a]
bs
       aux [a]
as     []     = (a -> Bool) -> [a] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all a -> Bool
forall a. C a => a -> Bool
isZero [a]
as
   in  (a -> Bool) -> [a] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all a -> Bool
forall a. C a => a -> Bool
isZero [a]
xPref  Bool -> Bool -> Bool
&&  [a] -> [a] -> Bool
forall a. (Eq a, C a) => [a] -> [a] -> Bool
aux [a]
xSuf [a]
y

instance (Eq a, ZeroTestable.C a) => Eq (T a) where
   == :: T a -> T a -> Bool
(==) = T a -> T a -> Bool
forall a. (Eq a, C a) => T a -> T a -> Bool
equivalent


identical :: (Eq a) => T a -> T a -> Bool
identical :: T a -> T a -> Bool
identical (Cons Int
xt [a]
xs) (Cons Int
yt [a]
ys) =
   Int
xtInt -> Int -> Bool
forall a. Eq a => a -> a -> Bool
==Int
yt Bool -> Bool -> Bool
&& [a]
xs [a] -> [a] -> Bool
forall a. Eq a => a -> a -> Bool
== [a]
ys


{- |
Check whether a Laurent polynomial has only the absolute term,
that is, it represents the constant polynomial.
-}
isAbsolute :: (ZeroTestable.C a) => T a -> Bool
isAbsolute :: T a -> Bool
isAbsolute (Cons Int
xt [a]
x) =
   [Bool] -> Bool
forall (t :: * -> *). Foldable t => t Bool -> Bool
and ((Int -> a -> Bool) -> [Int] -> [a] -> [Bool]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\Int
i a
xi -> a -> Bool
forall a. C a => a -> Bool
isZero a
xi Bool -> Bool -> Bool
|| Int
iInt -> Int -> Bool
forall a. Eq a => a -> a -> Bool
==Int
0) [Int
xt..] [a]
x)



{- * Transformations of arguments -}

{- | p(z) -> p(-z) -}
alternate :: Additive.C a => T a -> T a
alternate :: T a -> T a
alternate (Cons Int
xt [a]
x) =
   Int -> [a] -> T a
forall a. Int -> [a] -> T a
Cons Int
xt (((a -> a) -> a -> a) -> [a -> a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (a -> a) -> a -> a
forall a. a -> a
id (Int -> [a -> a] -> [a -> a]
forall a. Int -> [a] -> [a]
drop (Int -> Int -> Int
forall a. C a => a -> a -> a
mod Int
xt Int
2) ([a -> a] -> [a -> a]
forall a. [a] -> [a]
cycle [a -> a
forall a. a -> a
id,a -> a
forall a. C a => a -> a
NP.negate])) [a]
x)

{- | p(z) -> p(1\/z) -}
reverse :: T a -> T a
reverse :: T a -> T a
reverse (Cons Int
xt [a]
x) =
   Int -> [a] -> T a
forall a. Int -> [a] -> T a
Cons (Int
1 Int -> Int -> Int
forall a. C a => a -> a -> a
- Int
xt Int -> Int -> Int
forall a. C a => a -> a -> a
- [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [a]
x) ([a] -> [a]
forall a. [a] -> [a]
List.reverse [a]
x)

{- |
p(exp(i·x)) -> conjugate(p(exp(i·x)))

If you interpret @(p*)@ as a linear operator on the space of Laurent polynomials,
then @(adjoint p *)@ is the adjoint operator.
-}
adjoint :: Additive.C a => T (Complex.T a) -> T (Complex.T a)
adjoint :: T (T a) -> T (T a)
adjoint T (T a)
x =
   let (Cons Int
yt [T a]
y) = T (T a) -> T (T a)
forall a. T a -> T a
reverse T (T a)
x
   in  (Int -> [T a] -> T (T a)
forall a. Int -> [a] -> T a
Cons Int
yt ((T a -> T a) -> [T a] -> [T a]
forall a b. (a -> b) -> [a] -> [b]
map T a -> T a
forall a. C a => T a -> T a
Complex.conjugate [T a]
y))