Commit c3315921 authored by David Eichmann's avatar David Eichmann 🏋 Committed by Ben Gamari

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
parent 2693eb11
......@@ -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);
......
......@@ -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] ]
......
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]
......
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment