{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE GeneralizedNewtypeDeriving #-}
{-# LANGUAGE PatternSynonyms #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE ViewPatterns #-}
module Data.Poly.Laurent
( Laurent
, VLaurent
, ULaurent
, unLaurent
, toLaurent
, leading
, monomial
, scale
, pattern X
, (^-)
, eval
, subst
, deriv
, LaurentOverField(..)
) where
import Prelude hiding (quotRem, quot, rem, gcd)
import Control.Arrow (first)
import Control.DeepSeq (NFData(..))
import Data.Euclidean (GcdDomain(..), Euclidean(..), Field)
import Data.List (intersperse)
import Data.Semiring (Semiring(..), Ring())
import qualified Data.Semiring as Semiring
import qualified Data.Vector as V
import qualified Data.Vector.Generic as G
import qualified Data.Vector.Unboxed as U
import Data.Poly.Internal.Dense (Poly(..))
import qualified Data.Poly.Internal.Dense as Dense
import Data.Poly.Internal.Dense.Field ()
import Data.Poly.Internal.Dense.GcdDomain ()
import Data.Poly.Internal.PolyOverField
data Laurent v a = Laurent !Int !(Poly v a)
deriving (Eq, Ord)
unLaurent :: Laurent v a -> (Int, Poly v a)
unLaurent (Laurent off poly) = (off, poly)
toLaurent
:: (Eq a, Semiring a, G.Vector v a)
=> Int
-> Poly v a
-> Laurent v a
toLaurent off (Poly xs) = go 0
where
go k
| k >= G.length xs
= Laurent 0 zero
| G.unsafeIndex xs k == zero
= go (k + 1)
| otherwise
= Laurent (off + k) (Poly (G.unsafeDrop k xs))
{-# INLINE toLaurent #-}
toLaurentNum
:: (Eq a, Num a, G.Vector v a)
=> Int
-> Poly v a
-> Laurent v a
toLaurentNum off (Poly xs) = go 0
where
go k
| k >= G.length xs
= Laurent 0 0
| G.unsafeIndex xs k == 0
= go (k + 1)
| otherwise
= Laurent (off + k) (Poly (G.unsafeDrop k xs))
{-# INLINE toLaurentNum #-}
instance NFData (v a) => NFData (Laurent v a) where
rnf (Laurent off poly) = rnf off `seq` rnf poly
instance (Show a, G.Vector v a) => Show (Laurent v a) where
showsPrec d (Laurent off poly)
| G.null (unPoly poly)
= showString "0"
| otherwise
= showParen (d > 0)
$ foldl (.) id
$ intersperse (showString " + ")
$ G.ifoldl (\acc i c -> showCoeff (i + off) c : acc) []
$ unPoly poly
where
showCoeff 0 c = showsPrec 7 c
showCoeff 1 c = showsPrec 7 c . showString " * X"
showCoeff i c = showsPrec 7 c . showString (" * X^" ++ show i)
type VLaurent = Laurent V.Vector
type ULaurent = Laurent U.Vector
leading :: G.Vector v a => Laurent v a -> Maybe (Int, a)
leading (Laurent off poly) = first ((+ off) . fromIntegral) <$> Dense.leading poly
instance (Eq a, Num a, G.Vector v a) => Num (Laurent v a) where
Laurent off1 poly1 * Laurent off2 poly2 = toLaurentNum (off1 + off2) (poly1 * poly2)
Laurent off1 poly1 + Laurent off2 poly2 = case off1 `compare` off2 of
LT -> toLaurentNum off1 (poly1 + Dense.scale (fromIntegral $ off2 - off1) 1 poly2)
EQ -> toLaurentNum off1 (poly1 + poly2)
GT -> toLaurentNum off2 (Dense.scale (fromIntegral $ off1 - off2) 1 poly1 + poly2)
Laurent off1 poly1 - Laurent off2 poly2 = case off1 `compare` off2 of
LT -> toLaurentNum off1 (poly1 - Dense.scale (fromIntegral $ off2 - off1) 1 poly2)
EQ -> toLaurentNum off1 (poly1 - poly2)
GT -> toLaurentNum off2 (Dense.scale (fromIntegral $ off1 - off2) 1 poly1 - poly2)
negate (Laurent off poly) = Laurent off (negate poly)
abs = id
signum = const 1
fromInteger n = Laurent 0 (fromInteger n)
{-# INLINE (+) #-}
{-# INLINE (-) #-}
{-# INLINE negate #-}
{-# INLINE fromInteger #-}
{-# INLINE (*) #-}
instance (Eq a, Semiring a, G.Vector v a) => Semiring (Laurent v a) where
zero = Laurent 0 zero
one = Laurent 0 one
Laurent off1 poly1 `times` Laurent off2 poly2 =
toLaurent (off1 + off2) (poly1 `times` poly2)
Laurent off1 poly1 `plus` Laurent off2 poly2 = case off1 `compare` off2 of
LT -> toLaurent off1 (poly1 `plus` Dense.scale' (fromIntegral $ off2 - off1) one poly2)
EQ -> toLaurent off1 (poly1 `plus` poly2)
GT -> toLaurent off2 (Dense.scale' (fromIntegral $ off1 - off2) one poly1 `plus` poly2)
fromNatural n = Laurent 0 (fromNatural n)
{-# INLINE zero #-}
{-# INLINE one #-}
{-# INLINE plus #-}
{-# INLINE times #-}
{-# INLINE fromNatural #-}
instance (Eq a, Ring a, G.Vector v a) => Ring (Laurent v a) where
negate (Laurent off poly) = Laurent off (Semiring.negate poly)
monomial :: (Eq a, Semiring a, G.Vector v a) => Int -> a -> Laurent v a
monomial p c
| c == zero = Laurent 0 zero
| otherwise = Laurent p (Dense.monomial' 0 c)
{-# INLINE monomial #-}
scale :: (Eq a, Semiring a, G.Vector v a) => Int -> a -> Laurent v a -> Laurent v a
scale yp yc (Laurent off poly) = toLaurent (off + yp) (Dense.scale' 0 yc poly)
eval :: (Field a, G.Vector v a) => Laurent v a -> a -> a
eval (Laurent off poly) x = Dense.eval' poly x `times`
(if off >= 0 then x Semiring.^ off else quot one x Semiring.^ (- off))
{-# INLINE eval #-}
subst :: (Eq a, Semiring a, G.Vector v a, G.Vector w a) => Poly v a -> Laurent w a -> Laurent w a
subst = Dense.substitute' (scale 0)
{-# INLINE subst #-}
deriv :: (Eq a, Ring a, G.Vector v a) => Laurent v a -> Laurent v a
deriv (Laurent off (Poly xs)) =
toLaurent (off - 1) $ Dense.toPoly' $ G.imap (times . Semiring.fromIntegral . (+ off)) xs
{-# INLINE deriv #-}
pattern X :: (Eq a, Semiring a, G.Vector v a, Eq (v a)) => Laurent v a
pattern X <- ((==) var -> True)
where X = var
var :: forall a v. (Eq a, Semiring a, G.Vector v a, Eq (v a)) => Laurent v a
var
| (one :: a) == zero = Laurent 0 zero
| otherwise = Laurent 1 one
{-# INLINE var #-}
(^-)
:: (Eq a, Semiring a, G.Vector v a, Eq (v a))
=> Laurent v a
-> Int
-> Laurent v a
X^-n = monomial (negate n) one
_^-_ = error "(^-) can be applied only to X"
instance (Eq a, Ring a, GcdDomain a, Eq (v a), G.Vector v a) => GcdDomain (Laurent v a) where
divide (Laurent off1 poly1) (Laurent off2 poly2) =
toLaurent (off1 - off2) <$> divide poly1 poly2
{-# INLINE divide #-}
gcd (Laurent _ poly1) (Laurent _ poly2) =
toLaurent 0 (gcd poly1 poly2)
{-# INLINE gcd #-}
newtype LaurentOverField laurent = LaurentOverField { unLaurentOverField :: laurent }
deriving (Eq, NFData, Num, Ord, Ring, Semiring, Show)
instance (Eq a, Eq (v a), Field a, G.Vector v a) => GcdDomain (LaurentOverField (Laurent v a)) where
divide (LaurentOverField (Laurent off1 poly1)) (LaurentOverField (Laurent off2 poly2)) =
LaurentOverField . toLaurent (off1 - off2) . unPolyOverField <$> divide (PolyOverField poly1) (PolyOverField poly2)
gcd (LaurentOverField (Laurent _ poly1)) (LaurentOverField (Laurent _ poly2)) =
LaurentOverField (toLaurent 0 (unPolyOverField (gcd (PolyOverField poly1) (PolyOverField poly2))))
{-# INLINE gcd #-}