module Numeric.Domain.Euclidean (Euclidean(..), prs, normalize, gcd', leadingUnit, chineseRemainder) where
import Numeric.Additive.Group
import Numeric.Algebra.Class
import Numeric.Algebra.Unital
import Numeric.Decidable.Units
import Numeric.Decidable.Zero
import Numeric.Domain.Class
import Numeric.Natural (Natural)
import Numeric.Ring.Class
import Prelude (Eq (..), Integer, Maybe (..), abs)
import Prelude (fst, otherwise)
import Prelude (signum, snd, ($), (.))
import qualified Prelude as P
infixl 7 `quot`, `rem`
infix 7 `divide`
class (Ring r, DecidableZero r, DecidableUnits r, Domain r) => Euclidean r where
splitUnit :: r -> (r, r)
degree :: r -> Maybe Natural
divide :: r
-> r
-> (r,r)
quot :: r -> r -> r
quot a b = fst $ a `divide` b
rem :: r -> r -> r
rem a b = snd $ a `divide` b
gcd :: r -> r -> r
gcd a b = let (g,_,_):_ = euclid a b
in g
euclid :: r -> r -> [(r,r,r)]
euclid f g =
let (ug, g') = splitUnit g
Just t' = recipUnit ug
(uf, f') = splitUnit f
Just s = recipUnit uf
in step [(g', 0, t'), (f', s, 0)]
where
step acc@((r',s',t'):(r,s,t):_)
| isZero r' = P.tail acc
| otherwise =
let q = r `quot` r'
(ur, r'') = splitUnit $ r q * r'
Just u = recipUnit ur
s'' = (s q * s') * u
t'' = (t q * t') * u
in step ((r'', s'', t'') : acc)
step _ = P.error "cannot happen!"
#if (__GLASGOW_HASKELL__ > 708)
#endif
prs :: Euclidean r => r -> r -> [(r, r, r)]
prs f g = step [(g, 0, 1), (f, 1, 0)]
where
step acc@((r',s',t'):(r,s,t):_)
| isZero r' = P.tail acc
| otherwise =
let q = r `quot` r'
s'' = (s q * s')
t'' = (t q * t')
in step ((r q * r', s'', t'') : acc)
step _ = P.error "cannot happen!"
gcd' :: Euclidean r => [r] -> r
gcd' [] = one
gcd' [x] = leadingUnit x
gcd' [x,y] = gcd x y
gcd' (x:xs) = gcd x (gcd' xs)
normalize :: Euclidean r => r -> r
normalize = snd . splitUnit
leadingUnit :: Euclidean r => r -> r
leadingUnit = fst . splitUnit
instance Euclidean Integer where
splitUnit 0 = (1, 0)
splitUnit n = (signum n, abs n)
degree = Just . P.fromInteger . abs
divide = P.divMod
chineseRemainder :: Euclidean r
=> [(r, r)]
-> r
chineseRemainder mvs =
let (ms, _) = P.unzip mvs
m = product ms
in sum [((vi*s) `rem` mi)*n | (mi, vi) <- mvs
, let n = m `quot` mi
, let (_, s, _) : _ = euclid n mi
]