content: -- |
  -- Module:      Data.Poly.Semiring
  -- Copyright:   (c) 2019 Andrew Lelechenko
  -- Licence:     BSD3
  -- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>
  --
  -- Dense polynomials and a 'Semiring'-based interface.
  --
  -- @since 0.2.0.0
  
  {-# LANGUAGE CPP              #-}
  {-# LANGUAGE DataKinds        #-}
  {-# LANGUAGE FlexibleContexts #-}
  {-# LANGUAGE PatternSynonyms  #-}
  
  {-# OPTIONS_GHC -Wno-unrecognised-warning-flags -Wno-pattern-namespace-specifier #-}
  
  module Data.Poly.Semiring
    ( Poly
    , VPoly
    , UPoly
    , unPoly
    , leading
    , toPoly
    , monomial
    , scale
    , pattern X
    , eval
    , subst
    , deriv
    , integral
    , timesRing
  #ifdef SupportSparse
    , denseToSparse
    , sparseToDense
  #endif
    , dft
    , inverseDft
    , dftMult
    ) where
  
  import Data.Bits
  import Data.Euclidean (Field)
  import Data.Semiring (Semiring(..))
  import qualified Data.Vector.Generic as G
  
  import Data.Poly.Internal.Dense (Poly(..), VPoly, UPoly, leading, timesRing)
  import qualified Data.Poly.Internal.Dense as Dense
  import Data.Poly.Internal.Dense.Field ()
  import Data.Poly.Internal.Dense.DFT
  import Data.Poly.Internal.Dense.GcdDomain ()
  
  #ifdef SupportSparse
  import qualified Data.Poly.Internal.Multi as Sparse
  import qualified Data.Poly.Internal.Convert as Convert
  #endif
  
  -- | Make a 'Poly' from a vector of coefficients
  -- (first element corresponds to the constant term).
  --
  -- >>> :set -XOverloadedLists
  -- >>> toPoly [1,2,3] :: VPoly Integer
  -- 3 * X^2 + 2 * X + 1
  -- >>> toPoly [0,0,0] :: UPoly Int
  -- 0
  --
  -- @since 0.2.0.0
  toPoly :: (Eq a, Semiring a, G.Vector v a) => v a -> Poly v a
  toPoly = Dense.toPoly'
  
  -- | Create a monomial from a power and a coefficient.
  --
  -- @since 0.3.0.0
  monomial :: (Eq a, Semiring a, G.Vector v a) => Word -> a -> Poly v a
  monomial = Dense.monomial'
  
  -- | Multiply a polynomial by a monomial, expressed as a power and a coefficient.
  --
  -- >>> scale 2 3 (X^2 + 1) :: UPoly Int
  -- 3 * X^4 + 0 * X^3 + 3 * X^2 + 0 * X + 0
  --
  -- @since 0.3.0.0
  scale :: (Eq a, Semiring a, G.Vector v a) => Word -> a -> Poly v a -> Poly v a
  scale = Dense.scale'
  
  -- | The polynomial 'X'.
  --
  -- > X == monomial 1 one
  --
  -- @since 0.2.0.0
  pattern X :: (Eq a, Semiring a, G.Vector v a) => Poly v a
  pattern X = Dense.X'
  
  -- | Evaluate the polynomial at a given point.
  --
  -- >>> eval (X^2 + 1 :: UPoly Int) 3
  -- 10
  --
  -- @since 0.2.0.0
  eval :: (Semiring a, G.Vector v a) => Poly v a -> a -> a
  eval = Dense.eval'
  
  -- | Substitute another polynomial instead of 'X'.
  --
  -- >>> subst (X^2 + 1 :: UPoly Int) (X + 1 :: UPoly Int)
  -- 1 * X^2 + 2 * X + 2
  --
  -- @since 0.3.3.0
  subst :: (Eq a, Semiring a, G.Vector v a, G.Vector w a) => Poly v a -> Poly w a -> Poly w a
  subst = Dense.subst'
  
  -- | Take the derivative of the polynomial.
  --
  -- >>> deriv (X^3 + 3 * X) :: UPoly Int
  -- 3 * X^2 + 0 * X + 3
  --
  -- @since 0.2.0.0
  deriv :: (Eq a, Semiring a, G.Vector v a) => Poly v a -> Poly v a
  deriv = Dense.deriv'
  
  -- | Compute an indefinite integral of the polynomial,
  -- setting the constant term to zero.
  --
  -- >>> integral (3 * X^2 + 3) :: UPoly Double
  -- 1.0 * X^3 + 0.0 * X^2 + 3.0 * X + 0.0
  --
  -- @since 0.3.2.0
  integral :: (Eq a, Field a, G.Vector v a) => Poly v a -> Poly v a
  integral = Dense.integral'
  
  -- | Multiplication of polynomials using
  -- <https://en.wikipedia.org/wiki/Fast_Fourier_transform discrete Fourier transform>.
  -- It could be faster than '(*)' for large polynomials
  -- if multiplication of coefficients is particularly expensive.
  --
  -- @since 0.5.0.0
  dftMult
    :: (Eq a, Field a, G.Vector v a)
    => (Int -> a) -- ^ mapping from \( N = 2^n \) to a primitive root \( \sqrt[N]{1} \)
    -> Poly v a
    -> Poly v a
    -> Poly v a
  dftMult getPrimRoot (Poly xs) (Poly ys) =
    toPoly $ inverseDft primRoot $ G.zipWith times (dft primRoot xs') (dft primRoot ys')
    where
      nextPowerOf2 :: Int -> Int
      nextPowerOf2 0 = 1
      nextPowerOf2 1 = 1
      nextPowerOf2 x = 1 `unsafeShiftL` (finiteBitSize (0 :: Int) - countLeadingZeros (x - 1))
  
      padTo l vs = G.generate l (\k -> if k < G.length vs then vs G.! k else zero)
  
      zl = nextPowerOf2 (G.length xs + G.length ys)
      xs' = padTo zl xs
      ys' = padTo zl ys
      primRoot = getPrimRoot zl
  {-# INLINABLE dftMult #-}
  
  #ifdef SupportSparse
  -- | Convert from dense to sparse polynomials.
  --
  -- >>> :set -XFlexibleContexts
  -- >>> denseToSparse (1 `Data.Semiring.plus` Data.Poly.X^2) :: Data.Poly.Sparse.UPoly Int
  -- 1 * X^2 + 1
  --
  -- @since 0.5.0.0
  denseToSparse :: (Eq a, Semiring a, G.Vector v a, G.Vector v (Sparse.Monom 1 a)) => Dense.Poly v a -> Sparse.Poly v a
  denseToSparse = Convert.denseToSparse'
  
  -- | Convert from sparse to dense polynomials.
  --
  -- >>> :set -XFlexibleContexts
  -- >>> sparseToDense (1 `Data.Semiring.plus` Data.Poly.Sparse.X^2) :: Data.Poly.UPoly Int
  -- 1 * X^2 + 0 * X + 1
  --
  -- @since 0.5.0.0
  sparseToDense :: (Semiring a, G.Vector v a, G.Vector v (Sparse.Monom 1 a)) => Sparse.Poly v a -> Dense.Poly v a
  sparseToDense = Convert.sparseToDense'
  #endif
  