Expose GMP's `mpz_gcdext()` as internal primitive
authorHerbert Valerio Riedel <hvr@gnu.org>
Sun, 29 Sep 2013 16:25:07 +0000 (18:25 +0200)
committerHerbert Valerio Riedel <hvr@gnu.org>
Sun, 29 Sep 2013 16:25:07 +0000 (18:25 +0200)
The extended GCD computation is useful to have for implementing
algorithms such as the chinese reminder theorem.

Signed-off-by: Herbert Valerio Riedel <hvr@gnu.org>
GHC/Integer/GMP/Internals.hs
GHC/Integer/GMP/Prim.hs
GHC/Integer/Type.lhs
cbits/gmp-wrappers.cmm

index 4a7ff5d..b80840b 100644 (file)
@@ -1,6 +1,6 @@
 {-# LANGUAGE NoImplicitPrelude #-}
 
-module GHC.Integer.GMP.Internals (Integer(..), gcdInt, gcdInteger, lcmInteger, powInteger, powModInteger, recipModInteger)
+module GHC.Integer.GMP.Internals (Integer(..), gcdInt, gcdInteger, gcdExtInteger, lcmInteger, powInteger, powModInteger, recipModInteger)
     where
 
 import GHC.Integer.Type
index 59aa6f4..de9477f 100644 (file)
@@ -21,6 +21,7 @@ module GHC.Integer.GMP.Prim (
     divExactInteger#,
 
     gcdInteger#,
+    gcdExtInteger#,
     gcdIntegerInt#,
     gcdInt#,
 
@@ -124,6 +125,11 @@ foreign import prim "integer_cmm_divExactIntegerzh" divExactInteger#
 foreign import prim "integer_cmm_gcdIntegerzh" gcdInteger#
   :: Int# -> ByteArray# -> Int# -> ByteArray# -> (# Int#, ByteArray# #)
 
+-- | Extended greatest common divisor.
+--
+foreign import prim "integer_cmm_gcdExtIntegerzh" gcdExtInteger#
+  :: Int# -> ByteArray# -> Int# -> ByteArray# -> (# Int#, ByteArray#, Int#, ByteArray# #)
+
 -- | Greatest common divisor, where second argument is an ordinary {\tt Int\#}.
 --
 foreign import prim "integer_cmm_gcdIntegerIntzh" gcdIntegerInt#
index 554160c..3d4994a 100644 (file)
@@ -40,7 +40,7 @@ import GHC.Integer.GMP.Prim (
     plusInteger#, minusInteger#, timesInteger#,
     quotRemInteger#, quotInteger#, remInteger#,
     divModInteger#, divInteger#, modInteger#,
-    gcdInteger#, gcdIntegerInt#, gcdInt#, divExactInteger#,
+    gcdInteger#, gcdExtInteger#, gcdIntegerInt#, gcdInt#, divExactInteger#,
     decodeDouble#,
     int2Integer#, integer2Int#, word2Integer#, integer2Word#,
     andInteger#, orInteger#, xorInteger#, complementInteger#,
@@ -272,6 +272,17 @@ gcdInteger ia@(J# _ _) ib@(S# _) = gcdInteger ib ia
 gcdInteger (J# sa a) (J# sb b)
   = case gcdInteger# sa a sb b of (# sg, g #) -> J# sg g
 
+-- | For a and b, compute their greatest common divisor g and the
+-- coefficient s satisfying @a*s + b*t = g@.
+{-# NOINLINE gcdExtInteger #-}
+gcdExtInteger :: Integer -> Integer -> (# Integer, Integer #)
+gcdExtInteger a@(S# _)   b@(S# _) = gcdExtInteger (toBig a) (toBig b)
+gcdExtInteger a@(S# _) b@(J# _ _) = gcdExtInteger (toBig a) b
+gcdExtInteger a@(J# _ _) b@(S# _) = gcdExtInteger a (toBig b)
+gcdExtInteger (J# sa a) (J# sb b)
+  = case gcdExtInteger# sa a sb b of
+      (# sg, g, ss, s #) -> (# J# sg g, J# ss s #)
+
 {-# NOINLINE lcmInteger #-}
 lcmInteger :: Integer -> Integer -> Integer
 lcmInteger a b =      if a `eqInteger` S# 0# then S# 0#
index 5c7bb0b..68e6485 100644 (file)
@@ -36,6 +36,7 @@ import "integer-gmp" __gmpz_mul_2exp;
 import "integer-gmp" __gmpz_tstbit;
 import "integer-gmp" __gmpz_fdiv_q_2exp;
 import "integer-gmp" __gmpz_gcd;
+import "integer-gmp" __gmpz_gcdext;
 import "integer-gmp" __gmpn_gcd_1;
 import "integer-gmp" __gmpn_cmp;
 import "integer-gmp" __gmpz_tdiv_q;
@@ -417,6 +418,8 @@ GMP_TAKE2_RET1(integer_cmm_plusIntegerzh,           __gmpz_add)
 GMP_TAKE2_RET1(integer_cmm_minusIntegerzh,          __gmpz_sub)
 GMP_TAKE2_RET1(integer_cmm_timesIntegerzh,          __gmpz_mul)
 GMP_TAKE2_RET1(integer_cmm_gcdIntegerzh,            __gmpz_gcd)
+#define CMM_GMPZ_GCDEXT(g,s,a,b) __gmpz_gcdext(g,s,NULL,a,b)
+GMP_TAKE2_RET2(integer_cmm_gcdExtIntegerzh,         CMM_GMPZ_GCDEXT)
 GMP_TAKE2_RET1(integer_cmm_quotIntegerzh,           __gmpz_tdiv_q)
 GMP_TAKE2_RET1(integer_cmm_remIntegerzh,            __gmpz_tdiv_r)
 GMP_TAKE2_RET1(integer_cmm_divIntegerzh,            __gmpz_fdiv_q)