-- |
-- Module:      Data.Poly.Internal.Dense.Laurent
-- Copyright:   (c) 2020 Andrew Lelechenko
-- Licence:     BSD3
-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>
--
-- <https://en.wikipedia.org/wiki/Laurent_polynomial Laurent polynomials>.
--

{-# LANGUAGE FlexibleInstances          #-}
{-# LANGUAGE KindSignatures             #-}
{-# LANGUAGE PatternSynonyms            #-}
{-# LANGUAGE ScopedTypeVariables        #-}
{-# LANGUAGE ViewPatterns               #-}

module Data.Poly.Internal.Dense.Laurent
  ( Laurent
  , VLaurent
  , ULaurent
  , unLaurent
  , toLaurent
  , leading
  , monomial
  , scale
  , pattern X
  , (^-)
  , eval
  , subst
  , deriv
  ) where

import Prelude hiding (quotRem, quot, rem, gcd, lcm)
import Control.Arrow (first)
import Control.DeepSeq (NFData(..))
import Control.Exception
import Data.Euclidean (GcdDomain(..), Euclidean(..), Field)
import Data.Kind
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 ()

-- | <https://en.wikipedia.org/wiki/Laurent_polynomial Laurent polynomials>
-- of one variable with coefficients from @a@,
-- backed by a 'G.Vector' @v@ (boxed, unboxed, storable, etc.).
--
-- Use the pattern 'X' and the '^-' operator for construction:
--
-- >>> (X + 1) + (X^-1 - 1) :: VLaurent Integer
-- 1 * X + 0 + 1 * X^-1
-- >>> (X + 1) * (1 - X^-1) :: ULaurent Int
-- 1 * X + 0 + (-1) * X^-1
--
-- Polynomials are stored normalized, without leading
-- and trailing
-- zero coefficients, so 0 * X + 1 + 0 * X^-1 equals to 1.
--
-- The 'Ord' instance does not make much sense mathematically,
-- it is defined only for the sake of 'Data.Set.Set', 'Data.Map.Map', etc.
--
-- Due to being polymorphic by multiple axis, the performance of `Laurent` crucially
-- depends on specialisation of instances. Clients are strongly recommended
-- to compile with @ghc-options:@ @-fspecialise-aggressively@ and suggested to enable @-O2@.
--
-- @since 0.4.0.0
--
data Laurent (v :: Type -> Type) (a :: Type) = Laurent !Int !(Poly v a)
  deriving (Laurent v a -> Laurent v a -> Bool
(Laurent v a -> Laurent v a -> Bool)
-> (Laurent v a -> Laurent v a -> Bool) -> Eq (Laurent v a)
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
forall (v :: * -> *) a.
Eq (v a) =>
Laurent v a -> Laurent v a -> Bool
$c== :: forall (v :: * -> *) a.
Eq (v a) =>
Laurent v a -> Laurent v a -> Bool
== :: Laurent v a -> Laurent v a -> Bool
$c/= :: forall (v :: * -> *) a.
Eq (v a) =>
Laurent v a -> Laurent v a -> Bool
/= :: Laurent v a -> Laurent v a -> Bool
Eq, Eq (Laurent v a)
Eq (Laurent v a) =>
(Laurent v a -> Laurent v a -> Ordering)
-> (Laurent v a -> Laurent v a -> Bool)
-> (Laurent v a -> Laurent v a -> Bool)
-> (Laurent v a -> Laurent v a -> Bool)
-> (Laurent v a -> Laurent v a -> Bool)
-> (Laurent v a -> Laurent v a -> Laurent v a)
-> (Laurent v a -> Laurent v a -> Laurent v a)
-> Ord (Laurent v a)
Laurent v a -> Laurent v a -> Bool
Laurent v a -> Laurent v a -> Ordering
Laurent v a -> Laurent v a -> Laurent v a
forall a.
Eq a =>
(a -> a -> Ordering)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> a)
-> (a -> a -> a)
-> Ord a
forall (v :: * -> *) a. Ord (v a) => Eq (Laurent v a)
forall (v :: * -> *) a.
Ord (v a) =>
Laurent v a -> Laurent v a -> Bool
forall (v :: * -> *) a.
Ord (v a) =>
Laurent v a -> Laurent v a -> Ordering
forall (v :: * -> *) a.
Ord (v a) =>
Laurent v a -> Laurent v a -> Laurent v a
$ccompare :: forall (v :: * -> *) a.
Ord (v a) =>
Laurent v a -> Laurent v a -> Ordering
compare :: Laurent v a -> Laurent v a -> Ordering
$c< :: forall (v :: * -> *) a.
Ord (v a) =>
Laurent v a -> Laurent v a -> Bool
< :: Laurent v a -> Laurent v a -> Bool
$c<= :: forall (v :: * -> *) a.
Ord (v a) =>
Laurent v a -> Laurent v a -> Bool
<= :: Laurent v a -> Laurent v a -> Bool
$c> :: forall (v :: * -> *) a.
Ord (v a) =>
Laurent v a -> Laurent v a -> Bool
> :: Laurent v a -> Laurent v a -> Bool
$c>= :: forall (v :: * -> *) a.
Ord (v a) =>
Laurent v a -> Laurent v a -> Bool
>= :: Laurent v a -> Laurent v a -> Bool
$cmax :: forall (v :: * -> *) a.
Ord (v a) =>
Laurent v a -> Laurent v a -> Laurent v a
max :: Laurent v a -> Laurent v a -> Laurent v a
$cmin :: forall (v :: * -> *) a.
Ord (v a) =>
Laurent v a -> Laurent v a -> Laurent v a
min :: Laurent v a -> Laurent v a -> Laurent v a
Ord)

-- | Deconstruct a 'Laurent' polynomial into an offset (largest possible)
-- and a regular polynomial.
--
-- >>> unLaurent (2 * X + 1 :: ULaurent Int)
-- (0,2 * X + 1)
-- >>> unLaurent (1 + 2 * X^-1 :: ULaurent Int)
-- (-1,1 * X + 2)
-- >>> unLaurent (2 * X^2 + X :: ULaurent Int)
-- (1,2 * X + 1)
-- >>> unLaurent (0 :: ULaurent Int)
-- (0,0)
--
-- @since 0.4.0.0
unLaurent :: Laurent v a -> (Int, Poly v a)
unLaurent :: forall (v :: * -> *) a. Laurent v a -> (Int, Poly v a)
unLaurent (Laurent Int
off Poly v a
poly) = (Int
off, Poly v a
poly)

-- | Construct 'Laurent' polynomial from an offset and a regular polynomial.
-- One can imagine it as 'Data.Poly.Semiring.scale', but allowing negative offsets.
--
-- >>> toLaurent 2 (2 * Data.Poly.X + 1) :: ULaurent Int
-- 2 * X^3 + 1 * X^2
-- >>> toLaurent (-2) (2 * Data.Poly.X + 1) :: ULaurent Int
-- 2 * X^-1 + 1 * X^-2
--
-- @since 0.4.0.0
toLaurent
  :: (Eq a, Semiring a, G.Vector v a)
  => Int
  -> Poly v a
  -> Laurent v a
toLaurent :: forall a (v :: * -> *).
(Eq a, Semiring a, Vector v a) =>
Int -> Poly v a -> Laurent v a
toLaurent Int
off (Poly v a
xs) = Int -> Laurent v a
go Int
0
  where
    go :: Int -> Laurent v a
go Int
k
      | Int
k Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= v a -> Int
forall (v :: * -> *) a. Vector v a => v a -> Int
G.length v a
xs
      = Int -> Poly v a -> Laurent v a
forall (v :: * -> *) a. Int -> Poly v a -> Laurent v a
Laurent Int
0 Poly v a
forall a. Semiring a => a
zero
      | v a -> Int -> a
forall (v :: * -> *) a. Vector v a => v a -> Int -> a
G.unsafeIndex v a
xs Int
k a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
forall a. Semiring a => a
zero
      = Int -> Laurent v a
go (Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1)
      | Bool
otherwise
      = Int -> Poly v a -> Laurent v a
forall (v :: * -> *) a. Int -> Poly v a -> Laurent v a
Laurent (Int
off Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
k) (v a -> Poly v a
forall (v :: * -> *) a. v a -> Poly v a
Poly (Int -> v a -> v a
forall (v :: * -> *) a. Vector v a => Int -> v a -> v a
G.unsafeDrop Int
k v a
xs))
{-# INLINE toLaurent #-}

toLaurentNum
  :: (Eq a, Num a, G.Vector v a)
  => Int
  -> Poly v a
  -> Laurent v a
toLaurentNum :: forall a (v :: * -> *).
(Eq a, Num a, Vector v a) =>
Int -> Poly v a -> Laurent v a
toLaurentNum Int
off (Poly v a
xs) = Int -> Laurent v a
go Int
0
  where
    go :: Int -> Laurent v a
go Int
k
      | Int
k Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= v a -> Int
forall (v :: * -> *) a. Vector v a => v a -> Int
G.length v a
xs
      = Int -> Poly v a -> Laurent v a
forall (v :: * -> *) a. Int -> Poly v a -> Laurent v a
Laurent Int
0 Poly v a
0
      | v a -> Int -> a
forall (v :: * -> *) a. Vector v a => v a -> Int -> a
G.unsafeIndex v a
xs Int
k a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0
      = Int -> Laurent v a
go (Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1)
      | Bool
otherwise
      = Int -> Poly v a -> Laurent v a
forall (v :: * -> *) a. Int -> Poly v a -> Laurent v a
Laurent (Int
off Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
k) (v a -> Poly v a
forall (v :: * -> *) a. v a -> Poly v a
Poly (Int -> v a -> v a
forall (v :: * -> *) a. Vector v a => Int -> v a -> v a
G.unsafeDrop Int
k v a
xs))
{-# INLINE toLaurentNum #-}

instance NFData (v a) => NFData (Laurent v a) where
  rnf :: Laurent v a -> ()
rnf (Laurent Int
off Poly v a
poly) = Int -> ()
forall a. NFData a => a -> ()
rnf Int
off () -> () -> ()
forall a b. a -> b -> b
`seq` Poly v a -> ()
forall a. NFData a => a -> ()
rnf Poly v a
poly

instance (Show a, G.Vector v a) => Show (Laurent v a) where
  showsPrec :: Int -> Laurent v a -> ShowS
showsPrec Int
d (Laurent Int
off Poly v a
poly)
    | v a -> Bool
forall (v :: * -> *) a. Vector v a => v a -> Bool
G.null (Poly v a -> v a
forall (v :: * -> *) a. Poly v a -> v a
unPoly Poly v a
poly)
      = String -> ShowS
showString String
"0"
    | Bool
otherwise
      = Bool -> ShowS -> ShowS
showParen (Int
d Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
0)
      (ShowS -> ShowS) -> ShowS -> ShowS
forall a b. (a -> b) -> a -> b
$ (ShowS -> ShowS -> ShowS) -> ShowS -> [ShowS] -> ShowS
forall b a. (b -> a -> b) -> b -> [a] -> b
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl ShowS -> ShowS -> ShowS
forall b c a. (b -> c) -> (a -> b) -> a -> c
(.) ShowS
forall a. a -> a
id
      ([ShowS] -> ShowS) -> [ShowS] -> ShowS
forall a b. (a -> b) -> a -> b
$ ShowS -> [ShowS] -> [ShowS]
forall a. a -> [a] -> [a]
intersperse (String -> ShowS
showString String
" + ")
      ([ShowS] -> [ShowS]) -> [ShowS] -> [ShowS]
forall a b. (a -> b) -> a -> b
$ ([ShowS] -> Int -> a -> [ShowS]) -> [ShowS] -> v a -> [ShowS]
forall (v :: * -> *) b a.
Vector v b =>
(a -> Int -> b -> a) -> a -> v b -> a
G.ifoldl (\[ShowS]
acc Int
i a
c -> Int -> a -> ShowS
forall {a} {a}. (Eq a, Num a, Show a, Show a) => a -> a -> ShowS
showCoeff (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
off) a
c ShowS -> [ShowS] -> [ShowS]
forall a. a -> [a] -> [a]
: [ShowS]
acc) []
      (v a -> [ShowS]) -> v a -> [ShowS]
forall a b. (a -> b) -> a -> b
$ Poly v a -> v a
forall (v :: * -> *) a. Poly v a -> v a
unPoly Poly v a
poly
    where
      -- Negative powers should be displayed without surrounding brackets
      showCoeff :: a -> a -> ShowS
showCoeff a
0 a
c = Int -> a -> ShowS
forall a. Show a => Int -> a -> ShowS
showsPrec Int
7 a
c
      showCoeff a
1 a
c = Int -> a -> ShowS
forall a. Show a => Int -> a -> ShowS
showsPrec Int
7 a
c ShowS -> ShowS -> ShowS
forall b c a. (b -> c) -> (a -> b) -> a -> c
. String -> ShowS
showString String
" * X"
      showCoeff a
i a
c = Int -> a -> ShowS
forall a. Show a => Int -> a -> ShowS
showsPrec Int
7 a
c ShowS -> ShowS -> ShowS
forall b c a. (b -> c) -> (a -> b) -> a -> c
. String -> ShowS
showString (String
" * X^" String -> ShowS
forall a. [a] -> [a] -> [a]
++ a -> String
forall a. Show a => a -> String
show a
i)

-- | Laurent polynomials backed by boxed vectors.
--
-- @since 0.4.0.0
type VLaurent = Laurent V.Vector

-- | Laurent polynomials backed by unboxed vectors.
--
-- @since 0.4.0.0
type ULaurent = Laurent U.Vector

-- | Return the leading power and coefficient of a non-zero polynomial.
--
-- >>> leading ((2 * X + 1) * (2 * X^2 - 1) :: ULaurent Int)
-- Just (3,4)
-- >>> leading (0 :: ULaurent Int)
-- Nothing
--
-- @since 0.4.0.0
leading :: G.Vector v a => Laurent v a -> Maybe (Int, a)
leading :: forall (v :: * -> *) a. Vector v a => Laurent v a -> Maybe (Int, a)
leading (Laurent Int
off Poly v a
poly) = (Word -> Int) -> (Word, a) -> (Int, a)
forall b c d. (b -> c) -> (b, d) -> (c, d)
forall (a :: * -> * -> *) b c d.
Arrow a =>
a b c -> a (b, d) (c, d)
first ((Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
off) (Int -> Int) -> (Word -> Int) -> Word -> Int
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Word -> Int
forall a b. (Integral a, Num b) => a -> b
fromIntegral) ((Word, a) -> (Int, a)) -> Maybe (Word, a) -> Maybe (Int, a)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Poly v a -> Maybe (Word, a)
forall (v :: * -> *) a. Vector v a => Poly v a -> Maybe (Word, a)
Dense.leading Poly v a
poly
{-# INLINABLE leading #-}

-- | Note that 'abs' = 'id' and 'signum' = 'const' 1.
instance (Eq a, Num a, G.Vector v a) => Num (Laurent v a) where
  Laurent Int
off1 Poly v a
poly1 * :: Laurent v a -> Laurent v a -> Laurent v a
* Laurent Int
off2 Poly v a
poly2 = Int -> Poly v a -> Laurent v a
forall a (v :: * -> *).
(Eq a, Num a, Vector v a) =>
Int -> Poly v a -> Laurent v a
toLaurentNum (Int
off1 Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
off2) (Poly v a
poly1 Poly v a -> Poly v a -> Poly v a
forall a. Num a => a -> a -> a
* Poly v a
poly2)
  Laurent Int
off1 Poly v a
poly1 + :: Laurent v a -> Laurent v a -> Laurent v a
+ Laurent Int
off2 Poly v a
poly2 = case Int
off1 Int -> Int -> Ordering
forall a. Ord a => a -> a -> Ordering
`compare` Int
off2 of
    Ordering
LT -> Int -> Poly v a -> Laurent v a
forall a (v :: * -> *).
(Eq a, Num a, Vector v a) =>
Int -> Poly v a -> Laurent v a
toLaurentNum Int
off1 (Poly v a
poly1 Poly v a -> Poly v a -> Poly v a
forall a. Num a => a -> a -> a
+ Word -> a -> Poly v a -> Poly v a
forall a (v :: * -> *).
(Eq a, Num a, Vector v a) =>
Word -> a -> Poly v a -> Poly v a
Dense.scale (Int -> Word
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int -> Word) -> Int -> Word
forall a b. (a -> b) -> a -> b
$ Int
off2 Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
off1) a
1 Poly v a
poly2)
    Ordering
EQ -> Int -> Poly v a -> Laurent v a
forall a (v :: * -> *).
(Eq a, Num a, Vector v a) =>
Int -> Poly v a -> Laurent v a
toLaurentNum Int
off1 (Poly v a
poly1 Poly v a -> Poly v a -> Poly v a
forall a. Num a => a -> a -> a
+ Poly v a
poly2)
    Ordering
GT -> Int -> Poly v a -> Laurent v a
forall a (v :: * -> *).
(Eq a, Num a, Vector v a) =>
Int -> Poly v a -> Laurent v a
toLaurentNum Int
off2 (Word -> a -> Poly v a -> Poly v a
forall a (v :: * -> *).
(Eq a, Num a, Vector v a) =>
Word -> a -> Poly v a -> Poly v a
Dense.scale (Int -> Word
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int -> Word) -> Int -> Word
forall a b. (a -> b) -> a -> b
$ Int
off1 Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
off2) a
1 Poly v a
poly1 Poly v a -> Poly v a -> Poly v a
forall a. Num a => a -> a -> a
+ Poly v a
poly2)
  Laurent Int
off1 Poly v a
poly1 - :: Laurent v a -> Laurent v a -> Laurent v a
- Laurent Int
off2 Poly v a
poly2 = case Int
off1 Int -> Int -> Ordering
forall a. Ord a => a -> a -> Ordering
`compare` Int
off2 of
    Ordering
LT -> Int -> Poly v a -> Laurent v a
forall a (v :: * -> *).
(Eq a, Num a, Vector v a) =>
Int -> Poly v a -> Laurent v a
toLaurentNum Int
off1 (Poly v a
poly1 Poly v a -> Poly v a -> Poly v a
forall a. Num a => a -> a -> a
- Word -> a -> Poly v a -> Poly v a
forall a (v :: * -> *).
(Eq a, Num a, Vector v a) =>
Word -> a -> Poly v a -> Poly v a
Dense.scale (Int -> Word
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int -> Word) -> Int -> Word
forall a b. (a -> b) -> a -> b
$ Int
off2 Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
off1) a
1 Poly v a
poly2)
    Ordering
EQ -> Int -> Poly v a -> Laurent v a
forall a (v :: * -> *).
(Eq a, Num a, Vector v a) =>
Int -> Poly v a -> Laurent v a
toLaurentNum Int
off1 (Poly v a
poly1 Poly v a -> Poly v a -> Poly v a
forall a. Num a => a -> a -> a
- Poly v a
poly2)
    Ordering
GT -> Int -> Poly v a -> Laurent v a
forall a (v :: * -> *).
(Eq a, Num a, Vector v a) =>
Int -> Poly v a -> Laurent v a
toLaurentNum Int
off2 (Word -> a -> Poly v a -> Poly v a
forall a (v :: * -> *).
(Eq a, Num a, Vector v a) =>
Word -> a -> Poly v a -> Poly v a
Dense.scale (Int -> Word
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int -> Word) -> Int -> Word
forall a b. (a -> b) -> a -> b
$ Int
off1 Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
off2) a
1 Poly v a
poly1 Poly v a -> Poly v a -> Poly v a
forall a. Num a => a -> a -> a
- Poly v a
poly2)
  negate :: Laurent v a -> Laurent v a
negate (Laurent Int
off Poly v a
poly) = Int -> Poly v a -> Laurent v a
forall (v :: * -> *) a. Int -> Poly v a -> Laurent v a
Laurent Int
off (Poly v a -> Poly v a
forall a. Num a => a -> a
negate Poly v a
poly)
  abs :: Laurent v a -> Laurent v a
abs = Laurent v a -> Laurent v a
forall a. a -> a
id
  signum :: Laurent v a -> Laurent v a
signum = Laurent v a -> Laurent v a -> Laurent v a
forall a b. a -> b -> a
const Laurent v a
1
  fromInteger :: Integer -> Laurent v a
fromInteger Integer
n = Int -> Poly v a -> Laurent v a
forall (v :: * -> *) a. Int -> Poly v a -> Laurent v a
Laurent Int
0 (Integer -> Poly v a
forall a. Num a => Integer -> a
fromInteger Integer
n)
  {-# INLINE (+) #-}
  {-# INLINE (-) #-}
  {-# INLINE negate #-}
  {-# INLINE fromInteger #-}
  {-# INLINE (*) #-}

instance (Eq a, Semiring a, G.Vector v a) => Semiring (Laurent v a) where
  zero :: Laurent v a
zero = Int -> Poly v a -> Laurent v a
forall (v :: * -> *) a. Int -> Poly v a -> Laurent v a
Laurent Int
0 Poly v a
forall a. Semiring a => a
zero
  one :: Laurent v a
one  = Int -> Poly v a -> Laurent v a
forall (v :: * -> *) a. Int -> Poly v a -> Laurent v a
Laurent Int
0 Poly v a
forall a. Semiring a => a
one
  Laurent Int
off1 Poly v a
poly1 times :: Laurent v a -> Laurent v a -> Laurent v a
`times` Laurent Int
off2 Poly v a
poly2 =
    Int -> Poly v a -> Laurent v a
forall a (v :: * -> *).
(Eq a, Semiring a, Vector v a) =>
Int -> Poly v a -> Laurent v a
toLaurent (Int
off1 Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
off2) (Poly v a
poly1 Poly v a -> Poly v a -> Poly v a
forall a. Semiring a => a -> a -> a
`times` Poly v a
poly2)
  Laurent Int
off1 Poly v a
poly1 plus :: Laurent v a -> Laurent v a -> Laurent v a
`plus` Laurent Int
off2 Poly v a
poly2 = case Int
off1 Int -> Int -> Ordering
forall a. Ord a => a -> a -> Ordering
`compare` Int
off2 of
    Ordering
LT -> Int -> Poly v a -> Laurent v a
forall a (v :: * -> *).
(Eq a, Semiring a, Vector v a) =>
Int -> Poly v a -> Laurent v a
toLaurent Int
off1 (Poly v a
poly1 Poly v a -> Poly v a -> Poly v a
forall a. Semiring a => a -> a -> a
`plus` Word -> a -> Poly v a -> Poly v a
forall a (v :: * -> *).
(Eq a, Semiring a, Vector v a) =>
Word -> a -> Poly v a -> Poly v a
Dense.scale' (Int -> Word
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int -> Word) -> Int -> Word
forall a b. (a -> b) -> a -> b
$ Int
off2 Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
off1) a
forall a. Semiring a => a
one Poly v a
poly2)
    Ordering
EQ -> Int -> Poly v a -> Laurent v a
forall a (v :: * -> *).
(Eq a, Semiring a, Vector v a) =>
Int -> Poly v a -> Laurent v a
toLaurent Int
off1 (Poly v a
poly1 Poly v a -> Poly v a -> Poly v a
forall a. Semiring a => a -> a -> a
`plus` Poly v a
poly2)
    Ordering
GT -> Int -> Poly v a -> Laurent v a
forall a (v :: * -> *).
(Eq a, Semiring a, Vector v a) =>
Int -> Poly v a -> Laurent v a
toLaurent Int
off2 (Word -> a -> Poly v a -> Poly v a
forall a (v :: * -> *).
(Eq a, Semiring a, Vector v a) =>
Word -> a -> Poly v a -> Poly v a
Dense.scale' (Int -> Word
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int -> Word) -> Int -> Word
forall a b. (a -> b) -> a -> b
$ Int
off1 Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
off2) a
forall a. Semiring a => a
one Poly v a
poly1 Poly v a -> Poly v a -> Poly v a
forall a. Semiring a => a -> a -> a
`plus` Poly v a
poly2)
  fromNatural :: Natural -> Laurent v a
fromNatural Natural
n = Int -> Poly v a -> Laurent v a
forall (v :: * -> *) a. Int -> Poly v a -> Laurent v a
Laurent Int
0 (Natural -> Poly v a
forall a. Semiring a => Natural -> a
fromNatural Natural
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 v a -> Laurent v a
negate (Laurent Int
off Poly v a
poly) = Int -> Poly v a -> Laurent v a
forall (v :: * -> *) a. Int -> Poly v a -> Laurent v a
Laurent Int
off (Poly v a -> Poly v a
forall a. Ring a => a -> a
Semiring.negate Poly v a
poly)

-- | Create a monomial from a power and a coefficient.
--
-- @since 0.4.0.0
monomial :: (Eq a, Semiring a, G.Vector v a) => Int -> a -> Laurent v a
monomial :: forall a (v :: * -> *).
(Eq a, Semiring a, Vector v a) =>
Int -> a -> Laurent v a
monomial Int
p a
c
  | a
c a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
forall a. Semiring a => a
zero = Int -> Poly v a -> Laurent v a
forall (v :: * -> *) a. Int -> Poly v a -> Laurent v a
Laurent Int
0 Poly v a
forall a. Semiring a => a
zero
  | Bool
otherwise = Int -> Poly v a -> Laurent v a
forall (v :: * -> *) a. Int -> Poly v a -> Laurent v a
Laurent Int
p (Word -> a -> Poly v a
forall a (v :: * -> *).
(Eq a, Semiring a, Vector v a) =>
Word -> a -> Poly v a
Dense.monomial' Word
0 a
c)
{-# INLINE monomial #-}

-- | Multiply a polynomial by a monomial, expressed as a power and a coefficient.
--
-- >>> scale 2 3 (X^-2 + 1) :: ULaurent Int
-- 3 * X^2 + 0 * X + 3
--
-- @since 0.4.0.0
scale :: (Eq a, Semiring a, G.Vector v a) => Int -> a -> Laurent v a -> Laurent v a
scale :: forall a (v :: * -> *).
(Eq a, Semiring a, Vector v a) =>
Int -> a -> Laurent v a -> Laurent v a
scale Int
yp a
yc (Laurent Int
off Poly v a
poly) = Int -> Poly v a -> Laurent v a
forall a (v :: * -> *).
(Eq a, Semiring a, Vector v a) =>
Int -> Poly v a -> Laurent v a
toLaurent (Int
off Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
yp) (Word -> a -> Poly v a -> Poly v a
forall a (v :: * -> *).
(Eq a, Semiring a, Vector v a) =>
Word -> a -> Poly v a -> Poly v a
Dense.scale' Word
0 a
yc Poly v a
poly)
{-# INLINABLE scale #-}

-- | Evaluate the polynomial at a given point.
--
-- >>> eval (X^-2 + 1 :: ULaurent Double) 2
-- 1.25
--
-- @since 0.4.0.0
eval :: (Field a, G.Vector v a) => Laurent v a -> a -> a
eval :: forall a (v :: * -> *).
(Field a, Vector v a) =>
Laurent v a -> a -> a
eval (Laurent Int
off Poly v a
poly) a
x = Poly v a -> a -> a
forall a (v :: * -> *).
(Semiring a, Vector v a) =>
Poly v a -> a -> a
Dense.eval' Poly v a
poly a
x a -> a -> a
forall a. Semiring a => a -> a -> a
`times`
  (if Int
off Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
0 then a
x a -> Int -> a
forall a b. (Semiring a, Integral b) => a -> b -> a
Semiring.^ Int
off else a -> a -> a
forall a. Euclidean a => a -> a -> a
quot a
forall a. Semiring a => a
one a
x a -> Int -> a
forall a b. (Semiring a, Integral b) => a -> b -> a
Semiring.^ (- Int
off))
{-# INLINE eval #-}

-- | Substitute another polynomial instead of 'Data.Poly.X'.
--
-- >>> import Data.Poly (UPoly)
-- >>> subst (Data.Poly.X^2 + 1 :: UPoly Int) (X^-1 + 1 :: ULaurent Int)
-- 2 + 2 * X^-1 + 1 * X^-2
--
-- @since 0.4.0.0
subst :: (Eq a, Semiring a, G.Vector v a, G.Vector w a) => Poly v a -> Laurent w a -> Laurent w a
subst :: forall a (v :: * -> *) (w :: * -> *).
(Eq a, Semiring a, Vector v a, Vector w a) =>
Poly v a -> Laurent w a -> Laurent w a
subst = (a -> Laurent w a -> Laurent w a)
-> Poly v a -> Laurent w a -> Laurent w a
forall (v :: * -> *) a b.
(Vector v a, Semiring b) =>
(a -> b -> b) -> Poly v a -> b -> b
Dense.substitute' (Int -> a -> Laurent w a -> Laurent w a
forall a (v :: * -> *).
(Eq a, Semiring a, Vector v a) =>
Int -> a -> Laurent v a -> Laurent v a
scale Int
0)
{-# INLINE subst #-}

-- | Take the derivative of the polynomial.
--
-- >>> deriv (X^-1 + 3 * X) :: ULaurent Int
-- 3 + 0 * X^-1 + (-1) * X^-2
--
-- @since 0.4.0.0
deriv :: (Eq a, Ring a, G.Vector v a) => Laurent v a -> Laurent v a
deriv :: forall a (v :: * -> *).
(Eq a, Ring a, Vector v a) =>
Laurent v a -> Laurent v a
deriv (Laurent Int
off (Poly v a
xs)) =
  Int -> Poly v a -> Laurent v a
forall a (v :: * -> *).
(Eq a, Semiring a, Vector v a) =>
Int -> Poly v a -> Laurent v a
toLaurent (Int
off Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) (Poly v a -> Laurent v a) -> Poly v a -> Laurent v a
forall a b. (a -> b) -> a -> b
$ v a -> Poly v a
forall a (v :: * -> *).
(Eq a, Semiring a, Vector v a) =>
v a -> Poly v a
Dense.toPoly' (v a -> Poly v a) -> v a -> Poly v a
forall a b. (a -> b) -> a -> b
$ (Int -> a -> a) -> v a -> v a
forall (v :: * -> *) a b.
(Vector v a, Vector v b) =>
(Int -> a -> b) -> v a -> v b
G.imap (a -> a -> a
forall a. Semiring a => a -> a -> a
times (a -> a -> a) -> (Int -> a) -> Int -> a -> a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> a
forall a b. (Integral a, Ring b) => a -> b
Semiring.fromIntegral (Int -> a) -> (Int -> Int) -> Int -> a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
off)) v a
xs
{-# INLINE deriv #-}

-- | The polynomial 'X'.
--
-- > X == monomial 1 one
--
-- @since 0.4.0.0
pattern X :: (Eq a, Semiring a, G.Vector v a) => Laurent v a
pattern $mX :: forall {r} {a} {v :: * -> *}.
(Eq a, Semiring a, Vector v a) =>
Laurent v a -> ((# #) -> r) -> ((# #) -> r) -> r
$bX :: forall a (v :: * -> *).
(Eq a, Semiring a, Vector v a) =>
Laurent v a
X <- (isVar -> True)
  where X = Laurent v a
forall a (v :: * -> *).
(Eq a, Semiring a, Vector v a) =>
Laurent v a
var

var :: forall a v. (Eq a, Semiring a, G.Vector v a) => Laurent v a
var :: forall a (v :: * -> *).
(Eq a, Semiring a, Vector v a) =>
Laurent v a
var
  | (a
forall a. Semiring a => a
one :: a) a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
forall a. Semiring a => a
zero = Int -> Poly v a -> Laurent v a
forall (v :: * -> *) a. Int -> Poly v a -> Laurent v a
Laurent Int
0 Poly v a
forall a. Semiring a => a
zero
  | Bool
otherwise          = Int -> Poly v a -> Laurent v a
forall (v :: * -> *) a. Int -> Poly v a -> Laurent v a
Laurent Int
1 Poly v a
forall a. Semiring a => a
one
{-# INLINE var #-}

isVar :: forall v a. (Eq a, Semiring a, G.Vector v a) => Laurent v a -> Bool
isVar :: forall (v :: * -> *) a.
(Eq a, Semiring a, Vector v a) =>
Laurent v a -> Bool
isVar (Laurent Int
off (Poly v a
xs))
  | (a
forall a. Semiring a => a
one :: a) a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
forall a. Semiring a => a
zero = Int
off Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0 Bool -> Bool -> Bool
&& v a -> Bool
forall (v :: * -> *) a. Vector v a => v a -> Bool
G.null v a
xs
  | Bool
otherwise          = Int
off Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
1 Bool -> Bool -> Bool
&& v a -> Int
forall (v :: * -> *) a. Vector v a => v a -> Int
G.length v a
xs Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
1 Bool -> Bool -> Bool
&& v a -> a
forall (v :: * -> *) a. Vector v a => v a -> a
G.unsafeHead v a
xs a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
forall a. Semiring a => a
one
{-# INLINE isVar #-}

-- | Used to construct monomials with negative powers.
--
-- This operator can be applied only to monomials with unit coefficients,
-- but is instrumental to express Laurent polynomials
-- in a mathematical fashion:
--
-- >>> X^-3 :: ULaurent Int
-- 1 * X^-3
-- >>> X + 2 + 3 * (X^2)^-1 :: ULaurent Int
-- 1 * X + 2 + 0 * X^-1 + 3 * X^-2
--
-- @since 0.4.0.0
(^-)
  :: (Eq a, Num a, G.Vector v a)
  => Laurent v a
  -> Int
  -> Laurent v a
Laurent Int
off (Poly v a
xs) ^- :: forall a (v :: * -> *).
(Eq a, Num a, Vector v a) =>
Laurent v a -> Int -> Laurent v a
^- Int
n
  | v a -> Int
forall (v :: * -> *) a. Vector v a => v a -> Int
G.length v a
xs Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
1, v a -> a
forall (v :: * -> *) a. Vector v a => v a -> a
G.unsafeHead v a
xs a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
1
  = Int -> Poly v a -> Laurent v a
forall (v :: * -> *) a. Int -> Poly v a -> Laurent v a
Laurent (Int
off Int -> Int -> Int
forall a. Num a => a -> a -> a
* (-Int
n)) (v a -> Poly v a
forall (v :: * -> *) a. v a -> Poly v a
Poly v a
xs)
  | Bool
otherwise
  = PatternMatchFail -> Laurent v a
forall a e. Exception e => e -> a
throw (PatternMatchFail -> Laurent v a)
-> PatternMatchFail -> Laurent v a
forall a b. (a -> b) -> a -> b
$ String -> PatternMatchFail
PatternMatchFail String
"(^-) can be applied only to a monom with unit coefficient"

instance (Eq a, Ring a, GcdDomain a, G.Vector v a) => GcdDomain (Laurent v a) where
  divide :: Laurent v a -> Laurent v a -> Maybe (Laurent v a)
divide (Laurent Int
off1 Poly v a
poly1) (Laurent Int
off2 Poly v a
poly2) =
    Int -> Poly v a -> Laurent v a
forall a (v :: * -> *).
(Eq a, Semiring a, Vector v a) =>
Int -> Poly v a -> Laurent v a
toLaurent (Int
off1 Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
off2) (Poly v a -> Laurent v a)
-> Maybe (Poly v a) -> Maybe (Laurent v a)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Poly v a -> Poly v a -> Maybe (Poly v a)
forall a. GcdDomain a => a -> a -> Maybe a
divide Poly v a
poly1 Poly v a
poly2
  {-# INLINE divide #-}

  gcd :: Laurent v a -> Laurent v a -> Laurent v a
gcd (Laurent Int
_ Poly v a
poly1) (Laurent Int
_ Poly v a
poly2) =
    Int -> Poly v a -> Laurent v a
forall a (v :: * -> *).
(Eq a, Semiring a, Vector v a) =>
Int -> Poly v a -> Laurent v a
toLaurent Int
0 (Poly v a -> Poly v a -> Poly v a
forall a. GcdDomain a => a -> a -> a
gcd Poly v a
poly1 Poly v a
poly2)
  {-# INLINE gcd #-}

  lcm :: Laurent v a -> Laurent v a -> Laurent v a
lcm (Laurent Int
_ Poly v a
poly1) (Laurent Int
_ Poly v a
poly2) =
    Int -> Poly v a -> Laurent v a
forall a (v :: * -> *).
(Eq a, Semiring a, Vector v a) =>
Int -> Poly v a -> Laurent v a
toLaurent Int
0 (Poly v a -> Poly v a -> Poly v a
forall a. GcdDomain a => a -> a -> a
lcm Poly v a
poly1 Poly v a
poly2)
  {-# INLINE lcm #-}

  coprime :: Laurent v a -> Laurent v a -> Bool
coprime (Laurent Int
_ Poly v a
poly1) (Laurent Int
_ Poly v a
poly2) =
    Poly v a -> Poly v a -> Bool
forall a. GcdDomain a => a -> a -> Bool
coprime Poly v a
poly1 Poly v a
poly2
  {-# INLINE coprime #-}