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

{-# LANGUAGE CPP                        #-}
{-# LANGUAGE DataKinds                  #-}
{-# LANGUAGE FlexibleContexts           #-}
{-# LANGUAGE FlexibleInstances          #-}
{-# LANGUAGE GADTs                      #-}
{-# LANGUAGE QuantifiedConstraints      #-}
{-# LANGUAGE ScopedTypeVariables        #-}
{-# LANGUAGE TypeOperators              #-}
{-# LANGUAGE TypeFamilies               #-}
{-# LANGUAGE UndecidableInstances       #-}

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

module Data.Poly.Internal.Multi.GcdDomain
  () where

import Prelude hiding (gcd, lcm, (^))
import Control.Exception
import Data.Euclidean
import Data.Maybe
import Data.Proxy
import Data.Semiring (Semiring(..), Ring(), minus)
import Data.Type.Equality
import qualified Data.Vector.Generic as G
import qualified Data.Vector.Unboxed.Sized as SU
import GHC.TypeNats (KnownNat, type (+), SomeNat(..), natVal, sameNat, someNatVal)
import Unsafe.Coerce

import Data.Poly.Internal.Multi

instance {-# OVERLAPPING #-} (Eq a, Ring a, GcdDomain a, G.Vector v (SU.Vector 1 Word, a)) => GcdDomain (Poly v a) where
  divide :: Poly v a -> Poly v a -> Maybe (Poly v a)
divide Poly v a
xs Poly v a
ys
    | v (Vector 1 Word, a) -> Bool
forall (v :: * -> *) a. Vector v a => v a -> Bool
G.null (Poly v a -> v (Vector 1 Word, a)
forall (v :: * -> *) (n :: Nat) a.
MultiPoly v n a -> v (Vector n Word, a)
unMultiPoly Poly v a
ys) = ArithException -> Maybe (Poly v a)
forall a e. Exception e => e -> a
throw ArithException
DivideByZero
    | v (Vector 1 Word, a) -> Int
forall (v :: * -> *) a. Vector v a => v a -> Int
G.length (Poly v a -> v (Vector 1 Word, a)
forall (v :: * -> *) (n :: Nat) a.
MultiPoly v n a -> v (Vector n Word, a)
unMultiPoly Poly v a
ys) Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
1 = Poly v a -> (Vector 1 Word, a) -> Maybe (Poly v a)
forall a (v :: * -> *) (n :: Nat).
(GcdDomain a, Vector v (Vector n Word, a)) =>
MultiPoly v n a -> (Vector n Word, a) -> Maybe (MultiPoly v n a)
divideSingleton Poly v a
xs (v (Vector 1 Word, a) -> (Vector 1 Word, a)
forall (v :: * -> *) a. Vector v a => v a -> a
G.unsafeHead (Poly v a -> v (Vector 1 Word, a)
forall (v :: * -> *) (n :: Nat) a.
MultiPoly v n a -> v (Vector n Word, a)
unMultiPoly Poly v a
ys))
    | Bool
otherwise = Poly v a -> Poly v a -> Maybe (Poly v a)
forall a (v :: * -> *).
(Eq a, GcdDomain a, Ring a, Vector v (Vector 1 Word, a)) =>
Poly v a -> Poly v a -> Maybe (Poly v a)
divide1 Poly v a
xs Poly v a
ys

  gcd :: Poly v a -> Poly v a -> Poly v a
gcd Poly v a
xs Poly v a
ys
    | v (Vector 1 Word, a) -> Bool
forall (v :: * -> *) a. Vector v a => v a -> Bool
G.null (Poly v a -> v (Vector 1 Word, a)
forall (v :: * -> *) (n :: Nat) a.
MultiPoly v n a -> v (Vector n Word, a)
unMultiPoly Poly v a
xs) = Poly v a
ys
    | v (Vector 1 Word, a) -> Bool
forall (v :: * -> *) a. Vector v a => v a -> Bool
G.null (Poly v a -> v (Vector 1 Word, a)
forall (v :: * -> *) (n :: Nat) a.
MultiPoly v n a -> v (Vector n Word, a)
unMultiPoly Poly v a
ys) = Poly v a
xs
    | v (Vector 1 Word, a) -> Int
forall (v :: * -> *) a. Vector v a => v a -> Int
G.length (Poly v a -> v (Vector 1 Word, a)
forall (v :: * -> *) (n :: Nat) a.
MultiPoly v n a -> v (Vector n Word, a)
unMultiPoly Poly v a
xs) Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
1 = (Vector 1 Word, a) -> Poly v a -> Poly v a
forall a (v :: * -> *) (n :: Nat).
(Eq a, GcdDomain a, Vector v (Vector n Word, a)) =>
(Vector n Word, a) -> MultiPoly v n a -> MultiPoly v n a
gcdSingleton (v (Vector 1 Word, a) -> (Vector 1 Word, a)
forall (v :: * -> *) a. Vector v a => v a -> a
G.unsafeHead (Poly v a -> v (Vector 1 Word, a)
forall (v :: * -> *) (n :: Nat) a.
MultiPoly v n a -> v (Vector n Word, a)
unMultiPoly Poly v a
xs)) Poly v a
ys
    | v (Vector 1 Word, a) -> Int
forall (v :: * -> *) a. Vector v a => v a -> Int
G.length (Poly v a -> v (Vector 1 Word, a)
forall (v :: * -> *) (n :: Nat) a.
MultiPoly v n a -> v (Vector n Word, a)
unMultiPoly Poly v a
ys) Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
1 = (Vector 1 Word, a) -> Poly v a -> Poly v a
forall a (v :: * -> *) (n :: Nat).
(Eq a, GcdDomain a, Vector v (Vector n Word, a)) =>
(Vector n Word, a) -> MultiPoly v n a -> MultiPoly v n a
gcdSingleton (v (Vector 1 Word, a) -> (Vector 1 Word, a)
forall (v :: * -> *) a. Vector v a => v a -> a
G.unsafeHead (Poly v a -> v (Vector 1 Word, a)
forall (v :: * -> *) (n :: Nat) a.
MultiPoly v n a -> v (Vector n Word, a)
unMultiPoly Poly v a
ys)) Poly v a
xs
    | Bool
otherwise = Poly v a -> Poly v a -> Poly v a
forall a (v :: * -> *).
(Eq a, GcdDomain a, Ring a, Vector v (Vector 1 Word, a)) =>
Poly v a -> Poly v a -> Poly v a
gcd1 Poly v a
xs Poly v a
ys

  lcm :: Poly v a -> Poly v a -> Poly v a
lcm Poly v a
xs Poly v a
ys
    | v (Vector 1 Word, a) -> Bool
forall (v :: * -> *) a. Vector v a => v a -> Bool
G.null (Poly v a -> v (Vector 1 Word, a)
forall (v :: * -> *) (n :: Nat) a.
MultiPoly v n a -> v (Vector n Word, a)
unMultiPoly Poly v a
xs) Bool -> Bool -> Bool
|| v (Vector 1 Word, a) -> Bool
forall (v :: * -> *) a. Vector v a => v a -> Bool
G.null (Poly v a -> v (Vector 1 Word, a)
forall (v :: * -> *) (n :: Nat) a.
MultiPoly v n a -> v (Vector n Word, a)
unMultiPoly Poly v a
ys) = Poly v a
forall a. Semiring a => a
zero
    | Bool
otherwise = (Poly v a
xs Poly v a -> Poly v a -> Poly v a
forall a. GcdDomain a => a -> a -> a
`divide'` Poly v a -> Poly v a -> Poly v a
forall a. GcdDomain a => a -> a -> a
gcd Poly v a
xs Poly v a
ys) Poly v a -> Poly v a -> Poly v a
forall a. Semiring a => a -> a -> a
`times` Poly v a
ys

  coprime :: Poly v a -> Poly v a -> Bool
coprime Poly v a
x Poly v a
y = Maybe (Poly v a) -> Bool
forall a. Maybe a -> Bool
isJust (Poly v a
forall a. Semiring a => a
one Poly v a -> Poly v a -> Maybe (Poly v a)
forall a. GcdDomain a => a -> a -> Maybe a
`divide` Poly v a -> Poly v a -> Poly v a
forall a. GcdDomain a => a -> a -> a
gcd Poly v a
x Poly v a
y)

data IsSucc n where
  IsSucc :: KnownNat m => n :~: 1 + m -> IsSucc n

-- | This is unsafe when n ~ 0.
isSucc :: forall n. KnownNat n => IsSucc n
isSucc :: forall (n :: Nat). KnownNat n => IsSucc n
isSucc = case Nat -> SomeNat
someNatVal (Proxy n -> Nat
forall (n :: Nat) (proxy :: Nat -> *). KnownNat n => proxy n -> Nat
natVal (Proxy n
forall {k} (t :: k). Proxy t
Proxy :: Proxy n) Nat -> Nat -> Nat
forall a. Num a => a -> a -> a
- Nat
1) of
  SomeNat (Proxy n
_ :: Proxy m) -> (n :~: (1 + n)) -> IsSucc n
forall (m :: Nat) (n :: Nat).
KnownNat m =>
(n :~: (1 + m)) -> IsSucc n
IsSucc ((Any :~: Any) -> n :~: (1 + n)
forall a b. a -> b
unsafeCoerce Any :~: Any
forall {k} (a :: k). a :~: a
Refl :: n :~: 1 + m)

instance (Eq a, Ring a, GcdDomain a, KnownNat n, forall m. KnownNat m => G.Vector v (SU.Vector m Word, a), forall m. KnownNat m => Eq (v (SU.Vector m Word, a))) => GcdDomain (MultiPoly v n a) where
  divide :: MultiPoly v n a -> MultiPoly v n a -> Maybe (MultiPoly v n a)
divide MultiPoly v n a
xs MultiPoly v n a
ys
    | v (Vector n Word, a) -> Bool
forall (v :: * -> *) a. Vector v a => v a -> Bool
G.null (MultiPoly v n a -> v (Vector n Word, a)
forall (v :: * -> *) (n :: Nat) a.
MultiPoly v n a -> v (Vector n Word, a)
unMultiPoly MultiPoly v n a
ys) = ArithException -> Maybe (MultiPoly v n a)
forall a e. Exception e => e -> a
throw ArithException
DivideByZero
    | v (Vector n Word, a) -> Int
forall (v :: * -> *) a. Vector v a => v a -> Int
G.length (MultiPoly v n a -> v (Vector n Word, a)
forall (v :: * -> *) (n :: Nat) a.
MultiPoly v n a -> v (Vector n Word, a)
unMultiPoly MultiPoly v n a
ys) Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
1 = MultiPoly v n a -> (Vector n Word, a) -> Maybe (MultiPoly v n a)
forall a (v :: * -> *) (n :: Nat).
(GcdDomain a, Vector v (Vector n Word, a)) =>
MultiPoly v n a -> (Vector n Word, a) -> Maybe (MultiPoly v n a)
divideSingleton MultiPoly v n a
xs (v (Vector n Word, a) -> (Vector n Word, a)
forall (v :: * -> *) a. Vector v a => v a -> a
G.unsafeHead (MultiPoly v n a -> v (Vector n Word, a)
forall (v :: * -> *) (n :: Nat) a.
MultiPoly v n a -> v (Vector n Word, a)
unMultiPoly MultiPoly v n a
ys))
    -- Polynomials of zero variables are necessarily constants,
    -- so they have been dealt with above.
    | Just n :~: 1
Refl <- Proxy n -> Proxy 1 -> Maybe (n :~: 1)
forall (a :: Nat) (b :: Nat) (proxy1 :: Nat -> *)
       (proxy2 :: Nat -> *).
(KnownNat a, KnownNat b) =>
proxy1 a -> proxy2 b -> Maybe (a :~: b)
sameNat (Proxy n
forall {k} (t :: k). Proxy t
Proxy :: Proxy n) (Proxy 1
forall {k} (t :: k). Proxy t
Proxy :: Proxy 1)
    = Poly v a -> Poly v a -> Maybe (Poly v a)
forall a (v :: * -> *).
(Eq a, GcdDomain a, Ring a, Vector v (Vector 1 Word, a)) =>
Poly v a -> Poly v a -> Maybe (Poly v a)
divide1 MultiPoly v n a
Poly v a
xs MultiPoly v n a
Poly v a
ys
    | Bool
otherwise = case IsSucc n
forall (n :: Nat). KnownNat n => IsSucc n
isSucc :: IsSucc n of
      IsSucc n :~: (1 + m)
Refl -> VPoly (MultiPoly v m a) -> MultiPoly v n a
VPoly (MultiPoly v m a) -> MultiPoly v (1 + m) a
forall (v :: * -> *) (m :: Nat) a.
(Vector v (Vector (1 + m) Word, a), Vector v (Vector m Word, a)) =>
VPoly (MultiPoly v m a) -> MultiPoly v (1 + m) a
unsegregate (VPoly (MultiPoly v m a) -> MultiPoly v n a)
-> Maybe (VPoly (MultiPoly v m a)) -> Maybe (MultiPoly v n a)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> MultiPoly v (1 + m) a -> VPoly (MultiPoly v m a)
forall (v :: * -> *) (m :: Nat) a.
(Vector v (Vector (1 + m) Word, a), Vector v (Vector m Word, a)) =>
MultiPoly v (1 + m) a -> VPoly (MultiPoly v m a)
segregate MultiPoly v n a
MultiPoly v (1 + m) a
xs VPoly (MultiPoly v m a)
-> VPoly (MultiPoly v m a) -> Maybe (VPoly (MultiPoly v m a))
forall a. GcdDomain a => a -> a -> Maybe a
`divide` MultiPoly v (1 + m) a -> VPoly (MultiPoly v m a)
forall (v :: * -> *) (m :: Nat) a.
(Vector v (Vector (1 + m) Word, a), Vector v (Vector m Word, a)) =>
MultiPoly v (1 + m) a -> VPoly (MultiPoly v m a)
segregate MultiPoly v n a
MultiPoly v (1 + m) a
ys
  gcd :: MultiPoly v n a -> MultiPoly v n a -> MultiPoly v n a
gcd MultiPoly v n a
xs MultiPoly v n a
ys
    | v (Vector n Word, a) -> Bool
forall (v :: * -> *) a. Vector v a => v a -> Bool
G.null (MultiPoly v n a -> v (Vector n Word, a)
forall (v :: * -> *) (n :: Nat) a.
MultiPoly v n a -> v (Vector n Word, a)
unMultiPoly MultiPoly v n a
xs) = MultiPoly v n a
ys
    | v (Vector n Word, a) -> Bool
forall (v :: * -> *) a. Vector v a => v a -> Bool
G.null (MultiPoly v n a -> v (Vector n Word, a)
forall (v :: * -> *) (n :: Nat) a.
MultiPoly v n a -> v (Vector n Word, a)
unMultiPoly MultiPoly v n a
ys) = MultiPoly v n a
xs
    | v (Vector n Word, a) -> Int
forall (v :: * -> *) a. Vector v a => v a -> Int
G.length (MultiPoly v n a -> v (Vector n Word, a)
forall (v :: * -> *) (n :: Nat) a.
MultiPoly v n a -> v (Vector n Word, a)
unMultiPoly MultiPoly v n a
xs) Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
1 = (Vector n Word, a) -> MultiPoly v n a -> MultiPoly v n a
forall a (v :: * -> *) (n :: Nat).
(Eq a, GcdDomain a, Vector v (Vector n Word, a)) =>
(Vector n Word, a) -> MultiPoly v n a -> MultiPoly v n a
gcdSingleton (v (Vector n Word, a) -> (Vector n Word, a)
forall (v :: * -> *) a. Vector v a => v a -> a
G.unsafeHead (MultiPoly v n a -> v (Vector n Word, a)
forall (v :: * -> *) (n :: Nat) a.
MultiPoly v n a -> v (Vector n Word, a)
unMultiPoly MultiPoly v n a
xs)) MultiPoly v n a
ys
    | v (Vector n Word, a) -> Int
forall (v :: * -> *) a. Vector v a => v a -> Int
G.length (MultiPoly v n a -> v (Vector n Word, a)
forall (v :: * -> *) (n :: Nat) a.
MultiPoly v n a -> v (Vector n Word, a)
unMultiPoly MultiPoly v n a
ys) Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
1 = (Vector n Word, a) -> MultiPoly v n a -> MultiPoly v n a
forall a (v :: * -> *) (n :: Nat).
(Eq a, GcdDomain a, Vector v (Vector n Word, a)) =>
(Vector n Word, a) -> MultiPoly v n a -> MultiPoly v n a
gcdSingleton (v (Vector n Word, a) -> (Vector n Word, a)
forall (v :: * -> *) a. Vector v a => v a -> a
G.unsafeHead (MultiPoly v n a -> v (Vector n Word, a)
forall (v :: * -> *) (n :: Nat) a.
MultiPoly v n a -> v (Vector n Word, a)
unMultiPoly MultiPoly v n a
ys)) MultiPoly v n a
xs
    -- Polynomials of zero variables are necessarily constants,
    -- so they have been dealt with above.
    | Just n :~: 1
Refl <- Proxy n -> Proxy 1 -> Maybe (n :~: 1)
forall (a :: Nat) (b :: Nat) (proxy1 :: Nat -> *)
       (proxy2 :: Nat -> *).
(KnownNat a, KnownNat b) =>
proxy1 a -> proxy2 b -> Maybe (a :~: b)
sameNat (Proxy n
forall {k} (t :: k). Proxy t
Proxy :: Proxy n) (Proxy 1
forall {k} (t :: k). Proxy t
Proxy :: Proxy 1)
    = Poly v a -> Poly v a -> Poly v a
forall a (v :: * -> *).
(Eq a, GcdDomain a, Ring a, Vector v (Vector 1 Word, a)) =>
Poly v a -> Poly v a -> Poly v a
gcd1 MultiPoly v n a
Poly v a
xs MultiPoly v n a
Poly v a
ys
    | Bool
otherwise = case IsSucc n
forall (n :: Nat). KnownNat n => IsSucc n
isSucc :: IsSucc n of
      IsSucc n :~: (1 + m)
Refl -> VPoly (MultiPoly v m a) -> MultiPoly v (1 + m) a
forall (v :: * -> *) (m :: Nat) a.
(Vector v (Vector (1 + m) Word, a), Vector v (Vector m Word, a)) =>
VPoly (MultiPoly v m a) -> MultiPoly v (1 + m) a
unsegregate (VPoly (MultiPoly v m a) -> MultiPoly v (1 + m) a)
-> VPoly (MultiPoly v m a) -> MultiPoly v (1 + m) a
forall a b. (a -> b) -> a -> b
$ MultiPoly v (1 + m) a -> VPoly (MultiPoly v m a)
forall (v :: * -> *) (m :: Nat) a.
(Vector v (Vector (1 + m) Word, a), Vector v (Vector m Word, a)) =>
MultiPoly v (1 + m) a -> VPoly (MultiPoly v m a)
segregate MultiPoly v n a
MultiPoly v (1 + m) a
xs VPoly (MultiPoly v m a)
-> VPoly (MultiPoly v m a) -> VPoly (MultiPoly v m a)
forall a. GcdDomain a => a -> a -> a
`gcd` MultiPoly v (1 + m) a -> VPoly (MultiPoly v m a)
forall (v :: * -> *) (m :: Nat) a.
(Vector v (Vector (1 + m) Word, a), Vector v (Vector m Word, a)) =>
MultiPoly v (1 + m) a -> VPoly (MultiPoly v m a)
segregate MultiPoly v n a
MultiPoly v (1 + m) a
ys

divideSingleton
  :: (GcdDomain a, G.Vector v (SU.Vector n Word, a))
  => MultiPoly v n a
  -> (SU.Vector n Word, a)
  -> Maybe (MultiPoly v n a)
divideSingleton :: forall a (v :: * -> *) (n :: Nat).
(GcdDomain a, Vector v (Vector n Word, a)) =>
MultiPoly v n a -> (Vector n Word, a) -> Maybe (MultiPoly v n a)
divideSingleton (MultiPoly v (Vector n Word, a)
pcs) (Vector n Word
p, a
c) = v (Vector n Word, a) -> MultiPoly v n a
forall (v :: * -> *) (n :: Nat) a.
v (Vector n Word, a) -> MultiPoly v n a
MultiPoly (v (Vector n Word, a) -> MultiPoly v n a)
-> Maybe (v (Vector n Word, a)) -> Maybe (MultiPoly v n a)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> ((Vector n Word, a) -> Maybe (Vector n Word, a))
-> v (Vector n Word, a) -> Maybe (v (Vector n Word, a))
forall (m :: * -> *) (v :: * -> *) a b.
(Monad m, Vector v a, Vector v b) =>
(a -> m b) -> v a -> m (v b)
G.mapM (Vector n Word, a) -> Maybe (Vector n Word, a)
divideMonomial v (Vector n Word, a)
pcs
  where
    divideMonomial :: (Vector n Word, a) -> Maybe (Vector n Word, a)
divideMonomial (Vector n Word
p', a
c')
      | Vector n Bool -> Bool
forall (n :: Nat). Vector n Bool -> Bool
SU.and ((Word -> Word -> Bool)
-> Vector n Word -> Vector n Word -> Vector n Bool
forall a b c (n :: Nat).
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector n a -> Vector n b -> Vector n c
SU.zipWith Word -> Word -> Bool
forall a. Ord a => a -> a -> Bool
(>=) Vector n Word
p' Vector n Word
p)
      , Just a
c'' <- a
c' a -> a -> Maybe a
forall a. GcdDomain a => a -> a -> Maybe a
`divide` a
c
      = (Vector n Word, a) -> Maybe (Vector n Word, a)
forall a. a -> Maybe a
Just ((Word -> Word -> Word)
-> Vector n Word -> Vector n Word -> Vector n Word
forall a b c (n :: Nat).
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector n a -> Vector n b -> Vector n c
SU.zipWith (-) Vector n Word
p' Vector n Word
p, a
c'')
      | Bool
otherwise
      = Maybe (Vector n Word, a)
forall a. Maybe a
Nothing

gcdSingleton
  :: (Eq a, GcdDomain a, G.Vector v (SU.Vector n Word, a))
  => (SU.Vector n Word, a)
  -> MultiPoly v n a
  -> MultiPoly v n a
gcdSingleton :: forall a (v :: * -> *) (n :: Nat).
(Eq a, GcdDomain a, Vector v (Vector n Word, a)) =>
(Vector n Word, a) -> MultiPoly v n a -> MultiPoly v n a
gcdSingleton (Vector n Word, a)
pc (MultiPoly v (Vector n Word, a)
pcs) = (Vector n Word -> a -> MultiPoly v n a)
-> (Vector n Word, a) -> MultiPoly v n a
forall a b c. (a -> b -> c) -> (a, b) -> c
uncurry Vector n Word -> a -> MultiPoly v n a
forall a (v :: * -> *) (n :: Nat).
(Eq a, Semiring a, Vector v (Vector n Word, a)) =>
Vector n Word -> a -> MultiPoly v n a
monomial' ((Vector n Word, a) -> MultiPoly v n a)
-> (Vector n Word, a) -> MultiPoly v n a
forall a b. (a -> b) -> a -> b
$
  ((Vector n Word, a) -> (Vector n Word, a) -> (Vector n Word, a))
-> (Vector n Word, a) -> v (Vector n Word, a) -> (Vector n Word, a)
forall (v :: * -> *) b a.
Vector v b =>
(a -> b -> a) -> a -> v b -> a
G.foldl' (\(Vector n Word
accP, a
accC) (Vector n Word
p, a
c) -> ((Word -> Word -> Word)
-> Vector n Word -> Vector n Word -> Vector n Word
forall a b c (n :: Nat).
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector n a -> Vector n b -> Vector n c
SU.zipWith Word -> Word -> Word
forall a. Ord a => a -> a -> a
min Vector n Word
accP Vector n Word
p, a -> a -> a
forall a. GcdDomain a => a -> a -> a
gcd a
accC a
c)) (Vector n Word, a)
pc v (Vector n Word, a)
pcs

divide1
  :: (Eq a, GcdDomain a, Ring a, G.Vector v (SU.Vector 1 Word, a))
  => Poly v a
  -> Poly v a
  -> Maybe (Poly v a)
divide1 :: forall a (v :: * -> *).
(Eq a, GcdDomain a, Ring a, Vector v (Vector 1 Word, a)) =>
Poly v a -> Poly v a -> Maybe (Poly v a)
divide1 Poly v a
xs Poly v a
ys = case Poly v a -> Maybe (Word, a)
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 -> ArithException -> Maybe (Poly v a)
forall a e. Exception e => e -> a
throw ArithException
DivideByZero
  Just (Word
yp, a
yc) -> case Poly v a -> Maybe (Word, a)
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 -> Maybe (Poly v a)
forall a. a -> Maybe a
Just Poly v a
xs
    Just (Word
xp, a
xc)
      | Word
xp Word -> Word -> Bool
forall a. Ord a => a -> a -> Bool
< Word
yp -> Maybe (Poly v a)
forall a. Maybe a
Nothing
      | Bool
otherwise -> do
        a
zc <- a -> a -> Maybe a
forall a. GcdDomain a => a -> a -> Maybe a
divide a
xc a
yc
        let z :: Poly v a
z = v (Vector 1 Word, a) -> Poly v a
forall (v :: * -> *) (n :: Nat) a.
v (Vector n Word, a) -> MultiPoly v n a
MultiPoly (v (Vector 1 Word, a) -> Poly v a)
-> v (Vector 1 Word, a) -> Poly v a
forall a b. (a -> b) -> a -> b
$ (Vector 1 Word, a) -> v (Vector 1 Word, a)
forall (v :: * -> *) a. Vector v a => a -> v a
G.singleton (Word -> Vector 1 Word
forall a. Unbox a => a -> Vector 1 a
SU.singleton (Word
xp Word -> Word -> Word
forall a. Num a => a -> a -> a
- Word
yp), a
zc)
        Poly v a
rest <- Poly v a -> Poly v a -> Maybe (Poly v a)
forall a (v :: * -> *).
(Eq a, GcdDomain a, Ring a, Vector v (Vector 1 Word, a)) =>
Poly v a -> Poly v a -> Maybe (Poly v a)
divide1 (Poly v a
xs Poly v a -> Poly v a -> Poly v a
forall a. Ring a => a -> a -> a
`minus` Poly v a
z Poly v a -> Poly v a -> Poly v a
forall a. Semiring a => a -> a -> a
`times` Poly v a
ys) Poly v a
ys
        Poly v a -> Maybe (Poly v a)
forall a. a -> Maybe a
forall (f :: * -> *) a. Applicative f => a -> f a
pure (Poly v a -> Maybe (Poly v a)) -> Poly v a -> Maybe (Poly v a)
forall a b. (a -> b) -> a -> b
$ Poly v a
rest Poly v a -> Poly v a -> Poly v a
forall a. Semiring a => a -> a -> a
`plus` Poly v a
z

gcd1
  :: (Eq a, GcdDomain a, Ring a, G.Vector v (SU.Vector 1 Word, a))
  => Poly v a
  -> Poly v a
  -> Poly v a
gcd1 :: forall a (v :: * -> *).
(Eq a, GcdDomain a, Ring a, Vector v (Vector 1 Word, a)) =>
Poly v a -> Poly v a -> Poly v a
gcd1 x :: Poly v a
x@(MultiPoly v (Vector 1 Word, a)
xs) y :: Poly v a
y@(MultiPoly v (Vector 1 Word, a)
ys) =
  Poly v a -> Poly v a -> Poly v a
forall a. Semiring a => a -> a -> a
times Poly v a
xy (Poly v a -> Poly v a -> Poly v a
divide1' Poly v a
z (Vector 1 Word -> a -> Poly v a
forall a (v :: * -> *) (n :: Nat).
(Eq a, Semiring a, Vector v (Vector n Word, a)) =>
Vector n Word -> a -> MultiPoly v n a
monomial' Vector 1 Word
0 (v (Vector 1 Word, a) -> a
forall a (v :: * -> *) t.
(GcdDomain a, Vector v (t, a)) =>
v (t, a) -> a
content v (Vector 1 Word, a)
zs)))
  where
    z :: Poly v a
z@(MultiPoly v (Vector 1 Word, a)
zs) = Poly v a -> Poly v a -> Poly v a
forall a (v :: * -> *).
(Eq a, Ring a, GcdDomain a, Vector v (Vector 1 Word, a)) =>
Poly v a -> Poly v a -> Poly v a
gcdHelper Poly v a
x Poly v a
y
    xy :: Poly v a
xy = Vector 1 Word -> a -> Poly v a
forall a (v :: * -> *) (n :: Nat).
(Eq a, Semiring a, Vector v (Vector n Word, a)) =>
Vector n Word -> a -> MultiPoly v n a
monomial' Vector 1 Word
0 (a -> a -> a
forall a. GcdDomain a => a -> a -> a
gcd (v (Vector 1 Word, a) -> a
forall a (v :: * -> *) t.
(GcdDomain a, Vector v (t, a)) =>
v (t, a) -> a
content v (Vector 1 Word, a)
xs) (v (Vector 1 Word, a) -> a
forall a (v :: * -> *) t.
(GcdDomain a, Vector v (t, a)) =>
v (t, a) -> a
content v (Vector 1 Word, a)
ys))
    divide1' :: Poly v a -> Poly v a -> Poly v a
divide1' = (Poly v a -> Maybe (Poly v a) -> Poly v a
forall a. a -> Maybe a -> a
fromMaybe ([Char] -> Poly v a
forall a. HasCallStack => [Char] -> a
error [Char]
"gcd: violated internal invariant") (Maybe (Poly v a) -> Poly v a)
-> (Poly v a -> Maybe (Poly v a)) -> Poly v a -> Poly v a
forall b c a. (b -> c) -> (a -> b) -> a -> c
.) ((Poly v a -> Maybe (Poly v a)) -> Poly v a -> Poly v a)
-> (Poly v a -> Poly v a -> Maybe (Poly v a))
-> Poly v a
-> Poly v a
-> Poly v a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Poly v a -> Poly v a -> Maybe (Poly v a)
forall a (v :: * -> *).
(Eq a, GcdDomain a, Ring a, Vector v (Vector 1 Word, a)) =>
Poly v a -> Poly v a -> Maybe (Poly v a)
divide1

content :: (GcdDomain a, G.Vector v (t, a)) => v (t, a) -> a
content :: forall a (v :: * -> *) t.
(GcdDomain a, Vector v (t, a)) =>
v (t, a) -> a
content = (a -> (t, a) -> a) -> a -> v (t, a) -> a
forall (v :: * -> *) b a.
Vector v b =>
(a -> b -> a) -> a -> v b -> a
G.foldl' (\a
acc (t
_, a
t) -> a -> a -> a
forall a. GcdDomain a => a -> a -> a
gcd a
acc a
t) a
forall a. Semiring a => a
zero

gcdHelper
  :: (Eq a, Ring a, GcdDomain a, G.Vector v (SU.Vector 1 Word, a))
  => Poly v a
  -> Poly v a
  -> Poly v a
gcdHelper :: forall a (v :: * -> *).
(Eq a, Ring a, GcdDomain a, Vector v (Vector 1 Word, a)) =>
Poly v a -> Poly v a -> Poly v a
gcdHelper Poly v a
xs Poly v a
ys = case (Poly v a -> Maybe (Word, a)
forall (v :: * -> *) a.
Vector v (Vector 1 Word, a) =>
Poly v a -> Maybe (Word, a)
leading Poly v a
xs, Poly v a -> Maybe (Word, a)
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, Maybe (Word, a)
_) -> Poly v a
ys
  (Maybe (Word, a)
_, Maybe (Word, a)
Nothing) -> Poly v a
xs
  (Just (Word
xp, a
xc), Just (Word
yp, a
yc))
    | Word
yp Word -> Word -> Bool
forall a. Ord a => a -> a -> Bool
<= Word
xp
    , Just a
xy <- a
xc a -> a -> Maybe a
forall a. GcdDomain a => a -> a -> Maybe a
`divide` a
yc
    -> Poly v a -> Poly v a -> Poly v a
forall a (v :: * -> *).
(Eq a, Ring a, GcdDomain a, Vector v (Vector 1 Word, a)) =>
Poly v a -> Poly v a -> Poly v a
gcdHelper Poly v a
ys (Poly v a
xs Poly v a -> Poly v a -> Poly v a
forall a. Ring a => a -> a -> a
`minus` Poly v a
ys Poly v a -> Poly v a -> Poly v a
forall a. Semiring a => a -> a -> a
`times` Vector 1 Word -> a -> Poly v a
forall a (v :: * -> *) (n :: Nat).
(Eq a, Semiring a, Vector v (Vector n Word, a)) =>
Vector n Word -> a -> MultiPoly v n a
monomial' (Word -> Vector 1 Word
forall a. Unbox a => a -> Vector 1 a
SU.singleton (Word
xp Word -> Word -> Word
forall a. Num a => a -> a -> a
- Word
yp)) a
xy)
    | Word
xp Word -> Word -> Bool
forall a. Ord a => a -> a -> Bool
<= Word
yp
    , Just a
yx <- a
yc a -> a -> Maybe a
forall a. GcdDomain a => a -> a -> Maybe a
`divide` a
xc
    -> Poly v a -> Poly v a -> Poly v a
forall a (v :: * -> *).
(Eq a, Ring a, GcdDomain a, Vector v (Vector 1 Word, a)) =>
Poly v a -> Poly v a -> Poly v a
gcdHelper Poly v a
xs (Poly v a
ys Poly v a -> Poly v a -> Poly v a
forall a. Ring a => a -> a -> a
`minus` Poly v a
xs Poly v a -> Poly v a -> Poly v a
forall a. Semiring a => a -> a -> a
`times` Vector 1 Word -> a -> Poly v a
forall a (v :: * -> *) (n :: Nat).
(Eq a, Semiring a, Vector v (Vector n Word, a)) =>
Vector n Word -> a -> MultiPoly v n a
monomial' (Word -> Vector 1 Word
forall a. Unbox a => a -> Vector 1 a
SU.singleton (Word
yp Word -> Word -> Word
forall a. Num a => a -> a -> a
- Word
xp)) a
yx)
    | Word
yp Word -> Word -> Bool
forall a. Ord a => a -> a -> Bool
<= Word
xp
    -> Poly v a -> Poly v a -> Poly v a
forall a (v :: * -> *).
(Eq a, Ring a, GcdDomain a, Vector v (Vector 1 Word, a)) =>
Poly v a -> Poly v a -> Poly v a
gcdHelper Poly v a
ys (Poly v a
xs Poly v a -> Poly v a -> Poly v a
forall a. Semiring a => a -> a -> a
`times` Vector 1 Word -> a -> Poly v a
forall a (v :: * -> *) (n :: Nat).
(Eq a, Semiring a, Vector v (Vector n Word, a)) =>
Vector n Word -> a -> MultiPoly v n a
monomial' Vector 1 Word
0 a
gx Poly v a -> Poly v a -> Poly v a
forall a. Ring a => a -> a -> a
`minus` Poly v a
ys Poly v a -> Poly v a -> Poly v a
forall a. Semiring a => a -> a -> a
`times` Vector 1 Word -> a -> Poly v a
forall a (v :: * -> *) (n :: Nat).
(Eq a, Semiring a, Vector v (Vector n Word, a)) =>
Vector n Word -> a -> MultiPoly v n a
monomial' (Word -> Vector 1 Word
forall a. Unbox a => a -> Vector 1 a
SU.singleton (Word
xp Word -> Word -> Word
forall a. Num a => a -> a -> a
- Word
yp)) a
gy)
    | Bool
otherwise
    -> Poly v a -> Poly v a -> Poly v a
forall a (v :: * -> *).
(Eq a, Ring a, GcdDomain a, Vector v (Vector 1 Word, a)) =>
Poly v a -> Poly v a -> Poly v a
gcdHelper Poly v a
xs (Poly v a
ys Poly v a -> Poly v a -> Poly v a
forall a. Semiring a => a -> a -> a
`times` Vector 1 Word -> a -> Poly v a
forall a (v :: * -> *) (n :: Nat).
(Eq a, Semiring a, Vector v (Vector n Word, a)) =>
Vector n Word -> a -> MultiPoly v n a
monomial' Vector 1 Word
0 a
gy Poly v a -> Poly v a -> Poly v a
forall a. Ring a => a -> a -> a
`minus` Poly v a
xs Poly v a -> Poly v a -> Poly v a
forall a. Semiring a => a -> a -> a
`times` Vector 1 Word -> a -> Poly v a
forall a (v :: * -> *) (n :: Nat).
(Eq a, Semiring a, Vector v (Vector n Word, a)) =>
Vector n Word -> a -> MultiPoly v n a
monomial' (Word -> Vector 1 Word
forall a. Unbox a => a -> Vector 1 a
SU.singleton (Word
yp Word -> Word -> Word
forall a. Num a => a -> a -> a
- Word
xp)) a
gx)
    where
      g :: a
g = a -> a -> a
forall a. GcdDomain a => a -> a -> a
lcm a
xc a
yc
      gx :: a
gx = a -> a -> a
forall a. GcdDomain a => a -> a -> a
divide' a
g a
xc
      gy :: a
gy = a -> a -> a
forall a. GcdDomain a => a -> a -> a
divide' a
g a
yc

divide' :: GcdDomain a => a -> a -> a
divide' :: forall a. GcdDomain a => a -> a -> a
divide' = (a -> Maybe a -> a
forall a. a -> Maybe a -> a
fromMaybe ([Char] -> a
forall a. HasCallStack => [Char] -> a
error [Char]
"gcd: violated internal invariant") (Maybe a -> a) -> (a -> Maybe a) -> a -> a
forall b c a. (b -> c) -> (a -> b) -> a -> c
.) ((a -> Maybe a) -> a -> a) -> (a -> a -> Maybe a) -> a -> a -> a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. a -> a -> Maybe a
forall a. GcdDomain a => a -> a -> Maybe a
divide