-- |
-- Module:      Data.Poly.Internal.Dense.Field
-- Copyright:   (c) 2019 Andrew Lelechenko
-- Licence:     BSD3
-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>
--
-- 'Euclidean' instance with a 'Field' constraint on the coefficient type.
--

{-# LANGUAGE ConstraintKinds            #-}
{-# LANGUAGE FlexibleInstances          #-}
{-# LANGUAGE ScopedTypeVariables        #-}
{-# LANGUAGE TypeFamilies               #-}

{-# OPTIONS_GHC -fno-warn-orphans #-}

module Data.Poly.Internal.Dense.Field
  ( quotRemFractional
  ) where

import Prelude hiding (quotRem, quot, rem, gcd)
import Control.Exception
import Control.Monad
import Control.Monad.ST
import Data.Euclidean (Euclidean(..), Field)
import Data.Semiring (times, minus, zero, one)
import qualified Data.Vector.Generic as G
import qualified Data.Vector.Generic.Mutable as MG

import Data.Poly.Internal.Dense
import Data.Poly.Internal.Dense.GcdDomain ()

-- | Note that 'degree' 0 = 0.
--
-- @since 0.3.0.0
instance (Eq a, Field a, G.Vector v a) => Euclidean (Poly v a) where
  degree :: Poly v a -> Natural
degree (Poly v a
xs)
    | forall (v :: * -> *) a. Vector v a => v a -> Bool
G.null v a
xs = Natural
0
    | Bool
otherwise = forall a b. (Integral a, Num b) => a -> b
fromIntegral (forall (v :: * -> *) a. Vector v a => v a -> Int
G.length v a
xs forall a. Num a => a -> a -> a
- Int
1)

  quotRem :: Poly v a -> Poly v a -> (Poly v a, Poly v a)
quotRem (Poly v a
xs) (Poly v a
ys) = (forall a (v :: * -> *).
(Eq a, Semiring a, Vector v a) =>
v a -> Poly v a
toPoly' v a
qs, forall a (v :: * -> *).
(Eq a, Semiring a, Vector v a) =>
v a -> Poly v a
toPoly' v a
rs)
    where
      (v a
qs, v a
rs) = forall a (v :: * -> *).
(Eq a, Vector v a) =>
a
-> (a -> Bool)
-> (a -> a -> a)
-> (a -> a -> a)
-> (a -> a)
-> v a
-> v a
-> (v a, v a)
quotientAndRemainder forall a. Semiring a => a
zero (forall a. Eq a => a -> a -> Bool
== forall a. Semiring a => a
one) forall a. Ring a => a -> a -> a
minus forall a. Semiring a => a -> a -> a
times (forall a. Semiring a => a
one forall a. Euclidean a => a -> a -> a
`quot`) v a
xs v a
ys
  {-# INLINE quotRem #-}

  rem :: Poly v a -> Poly v a -> Poly v a
rem (Poly v a
xs) (Poly v a
ys) = forall a (v :: * -> *).
(Eq a, Semiring a, Vector v a) =>
v a -> Poly v a
toPoly' forall a b. (a -> b) -> a -> b
$ forall a (v :: * -> *).
(Eq a, Field a, Vector v a) =>
v a -> v a -> v a
remainder v a
xs v a
ys
  {-# INLINE rem #-}

-- | Polynomial division with remainder.
--
-- >>> quotRemFractional (X^3 + 2) (X^2 - 1 :: UPoly Double)
-- (1.0 * X + 0.0,1.0 * X + 2.0)
--
-- @since 0.5.0.0
quotRemFractional :: (Eq a, Fractional a, G.Vector v a) => Poly v a -> Poly v a -> (Poly v a, Poly v a)
quotRemFractional :: forall a (v :: * -> *).
(Eq a, Fractional a, Vector v a) =>
Poly v a -> Poly v a -> (Poly v a, Poly v a)
quotRemFractional (Poly v a
xs) (Poly v a
ys) = (forall a (v :: * -> *).
(Eq a, Num a, Vector v a) =>
v a -> Poly v a
toPoly v a
qs, forall a (v :: * -> *).
(Eq a, Num a, Vector v a) =>
v a -> Poly v a
toPoly v a
rs)
  where
    (v a
qs, v a
rs) = forall a (v :: * -> *).
(Eq a, Vector v a) =>
a
-> (a -> Bool)
-> (a -> a -> a)
-> (a -> a -> a)
-> (a -> a)
-> v a
-> v a
-> (v a, v a)
quotientAndRemainder a
0 (forall a. Eq a => a -> a -> Bool
== a
1) (-) forall a. Num a => a -> a -> a
(*) forall a. Fractional a => a -> a
recip v a
xs v a
ys
{-# INLINE quotRemFractional #-}

quotientAndRemainder
  :: (Eq a, G.Vector v a)
  => a             -- ^ zero
  -> (a -> Bool)   -- ^ is one?
  -> (a -> a -> a) -- ^ subtract
  -> (a -> a -> a) -- ^ multiply
  -> (a -> a)      -- ^ invert
  -> v a           -- ^ dividend
  -> v a           -- ^ divisor
  -> (v a, v a)
quotientAndRemainder :: forall a (v :: * -> *).
(Eq a, Vector v a) =>
a
-> (a -> Bool)
-> (a -> a -> a)
-> (a -> a -> a)
-> (a -> a)
-> v a
-> v a
-> (v a, v a)
quotientAndRemainder a
zer a -> Bool
isOne a -> a -> a
sub a -> a -> a
mul a -> a
inv v a
xs v a
ys
  | Int
lenXs forall a. Ord a => a -> a -> Bool
< Int
lenYs = (forall (v :: * -> *) a. Vector v a => v a
G.empty, v a
xs)
  | Int
lenYs forall a. Eq a => a -> a -> Bool
== Int
0 = forall a e. Exception e => e -> a
throw ArithException
DivideByZero
  | Int
lenYs forall a. Eq a => a -> a -> Bool
== Int
1 = let invY :: a
invY = a -> a
inv (forall (v :: * -> *) a. Vector v a => v a -> a
G.unsafeHead v a
ys) in
                 (forall (v :: * -> *) a b.
(Vector v a, Vector v b) =>
(a -> b) -> v a -> v b
G.map (a -> a -> a
`mul` a
invY) v a
xs, forall (v :: * -> *) a. Vector v a => v a
G.empty)
  | Bool
otherwise = forall a. (forall s. ST s a) -> a
runST forall a b. (a -> b) -> a -> b
$ do
    Mutable v s a
qs <- forall (m :: * -> *) (v :: * -> * -> *) a.
(PrimMonad m, MVector v a) =>
Int -> m (v (PrimState m) a)
MG.unsafeNew Int
lenQs
    Mutable v s a
rs <- forall (m :: * -> *) (v :: * -> * -> *) a.
(PrimMonad m, MVector v a) =>
Int -> m (v (PrimState m) a)
MG.unsafeNew Int
lenXs
    forall (m :: * -> *) (v :: * -> *) a.
(PrimMonad m, Vector v a) =>
Mutable v (PrimState m) a -> v a -> m ()
G.unsafeCopy Mutable v s a
rs v a
xs
    let yLast :: a
yLast = forall (v :: * -> *) a. Vector v a => v a -> a
G.unsafeLast v a
ys
        invYLast :: a
invYLast = a -> a
inv a
yLast
    forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [Int
lenQs forall a. Num a => a -> a -> a
- Int
1, Int
lenQs forall a. Num a => a -> a -> a
- Int
2 .. Int
0] forall a b. (a -> b) -> a -> b
$ \Int
i -> do
      a
r <- forall (m :: * -> *) (v :: * -> * -> *) a.
(PrimMonad m, MVector v a) =>
v (PrimState m) a -> Int -> m a
MG.unsafeRead Mutable v s a
rs (Int
lenYs forall a. Num a => a -> a -> a
- Int
1 forall a. Num a => a -> a -> a
+ Int
i)
      let q :: a
q = if a -> Bool
isOne a
yLast then a
r else a
r a -> a -> a
`mul` a
invYLast
      forall (m :: * -> *) (v :: * -> * -> *) a.
(PrimMonad m, MVector v a) =>
v (PrimState m) a -> Int -> a -> m ()
MG.unsafeWrite Mutable v s a
qs Int
i a
q
      forall (m :: * -> *) (v :: * -> * -> *) a.
(PrimMonad m, MVector v a) =>
v (PrimState m) a -> Int -> a -> m ()
MG.unsafeWrite Mutable v s a
rs (Int
lenYs forall a. Num a => a -> a -> a
- Int
1 forall a. Num a => a -> a -> a
+ Int
i) a
zer
      forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [Int
0 .. Int
lenYs forall a. Num a => a -> a -> a
- Int
2] forall a b. (a -> b) -> a -> b
$ \Int
k -> do
        let y :: a
y = forall (v :: * -> *) a. Vector v a => v a -> Int -> a
G.unsafeIndex v a
ys Int
k
        forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (a
y forall a. Eq a => a -> a -> Bool
/= a
zer) forall a b. (a -> b) -> a -> b
$
          forall (m :: * -> *) (v :: * -> * -> *) a.
(PrimMonad m, MVector v a) =>
v (PrimState m) a -> (a -> a) -> Int -> m ()
MG.unsafeModify Mutable v s a
rs (\a
c -> a
c a -> a -> a
`sub` (a
q a -> a -> a
`mul` a
y)) (Int
i forall a. Num a => a -> a -> a
+ Int
k)
    let rs' :: Mutable v s a
rs' = forall (v :: * -> * -> *) a s.
MVector v a =>
Int -> Int -> v s a -> v s a
MG.unsafeSlice Int
0 Int
lenYs Mutable v s a
rs
    (,) forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> forall (m :: * -> *) (v :: * -> *) a.
(PrimMonad m, Vector v a) =>
Mutable v (PrimState m) a -> m (v a)
G.unsafeFreeze Mutable v s a
qs forall (f :: * -> *) a b. Applicative f => f (a -> b) -> f a -> f b
<*> forall (m :: * -> *) (v :: * -> *) a.
(PrimMonad m, Vector v a) =>
Mutable v (PrimState m) a -> m (v a)
G.unsafeFreeze Mutable v s a
rs'
  where
    lenXs :: Int
lenXs = forall (v :: * -> *) a. Vector v a => v a -> Int
G.length v a
xs
    lenYs :: Int
lenYs = forall (v :: * -> *) a. Vector v a => v a -> Int
G.length v a
ys
    lenQs :: Int
lenQs = Int
lenXs forall a. Num a => a -> a -> a
- Int
lenYs forall a. Num a => a -> a -> a
+ Int
1
{-# INLINABLE quotientAndRemainder #-}

remainder
  :: (Eq a, Field a, G.Vector v a)
  => v a
  -> v a
  -> v a
remainder :: forall a (v :: * -> *).
(Eq a, Field a, Vector v a) =>
v a -> v a -> v a
remainder v a
xs v a
ys
  | forall (v :: * -> *) a. Vector v a => v a -> Bool
G.null v a
ys = forall a e. Exception e => e -> a
throw ArithException
DivideByZero
  | Bool
otherwise = forall a. (forall s. ST s a) -> a
runST forall a b. (a -> b) -> a -> b
$ do
    Mutable v s a
rs <- forall (m :: * -> *) (v :: * -> *) a.
(PrimMonad m, Vector v a) =>
v a -> m (Mutable v (PrimState m) a)
G.thaw v a
xs
    Mutable v s a
ys' <- forall (m :: * -> *) (v :: * -> *) a.
(PrimMonad m, Vector v a) =>
v a -> m (Mutable v (PrimState m) a)
G.unsafeThaw v a
ys
    forall a (v :: * -> *) s.
(Eq a, Field a, Vector v a) =>
Mutable v s a -> Mutable v s a -> ST s ()
remainderM Mutable v s a
rs Mutable v s a
ys'
    forall (m :: * -> *) (v :: * -> *) a.
(PrimMonad m, Vector v a) =>
Mutable v (PrimState m) a -> m (v a)
G.unsafeFreeze forall a b. (a -> b) -> a -> b
$ forall (v :: * -> * -> *) a s.
MVector v a =>
Int -> Int -> v s a -> v s a
MG.unsafeSlice Int
0 (forall (v :: * -> *) a. Vector v a => v a -> Int
G.length v a
xs forall a. Ord a => a -> a -> a
`min` forall (v :: * -> *) a. Vector v a => v a -> Int
G.length v a
ys) Mutable v s a
rs
{-# INLINABLE remainder #-}

remainderM
  :: (Eq a, Field a, G.Vector v a)
  => G.Mutable v s a
  -> G.Mutable v s a
  -> ST s ()
remainderM :: forall a (v :: * -> *) s.
(Eq a, Field a, Vector v a) =>
Mutable v s a -> Mutable v s a -> ST s ()
remainderM Mutable v s a
xs Mutable v s a
ys
  | Int
lenXs forall a. Ord a => a -> a -> Bool
< Int
lenYs = forall (f :: * -> *) a. Applicative f => a -> f a
pure ()
  | Int
lenYs forall a. Eq a => a -> a -> Bool
== Int
0 = forall a e. Exception e => e -> a
throw ArithException
DivideByZero
  | Int
lenYs forall a. Eq a => a -> a -> Bool
== Int
1 = forall (m :: * -> *) (v :: * -> * -> *) a.
(PrimMonad m, MVector v a) =>
v (PrimState m) a -> a -> m ()
MG.set Mutable v s a
xs forall a. Semiring a => a
zero
  | Bool
otherwise = do
    a
yLast <- forall (m :: * -> *) (v :: * -> * -> *) a.
(PrimMonad m, MVector v a) =>
v (PrimState m) a -> Int -> m a
MG.unsafeRead Mutable v s a
ys (Int
lenYs forall a. Num a => a -> a -> a
- Int
1)
    let invYLast :: a
invYLast = forall a. Semiring a => a
one forall a. Euclidean a => a -> a -> a
`quot` a
yLast
    forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [Int
lenQs forall a. Num a => a -> a -> a
- Int
1, Int
lenQs forall a. Num a => a -> a -> a
- Int
2 .. Int
0] forall a b. (a -> b) -> a -> b
$ \Int
i -> do
      a
r <- forall (m :: * -> *) (v :: * -> * -> *) a.
(PrimMonad m, MVector v a) =>
v (PrimState m) a -> Int -> m a
MG.unsafeRead Mutable v s a
xs (Int
lenYs forall a. Num a => a -> a -> a
- Int
1 forall a. Num a => a -> a -> a
+ Int
i)
      forall (m :: * -> *) (v :: * -> * -> *) a.
(PrimMonad m, MVector v a) =>
v (PrimState m) a -> Int -> a -> m ()
MG.unsafeWrite Mutable v s a
xs (Int
lenYs forall a. Num a => a -> a -> a
- Int
1 forall a. Num a => a -> a -> a
+ Int
i) forall a. Semiring a => a
zero
      let q :: a
q = if a
yLast forall a. Eq a => a -> a -> Bool
== forall a. Semiring a => a
one then a
r else a
r forall a. Semiring a => a -> a -> a
`times` a
invYLast
      forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [Int
0 .. Int
lenYs forall a. Num a => a -> a -> a
- Int
2] forall a b. (a -> b) -> a -> b
$ \Int
k -> do
        a
y <- forall (m :: * -> *) (v :: * -> * -> *) a.
(PrimMonad m, MVector v a) =>
v (PrimState m) a -> Int -> m a
MG.unsafeRead Mutable v s a
ys Int
k
        forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (a
y forall a. Eq a => a -> a -> Bool
/= forall a. Semiring a => a
zero) forall a b. (a -> b) -> a -> b
$
          forall (m :: * -> *) (v :: * -> * -> *) a.
(PrimMonad m, MVector v a) =>
v (PrimState m) a -> (a -> a) -> Int -> m ()
MG.unsafeModify Mutable v s a
xs (\a
c -> a
c forall a. Ring a => a -> a -> a
`minus` a
q forall a. Semiring a => a -> a -> a
`times` a
y) (Int
i forall a. Num a => a -> a -> a
+ Int
k)
  where
    lenXs :: Int
lenXs = forall (v :: * -> * -> *) a s. MVector v a => v s a -> Int
MG.length Mutable v s a
xs
    lenYs :: Int
lenYs = forall (v :: * -> * -> *) a s. MVector v a => v s a -> Int
MG.length Mutable v s a
ys
    lenQs :: Int
lenQs = Int
lenXs forall a. Num a => a -> a -> a
- Int
lenYs forall a. Num a => a -> a -> a
+ Int
1
{-# INLINABLE remainderM #-}