-- |
-- Module:      Data.Poly.Internal.Multi.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 DataKinds                  #-}
{-# LANGUAGE FlexibleContexts           #-}
{-# LANGUAGE FlexibleInstances          #-}
{-# LANGUAGE ScopedTypeVariables        #-}
{-# LANGUAGE TypeFamilies               #-}
{-# LANGUAGE UndecidableInstances       #-}

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

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

import Prelude hiding (quotRem, quot, rem, div, gcd)
import Control.Arrow
import Control.Exception
import Data.Euclidean (Euclidean(..), Field)
import Data.Semiring (Semiring(..), minus)
import qualified Data.Vector.Generic as G
import qualified Data.Vector.Unboxed.Sized as SU

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

-- | Note that 'degree' 0 = 0.
instance (Eq a, Field a, G.Vector v (SU.Vector 1 Word, a)) => Euclidean (Poly v a) where
  degree :: Poly v a -> Natural
degree (MultiPoly v (Vector 1 Word, a)
xs)
    | forall (v :: * -> *) a. Vector v a => v a -> Bool
G.null v (Vector 1 Word, a)
xs = Natural
0
    | Bool
otherwise = forall a b. (Integral a, Num b) => a -> b
fromIntegral (forall (n :: Natural) a. Unbox a => Vector (1 + n) a -> a
SU.head (forall a b. (a, b) -> a
fst (forall (v :: * -> *) a. Vector v a => v a -> a
G.unsafeLast v (Vector 1 Word, a)
xs)))

  quotRem :: Poly v a -> Poly v a -> (Poly v a, Poly v a)
quotRem = forall (v :: * -> *) a.
Vector v (Vector 1 Word, a) =>
Poly v a
-> (Poly v a -> Poly v a -> Poly v a)
-> (Poly v a -> Poly v a -> Poly v a)
-> (Poly v a -> Poly v a -> Poly v a)
-> (a -> a -> a)
-> Poly v a
-> Poly v a
-> (Poly v a, Poly v a)
quotientRemainder forall a. Semiring a => a
zero forall a. Semiring a => a -> a -> a
plus forall a. Ring a => a -> a -> a
minus forall a. Semiring a => a -> a -> a
times forall a. Euclidean a => a -> a -> a
quot

-- | Polynomial division with remainder.
--
-- >>> quotRemFractional (X^3 + 2) (X^2 - 1 :: UPoly Double)
-- (1.0 * X,1.0 * X + 2.0)
--
-- @since 0.5.0.0
quotRemFractional :: (Eq a, Fractional a, G.Vector v (SU.Vector 1 Word, a)) => Poly v a -> Poly v a -> (Poly v a, Poly v a)
quotRemFractional :: forall a (v :: * -> *).
(Eq a, Fractional a, Vector v (Vector 1 Word, a)) =>
Poly v a -> Poly v a -> (Poly v a, Poly v a)
quotRemFractional = forall (v :: * -> *) a.
Vector v (Vector 1 Word, a) =>
Poly v a
-> (Poly v a -> Poly v a -> Poly v a)
-> (Poly v a -> Poly v a -> Poly v a)
-> (Poly v a -> Poly v a -> Poly v a)
-> (a -> a -> a)
-> Poly v a
-> Poly v a
-> (Poly v a, Poly v a)
quotientRemainder Poly v a
0 forall a. Num a => a -> a -> a
(+) (-) forall a. Num a => a -> a -> a
(*) forall a. Fractional a => a -> a -> a
(/)
{-# INLINE quotRemFractional #-}

quotientRemainder
  :: G.Vector v (SU.Vector 1 Word, a)
  => Poly v a                           -- ^ zero
  -> (Poly v a -> Poly v a -> Poly v a) -- ^ add
  -> (Poly v a -> Poly v a -> Poly v a) -- ^ subtract
  -> (Poly v a -> Poly v a -> Poly v a) -- ^ multiply
  -> (a -> a -> a)                      -- ^ divide
  -> Poly v a                           -- ^ dividend
  -> Poly v a                           -- ^ divisor
  -> (Poly v a, Poly v a)
quotientRemainder :: forall (v :: * -> *) a.
Vector v (Vector 1 Word, a) =>
Poly v a
-> (Poly v a -> Poly v a -> Poly v a)
-> (Poly v a -> Poly v a -> Poly v a)
-> (Poly v a -> Poly v a -> Poly v a)
-> (a -> a -> a)
-> Poly v a
-> Poly v a
-> (Poly v a, Poly v a)
quotientRemainder Poly v a
zer Poly v a -> Poly v a -> Poly v a
add Poly v a -> Poly v a -> Poly v a
sub Poly v a -> Poly v a -> Poly v a
mul a -> a -> a
div Poly v a
ts Poly v a
ys = case forall (v :: * -> *) a.
Vector v (Vector 1 Word, a) =>
Poly v a -> Maybe (Word, a)
leading Poly v a
ys of
  Maybe (Word, a)
Nothing -> forall a e. Exception e => e -> a
throw ArithException
DivideByZero
  Just (Word
yp, a
yc) -> Poly v a -> (Poly v a, Poly v a)
go Poly v a
ts
    where
      go :: Poly v a -> (Poly v a, Poly v a)
go Poly v a
xs = case forall (v :: * -> *) a.
Vector v (Vector 1 Word, a) =>
Poly v a -> Maybe (Word, a)
leading Poly v a
xs of
        Maybe (Word, a)
Nothing -> (Poly v a
zer, Poly v a
zer)
        Just (Word
xp, a
xc) -> case Word
xp forall a. Ord a => a -> a -> Ordering
`compare` Word
yp of
          Ordering
LT -> (Poly v a
zer, Poly v a
xs)
          Ordering
EQ -> (Poly v a
zs, Poly v a
xs')
          Ordering
GT -> forall (a :: * -> * -> *) b c d.
Arrow a =>
a b c -> a (b, d) (c, d)
first (Poly v a -> Poly v a -> Poly v a
`add` Poly v a
zs) forall a b. (a -> b) -> a -> b
$ Poly v a -> (Poly v a, Poly v a)
go Poly v a
xs'
          where
            zs :: Poly v a
zs = forall (v :: * -> *) (n :: Natural) a.
v (Vector n Word, a) -> MultiPoly v n a
MultiPoly forall a b. (a -> b) -> a -> b
$ forall (v :: * -> *) a. Vector v a => a -> v a
G.singleton (forall a. Unbox a => a -> Vector 1 a
SU.singleton (Word
xp forall a. Num a => a -> a -> a
- Word
yp), a
xc a -> a -> a
`div` a
yc)
            xs' :: Poly v a
xs' = Poly v a
xs Poly v a -> Poly v a -> Poly v a
`sub` (Poly v a
zs Poly v a -> Poly v a -> Poly v a
`mul` Poly v a
ys)