Re-implement `testPrimeInteger` predicate (#9281)
authorHerbert Valerio Riedel <hvr@gnu.org>
Fri, 28 Nov 2014 13:59:50 +0000 (14:59 +0100)
committerHerbert Valerio Riedel <hvr@gnu.org>
Fri, 28 Nov 2014 14:05:44 +0000 (15:05 +0100)
This also adds `testPrimeWord#` and `testPrimeBigNat` predicates.

`testPrimeInteger` has been available since `integer-gmp-0.5.1`
(added via f49735486533842cc84df70cafc8d565dffd75db).
The `nextPrimeInteger` function is still missing though.

libraries/integer-gmp2/cbits/wrappers.c
libraries/integer-gmp2/src/GHC/Integer/GMP/Internals.hs

index 9aac390..6c188a3 100644 (file)
@@ -449,3 +449,28 @@ integer_gmp_rscan_nzbyte(const uint8_t *srcptr,
 
   return 0;
 }
+
+/* wrapper around mpz_probab_prime_p */
+HsInt
+integer_gmp_test_prime(const mp_limb_t s[], const mp_size_t sn, const HsInt rep)
+{
+  if (!sn) return 0;
+
+  const mpz_t sz = {{
+      ._mp_alloc = sn,
+      ._mp_size  = sn,
+      ._mp_d = (mp_limb_t*)s
+    }};
+
+  // int mpz_probab_prime_p (const mpz_t n, int reps)
+  return mpz_probab_prime_p(sz, rep);
+}
+
+/* wrapper around mpz_probab_prime_p */
+HsInt
+integer_gmp_test_prime1(const mp_limb_t limb, const HsInt rep)
+{
+  if (!limb) return 0;
+
+  return integer_gmp_test_prime(&limb, 1, rep);
+}
index 77d73bf..480866b 100644 (file)
@@ -120,6 +120,11 @@ module GHC.Integer.GMP.Internals
     , gcdInt
     , gcdWord
 
+      -- * Primality tests
+    , testPrimeInteger
+    , testPrimeBigNat
+    , testPrimeWord#
+
       -- * Import/export functions
       -- ** Compute size of serialisation
     , sizeInBaseBigNat
@@ -280,3 +285,41 @@ exportWordToMutableByteArray (W# w#) = c_mpn_export1ToMutableByteArray# w#
 foreign import ccall unsafe "integer_gmp_mpn_export1"
   c_mpn_export1ToMutableByteArray# :: GmpLimb# -> MutableByteArray# RealWorld
                                    -> Word# -> Int# -> IO Word
+
+
+-- | Probalistic Miller-Rabin primality test.
+--
+-- \"@'testPrimeInteger' /n/ /k/@\" determines whether @/n/@ is prime
+-- and returns one of the following results:
+--
+-- * @2#@ is returned if @/n/@ is definitely prime,
+--
+-- * @1#@ if @/n/@ is a /probable prime/, or
+--
+-- * @0#@ if @/n/@ is definitely not a prime.
+--
+-- The @/k/@ argument controls how many test rounds are performed for
+-- determining a /probable prime/. For more details, see
+-- <http://gmplib.org/manual/Number-Theoretic-Functions.html#index-mpz_005fprobab_005fprime_005fp-360 GMP documentation for `mpz_probab_prime_p()`>.
+--
+-- /Since: 0.5.1.0/
+{-# NOINLINE testPrimeInteger #-}
+testPrimeInteger :: Integer -> Int# -> Int#
+testPrimeInteger (S# i#) = testPrimeWord# (int2Word# (absI# i#))
+testPrimeInteger (Jp# n) = testPrimeBigNat n
+testPrimeInteger (Jn# n) = testPrimeBigNat n
+
+-- | Version of 'testPrimeInteger' operating on 'BigNat's
+--
+-- /Since: 1.0.0.0/
+testPrimeBigNat :: BigNat -> Int# -> Int#
+testPrimeBigNat bn@(BN# ba#) = c_integer_gmp_test_prime# ba# (sizeofBigNat# bn)
+
+foreign import ccall unsafe "integer_gmp_test_prime"
+  c_integer_gmp_test_prime# :: ByteArray# -> GmpSize# -> Int# -> Int#
+
+-- | Version of 'testPrimeInteger' operating on 'Word#'s
+--
+-- /Since: 1.0.0.0/
+foreign import ccall unsafe "integer_gmp_test_prime1"
+  testPrimeWord# :: GmpLimb# -> Int# -> Int#