Correct limb length and assertion for gcdExtInteger

Reviewers: hvr, bgamari, monoidal

Reviewed By: monoidal

Subscribers: monoidal, rwbarton, thomie, carter

GHC Trac Issues: #15350

Differential Revision: https://phabricator.haskell.org/D5042

3 files changed, 46 insertions, 3 deletions

diff --git a/libraries/integer-gmp/cbits/wrappers.c b/libraries/integer-gmp/cbits/wrappers.c
index 8f147ad0c1..11e5179323 100644
--- a/libraries/integer-gmp/cbits/wrappers.c
+++ b/libraries/integer-gmp/cbits/wrappers.c
@@ -286,9 +286,9 @@ integer_gmp_mpn_gcd(mp_limb_t r[],
  * reconstructed).
  *
  * g must have space for exactly gn=min(xn,yn) limbs.
- * s must have space for at least xn limbs.
+ * s must have space for at least yn limbs.
  *
- * return value: signed 'sn' of {sp,sn}
+ * return value: signed 'sn' of {sp,sn} where |sn| >= 1
  */
 mp_size_t
 integer_gmp_gcdext(mp_limb_t s0[], mp_limb_t g0[],
@@ -305,15 +305,25 @@ integer_gmp_gcdext(mp_limb_t s0[], mp_limb_t g0[],
 
   mpz_gcdext (g, s, NULL, x, y);
 
+  // g must be positive (0 <= gn).
+  // According to the docs for mpz_gcdext(), we have:
+  //       g < min(|y|/2|s|, |x|/2|t|)
+  // -->   g < min(|y|, |x|)
+  // -->   gn <= min(yn, xn)
+  // <->   gn <= gn0
   const mp_size_t gn = g[0]._mp_size;
   assert(0 <= gn && gn <= gn0);
   memset(g0, 0, gn0*sizeof(mp_limb_t));
   memcpy(g0, g[0]._mp_d, gn*sizeof(mp_limb_t));
   mpz_clear (g);
 
+  // According to the docs for mpz_gcdext(), we have:
+  //       |s| < |y| / 2g
+  // -->   |s| < |y|         (note g > 0)
+  // -->   sn <= yn
   const mp_size_t ssn = s[0]._mp_size;
   const mp_size_t sn  = mp_size_abs(ssn);
-  assert(sn <= mp_size_abs(xn));
+  assert(sn <= mp_size_abs(yn));
   memcpy(s0, s[0]._mp_d, sn*sizeof(mp_limb_t));
   mpz_clear (s);
 
diff --git a/testsuite/tests/lib/integer/integerGmpInternals.hs b/testsuite/tests/lib/integer/integerGmpInternals.hs
index c90df5c90d..e45c6f4a48 100644
--- a/testsuite/tests/lib/integer/integerGmpInternals.hs
+++ b/testsuite/tests/lib/integer/integerGmpInternals.hs
@@ -79,6 +79,8 @@ main = do
     print $ powModInteger b e m
     print $ powModInteger b e (m-1)
     print $ powModSecInteger b e (m-1)
+
+    putStrLn "\n# gcdExtInteger"
     print $ gcdExtInteger b e
     print $ gcdExtInteger e b
     print $ gcdExtInteger x y
@@ -86,10 +88,27 @@ main = do
     print $ gcdExtInteger x (-y)
     print $ gcdExtInteger (-x) y
     print $ gcdExtInteger (-x) (-y)
+
+    -- see Trac #15350
+    do
+        let a = 2
+            b = 2^65 + 1
+        print $ gcdExtInteger a b
+        print $ gcdExtInteger a (-b)
+        print $ gcdExtInteger (-a) b
+        print $ gcdExtInteger (-a) (-b)
+        print $ gcdExtInteger b a
+        print $ gcdExtInteger b (-a)
+        print $ gcdExtInteger (-b) a
+        print $ gcdExtInteger (-b) (-a)
+
+    putStrLn "\n# powInteger"
     print $ powInteger 12345 0
     print $ powInteger 12345 1
     print $ powInteger 12345 30
     print $ [ (x,i) | x <- [-7..71], let i = recipModInteger x (2*3*11*11*17*17), i /= 0 ]
+
+    putStrLn "\n# nextPrimeInteger"
     print $ I.nextPrimeInteger b
     print $ I.nextPrimeInteger e
     print $ [ k | k <- [ 0 .. 200 ], S# (I.testPrimeInteger k 25#) `elem` [1,2] ]
diff --git a/testsuite/tests/lib/integer/integerGmpInternals.stdout b/testsuite/tests/lib/integer/integerGmpInternals.stdout
index d5c1374e39..cff835bd27 100644
--- a/testsuite/tests/lib/integer/integerGmpInternals.stdout
+++ b/testsuite/tests/lib/integer/integerGmpInternals.stdout
@@ -1,6 +1,8 @@
 1527229998585248450016808958343740453059
 682382427572745901624116300491295556924
 682382427572745901624116300491295556924
+
+# gcdExtInteger
 (1,-238164827888328100873319793437342927637138278785737103723156342382925)
 (1,302679100340807588460107986194035692812415103244388831792688023418704)
 (92889294,115110207004456909698806038261)
@@ -8,10 +10,22 @@
 (92889294,115110207004456909698806038261)
 (92889294,-115110207004456909698806038261)
 (92889294,-115110207004456909698806038261)
+(1,-18446744073709551616)
+(1,-18446744073709551616)
+(1,18446744073709551616)
+(1,18446744073709551616)
+(1,1)
+(1,1)
+(1,-1)
+(1,-1)
+
+# powInteger
 1
 12345
 555562377826831043419246079513769804614412256811161773362797946971665712715296306339052301636736176350153982639312744140625
 [(-7,149867),(-5,167851),(-1,209813),(1,1),(5,41963),(7,59947),(13,177535),(19,143557),(23,182447),(25,134281),(29,7235),(31,33841),(35,95915),(37,113413),(41,61409),(43,24397),(47,174101),(49,158431),(53,193979),(59,188477),(61,185737),(65,35507),(67,118999),(71,186173)]
+
+# nextPrimeInteger
 2988348162058574136915891421498819466320163312926952423791023078876343
 2351399303373464486466122544523690094744975233415544072992656881240451
 [2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97,101,103,107,109,113,127,131,137,139,149,151,157,163,167,173,179,181,191,193,197,199]