content: -- |
  -- 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 Data.Bifunctor
  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 (Monom 1 a)) => Euclidean (Poly v a) where
    degree (MultiPoly xs)
      | G.null xs = 0
      | otherwise = fromIntegral (SU.head (monomPower (G.unsafeLast xs)))
  
    quotRem = quotientRemainder zero plus minus times 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 (Monom 1 a)) => Poly v a -> Poly v a -> (Poly v a, Poly v a)
  quotRemFractional = quotientRemainder 0 (+) (-) (*) (/)
  {-# INLINE quotRemFractional #-}
  
  quotientRemainder
    :: G.Vector v (Monom 1 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 zer add sub mul div ts ys = case leading ys of
    Nothing -> throw DivideByZero
    Just (yp, yc) -> go ts
      where
        go xs = case leading xs of
          Nothing -> (zer, zer)
          Just (xp, xc) -> case xp `compare` yp of
            LT -> (zer, xs)
            EQ -> (zs, xs')
            GT -> first (`add` zs) $ go xs'
            where
              zs = MultiPoly $ G.singleton (Monom (SU.singleton (xp - yp)) (xc `div` yc))
              xs' = xs `sub` (zs `mul` ys)
  