On the improvement of the the primality test function
Adrien Prost-Boucle
adrien.prost-boucle at laposte.net
Thu Oct 13 20:41:23 UTC 2016
... The 3 attached files have been wiped out, it seems.
I realize now that I should have posted on -devel mailing list?
Also C99 is required to compile my test code.
===================================================
Here is the patch for file src/mpz/pprime_p.c
===================================================
--- mpz/pprime_p.c.orig 2016-06-18 22:00:00.000000000 +0200
+++ mpz/pprime_p.c 2016-10-12 21:35:38.561167755 +0200
@@ -41,6 +41,162 @@
static int isprime (unsigned long int);
+
+// For 64-bit systems where the compiler has no support for 128-bit
+// A utility function to do multiply-modulo with 64-bit operands
+#if 0
+uint64_t mulmod64(uint64_t a, uint64_t b, uint64_t mod) {
+ uint64_t acc;
+
+ // Compute a * b into a 2-words accumulator
+ uint64_t acc_h, acc_l;
+ acc_l = (a & 0xFFFFFFFF) * (b & 0xFFFFFFFF);
+ acc_h = (a >> 32) * (b >> 32);
+ acc = (a & 0xFFFFFFFF) * (b >> 32);
+ acc_h += acc >> 32;
+ acc_l += acc << 32;
+ if(acc_l < acc << 32) acc_h++;
+ acc = (a >> 32) * (b & 0xFFFFFFFF);
+ acc_h += acc >> 32;
+ acc_l += acc << 32;
+ if(acc_l < acc << 32) acc_h++;
+
+ // Early exit condition
+ if(acc_h==0) return acc_l % mod;
+
+ // Work on a shifted version of the modulo, on 2 words
+ uint64_t mod_h = 0, mod_l = mod;
+ uint64_t mask_h = ((uint64_t)1) << 63;
+ // Shift the 2-words modulo until reaching the size of a * b
+ while(mod_h < acc_h && mod_h < mask_h) {
+ mod_h = (mod_h << 1) | (mod_l >> 63);
+ mod_l <<= 1;
+ }
+
+ // Iteratively subtract the 2-words modulo until reaching only one word size
+ do {
+ // Subtract the modulo if possible
+ if(acc_h > mod_h || (acc_h == mod_h && acc_l >= mod_l)) {
+ uint64_t prev_l = acc_l;
+ acc_l -= mod_l;
+ if(acc_l > prev_l) acc_h--;
+ acc_h -= mod_h;
+ }
+ // Shift the modulo
+ mod_l = (mod_h << 63) | (mod_l >> 1);
+ mod_h >>= 1;
+ } while(acc_h != 0);
+
+ return acc_l % mod;
+}
+#endif
+
+// Use umul_ppmm() macro and friends
+// FIXME Unsure about operand size for these macros
+#if 1
+static inline uint64_t mulmod64(uint64_t a, uint64_t b, uint64_t mod) {
+ // Compute a * b into a 2-words accumulator
+ uint64_t acc_h, acc_l;
+ umul_ppmm(acc_h, acc_l, a, b);
+ // Compute the remainder
+ uint64_t res_q, res_r;
+ udiv_qrnnd(res_q, res_r, acc_h, acc_l, mod);
+ return res_r;
+}
+#endif
+
+// Fast exponentiation modulo a given number
+static uint32_t powmod32(uint32_t val, uint32_t p, uint32_t mod) {
+ uint64_t accpow = val % mod;
+ uint64_t res = 1;
+ for( ; p>0; p>>=1) {
+ if(p % 2 != 0) res = (res * accpow) % mod;
+ accpow = (accpow * accpow) % mod;
+ }
+ return res;
+}
+static uint64_t powmod64(uint64_t val, uint64_t p, uint64_t mod) {
+
+ // FIXME Check if the compiler has support for 128-bit integers
+ #if 0
+
+ unsigned __int128 accpow = val % mod;
+ unsigned __int128 res = 1;
+ for( ; p>0; p>>=1) {
+ if(p % 2 != 0) res = (res * accpow) % mod;
+ accpow = (accpow * accpow) % mod;
+ }
+
+ #else
+
+ uint64_t accpow = val % mod;
+ uint64_t res = 1;
+ for( ; p>0; p>>=1) {
+ if(p % 2 != 0) res = mulmod64(res, accpow, mod);
+ accpow = mulmod64(accpow, accpow, mod);
+ }
+
+ #endif
+
+ return res;
+}
+
+/* Test: Strong probable-prime base a (a-SPRP), for n > 1
+ *
+ * See http://primes.utm.edu/prove/prove2_3.html
+ *
+ * Citing:
+ * Write n-1 = 2^s * d where d is odd and s is non-negative:
+ * n is a strong probable-prime base a (an a-SPRP) if either a^d = 1 (mod n)
+ * or (a^d)^(2^r) = -1 (mod n) for some non-negative r less than s.
+ *
+ * Note: this code requires n > 1 and n odd
+ */
+static int IsPrime32_TestSPRP(uint32_t N, uint32_t a) {
+ uint32_t d = N - 1;
+ unsigned s = 0;
+ while(d % 2 == 0) { s++; d/=2; }
+
+ // Check if a^d = 1 (mod n)
+ a = powmod32(a, d, N);
+ if(a == 1 || a == N-1) return 1;
+
+ // Check if (a^d)^(2^r) = -1 (mod n) for some non-negative r less than s
+ while(--s != 0) {
+ a = ((uint64_t)a * a) % N;
+ if(a == N-1) return 1;
+ }
+
+ return 0;
+}
+static int IsPrime64_TestSPRP(uint64_t N, uint64_t a) {
+ // Fallback to 32-bit code if possible
+ if(N >> 32 == 0 && a >> 32 == 0) return IsPrime32_TestSPRP(N, a);
+
+ uint64_t d = N - 1;
+ unsigned s = 0;
+ while(d % 2 == 0) { s++; d/=2; }
+
+ // Check if a^d = 1 (mod n)
+ a = powmod64(a, d, N);
+ if(a == 1 || a == N-1) return 1;
+
+ // Check if (a^d)^(2^r) = -1 (mod n) for some non-negative r less than s
+ while(--s != 0) {
+ // FIXME If 128-bit capable compiler, remove call to mulmod64()
+ a = mulmod64(a, a, N);
+ if(a == N-1) return 1;
+ }
+
+ return 0;
+}
+
+// FIXME Just to make that test file compile
+extern int rand();
+extern int srand(int seed);
+extern int time(void* p);
+
+
/* MPN_MOD_OR_MODEXACT_1_ODD can be used instead of mpn_mod_1 for the trial
division. It gives a result which is not the actual remainder r but a
value congruent to r*2^n mod d. Since all the primes being tested are
@@ -52,24 +208,112 @@
mp_limb_t r;
mpz_t n2;
- /* Handle small and negative n. */
- if (mpz_cmp_ui (n, 1000000L) <= 0)
- {
- int is_prime;
- if (mpz_cmpabs_ui (n, 1000000L) <= 0)
- {
- is_prime = isprime (mpz_get_ui (n));
- return is_prime ? 2 : 0;
- }
- /* Negative number. Negate and fall out. */
- PTR(n2) = PTR(n);
- SIZ(n2) = -SIZ(n);
- n = n2;
- }
+ // Negative number: negate and fall out
+ if(mpz_sgn(n) < 0) {
+ PTR(n2) = PTR(n);
+ SIZ(n2) = -SIZ(n);
+ n = n2;
+ }
+
+ // Separately handle numbers that fit a machine word
+ if(mpz_fits_ulong_p(n) != 0) {
+ unsigned long val = mpz_get_ui(n);
+ int isprime;
+
+ // Use pre-flagged prime numbers for all primes under 63
+ const unsigned long primeflags =
+ (1ULL << 2) | (1ULL << 3) | (1ULL << 5) | (1ULL << 7) |
+ (1ULL << 11) | (1ULL << 13) | (1ULL << 17) | (1ULL << 19) |
+ (1ULL << 23) | (1ULL << 29) |
+ (1ULL << 31) | (1ULL << 37) |
+ (1ULL << 41) | (1ULL << 43) | (1ULL << 47) |
+ (1ULL << 53) | (1ULL << 59) |
+ (1ULL << 61);
+
+ // Handle all numbers lower than 64
+ if(val <= 63) {
+ if(((primeflags >> val) & 1) != 0) return 2;
+ return 0;
+ }
+
+ // From here, even values are not prime
+ if(val % 2 == 0) return 0;
+
+ // Try division by small primes
+ if(
+ (val % 3 == 0) || (val % 5 == 0) || (val % 7 == 0) |
+ (val % 11 == 0) || (val % 13 == 0) || (val % 17 == 0) || (val % 19 == 0) |
+ (val % 23 == 0) || (val % 29 == 0) ||
+ (val % 31 == 0) || (val % 37 == 0) ||
+ (val % 41 == 0) || (val % 43 == 0) || (val % 47 == 0) ||
+ (val % 53 == 0) || (val % 59 == 0) ||
+ (val % 61 == 0)
+ ) return 0;
+
+ // For small values, check all odd numbers up to the square root
+ // Begin at first unchecked prime
+ if(val <= 100000) {
+ unsigned r=1, d, q = val;
+ for(d=67; r!=0; d+=2) {
+ q = val / d;
+ r = val - d * q;
+ if(d > q) return 2;
+ }
+ return 0;
+ }
+
+ // Launch strong probable-prime tests
+ isprime = IsPrime64_TestSPRP(val, 2);
+ if(isprime==0) return 0;
+ isprime = IsPrime64_TestSPRP(val, 3);
+ if(isprime==0) return 0;
+
+ if(val < 1373653) return 2;
+
+ isprime = IsPrime64_TestSPRP(val, 5);
+ if(isprime==0) return 0;
+
+ if(val < 25326001) return 2;
+
+ if(val == 3215031751) return 0;
+
+ isprime = IsPrime64_TestSPRP(val, 7);
+ if(isprime==0) return 0;
+
+ if(val < 118670087467ULL) return 2;
+
+ isprime = IsPrime64_TestSPRP(val, 11);
+ if(isprime==0) return 0;
+
+ if(val < 2152302898747ULL) return 2;
+
+ isprime = IsPrime64_TestSPRP(val, 13);
+ if(isprime==0) return 0;
+
+ if(val < 3474749660383ULL) return 2;
+
+ isprime = IsPrime64_TestSPRP(val, 17);
+ if(isprime==0) return 0;
+
+ if(val < 341550071728321ULL) return 2;
+
+ // Finally, apply extra random strong probable-prime tests
+ // Follow algorithm from mpz_millerrabin()
+ srand(time(NULL));
+ for(unsigned r=7; r<reps; r++) {
+ // Randomly generate a number between 19 and n
+ unsigned long a = ((unsigned long)rand()) % (val - 19) + 19;
+ // Perform the test
+ isprime = IsPrime64_TestSPRP(val, a);
+ if(isprime==0) return 0;
+ }
+
+ // Probably prime
+ return 1;
+ } // End of handling when the numbers that fit a machine word
/* If n is now even, it is not a prime. */
- if ((mpz_get_ui (n) & 1) == 0)
- return 0;
+ if(mpz_even_p(n) != 0) return 0;
#if defined (PP)
/* Check if n has small factors. */
===================================================
===================================================
Here is the benchmarking program primes.c
===================================================
/*
Little program that deals with prime numbers
*/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <stdint.h>
#include <inttypes.h>
#include <stdbool.h>
#include <math.h>
#include <limits.h>
#include <time.h>
#include <gmp.h>
//==================================================
// Utility functions
//==================================================
// A utility function
double TimeSpec_ToDouble(struct timespec *ts) {
return (double)ts->tv_sec + ((double)ts->tv_nsec)/1.0e9;
}
double Time_GetReal() {
struct timespec newtime;
clock_gettime(CLOCK_REALTIME, &newtime);
return TimeSpec_ToDouble(&newtime);
}
// For 64-bit systems where the compiler has no support for 128-bit
// A utility function to do multiply-modulo with 64-bit operands
#if 0
uint64_t mulmod64(uint64_t a, uint64_t b, uint64_t mod) {
uint64_t acc;
// Compute a * b into a 2-words accumulator
uint64_t acc_h, acc_l;
acc_l = (a & 0xFFFFFFFF) * (b & 0xFFFFFFFF);
acc_h = (a >> 32) * (b >> 32);
acc = (a & 0xFFFFFFFF) * (b >> 32);
acc_h += acc >> 32;
acc_l += acc << 32;
if(acc_l < acc << 32) acc_h++;
acc = (a >> 32) * (b & 0xFFFFFFFF);
acc_h += acc >> 32;
acc_l += acc << 32;
if(acc_l < acc << 32) acc_h++;
// Early exit condition
if(acc_h==0) return acc_l % mod;
// Work on a shifted version of the modulo, on 2 words
uint64_t mod_h = 0, mod_l = mod;
uint64_t mask_h = ((uint64_t)1) << 63;
// Shift the 2-words modulo until reaching the size of a * b
while(mod_h < acc_h && mod_h < mask_h) {
mod_h = (mod_h << 1) | (mod_l >> 63);
mod_l <<= 1;
}
// Iteratively subtract the 2-words modulo until reaching only
one word size
do {
// Subtract the modulo if possible
if(acc_h > mod_h || (acc_h == mod_h && acc_l >= mod_l))
{
uint64_t prev_l = acc_l;
acc_l -= mod_l;
if(acc_l > prev_l) acc_h--;
acc_h -= mod_h;
}
// Shift the modulo
mod_l = (mod_h << 63) | (mod_l >> 1);
mod_h >>= 1;
} while(acc_h != 0);
return acc_l % mod;
}
#endif
#if 1
// From GMP
typedef long long int DItype;
typedef unsigned long long int UDItype;
// Fom GMP with: defined (__amd64__) && W_TYPE_SIZE == 64
#define umul_ppmm(w1, w0, u, v) \
__asm__ ("mulq %3"
\
: "=a" (w0), "=d" (w1)
\
: "%0" ((UDItype)(u)), "rm" ((UDItype)(v)))
#define udiv_qrnnd(q, r, n1, n0, dx) /* d renamed to dx avoiding "=d"
*/\
__asm__ ("divq %4" /* stringification in K&R C */
\
: "=a" (q), "=d" (r)
\
: "0" ((UDItype)(n0)), "1" ((UDItype)(n1)), "rm"
((UDItype)(dx)))
static inline uint64_t mulmod64(uint64_t a, uint64_t b, uint64_t mod) {
// Compute a * b into a 2-words accumulator
uint64_t acc_h, acc_l;
umul_ppmm(acc_h, acc_l, a, b);
// Compute the remainder
uint64_t res_q, res_r;
udiv_qrnnd(res_q, res_r, acc_h, acc_l, mod);
return res_r;
}
#endif
// Fast exponentiation modulo a given number
uint32_t powmod32(uint32_t val, uint32_t p, uint32_t mod) {
uint64_t accpow = val % mod;
uint64_t res = 1;
for( ; p>0; p>>=1) {
if(p % 2 != 0) res = (res * accpow) % mod;
accpow = (accpow * accpow) % mod;
}
return res;
}
uint64_t powmod64(uint64_t val, uint64_t p, uint64_t mod) {
uint64_t accpow = val % mod;
uint64_t res = 1;
for( ; p>0; p>>=1) {
if(p % 2 != 0) res = mulmod64(res, accpow, mod);
accpow = mulmod64(accpow, accpow, mod);
}
return res;
}
// Factorial modulo a given number, naive: complexity O(N)
uint32_t factmod32_naive(uint32_t val, uint32_t mod) {
uint64_t res = 1;
for(uint32_t i=2; i<=val; i++) res = (res * i) % mod;
return res;
}
//==================================================
// Check whether a number is prime
//==================================================
/* Test: Strong probable-prime base a (a-SPRP), for n > 1
*
* See http://primes.utm.edu/prove/prove2_3.html
*
* Citing:
* Write n-1 = 2^s * d where d is odd and s is non-negative:
* n is a strong probable-prime base a (an a-SPRP) if either a^d = 1
(mod n)
* or (a^d)^(2^r) = -1 (mod n) for some non-negative r less than s.
*
* Note: this code requires n > 1 and n odd
*/
bool IsPrime32_TestSPRP(uint32_t N, uint32_t a) {
uint32_t d = N - 1;
unsigned s = 0;
while(d % 2 == 0) { s++; d/=2; }
// Check if a^d = 1 (mod n)
a = powmod32(a, d, N);
if(a == 1 || a == N-1) return true;
// Check if (a^d)^(2^r) = -1 (mod n) for some non-negative r
less than s
while(--s != 0) {
a = ((uint64_t)a * a) % N;
if(a == N-1) return true;
}
return false;
}
bool IsPrime64_TestSPRP(uint64_t N, uint64_t a) {
// Fallback to 32-bit code if possible
if(N >> 32 == 0 && a >> 32 == 0) return IsPrime32_TestSPRP(N,
a);
uint64_t d = N - 1;
unsigned s = 0;
while(d % 2 == 0) { s++; d/=2; }
// Check if a^d = 1 (mod n)
a = powmod64(a, d, N);
if(a == 1 || a == N-1) return true;
// Check if (a^d)^(2^r) = -1 (mod n) for some non-negative r
less than s
while(--s != 0) {
a = mulmod64(a, a, N);
if(a == N-1) return true;
}
return false;
}
//==================================================
// Primality test - For comparison with GMP lib
//==================================================
/* Making GMP primatlity test faster for small numbers, and returning
proven primality more often
* See http://primes.utm.edu/prove/prove2_3.html
*
* FIXME This is 2-10 times faster than GMP on numbers within 32-bits
range, above 1000000
* FIXME Test in full 64-bits range because Miller-Rabin test uses 64-
bit implementation and should be faster
*/
// GMP Field access macros
#define SIZ(x) ((x)->_mp_size)
#define PTR(x) ((x)->_mp_d)
#define numberof(x) (sizeof (x) / sizeof ((x)[0]))
// Return 0 if not prime, 1 if probably prime, 2 if certainly prime
int IsPrimeGMP_Probable(mpz_t n, int reps) {
mpz_t n2;
// Negative number: negate and fall out
if(mpz_sgn(n) < 0) {
PTR(n2) = PTR(n);
SIZ(n2) = -SIZ(n);
n = n2;
}
// Separately handle numbers that fit a machine word
if(mpz_fits_ulong_p(n) != 0) {
unsigned long val = mpz_get_ui(n);
int isprime;
// Use pre-flagged prime numbers for all primes under
63
const unsigned long primeflags =
(1ULL << 2) | (1ULL << 3) | (1ULL << 5) |
(1ULL << 7) |
(1ULL << 11) | (1ULL << 13) | (1ULL << 17) |
(1ULL << 19) |
(1ULL << 23) | (1ULL << 29) |
(1ULL << 31) | (1ULL << 37) |
(1ULL << 41) | (1ULL << 43) | (1ULL << 47) |
(1ULL << 53) | (1ULL << 59) |
(1ULL << 61);
// Handle all numbers lower than 64
if(val <= 63) {
if(((primeflags >> val) & 1) != 0) return 2;
return 0;
}
// From here, even values are not prime
if(val % 2 == 0) return 0;
// Try division by small primes
if(
(val % 3 == 0) || (val % 5 == 0) || (val % 7
== 0) |
(val % 11 == 0) || (val % 13 == 0) || (val % 17
== 0) || (val % 19 == 0) |
(val % 23 == 0) || (val % 29 == 0) ||
(val % 31 == 0) || (val % 37 == 0) ||
(val % 41 == 0) || (val % 43 == 0) || (val % 47
== 0) ||
(val % 53 == 0) || (val % 59 == 0) ||
(val % 61 == 0)
) return 0;
// Mimicking GMP: For small values, check all odd
numbers up to the square root
// Begin at first unchecked prime
if(val <= 100000) {
unsigned r=1, d, q = val;
for(d=67; r!=0; d+=2) {
// FIXME There may be a macro or an ASM
instruction to do div and rem simultaneously
q = val / d;
r = val - d * q;
if(d > q) return 2;
}
return 0;
}
// Launch strong probable-prime tests
isprime = IsPrime64_TestSPRP(val, 2);
if(isprime==false) return 0;
isprime = IsPrime64_TestSPRP(val, 3);
if(isprime==false) return 0;
if(val < 1373653) return 2;
isprime = IsPrime64_TestSPRP(val, 5);
if(isprime==false) return 0;
if(val < 25326001) return 2;
if(val == 3215031751) return 0;
isprime = IsPrime64_TestSPRP(val, 7);
if(isprime==false) return 0;
if(val < 118670087467ULL) return 2;
isprime = IsPrime64_TestSPRP(val, 11);
if(isprime==false) return 0;
if(val < 2152302898747ULL) return 2;
isprime = IsPrime64_TestSPRP(val, 13);
if(isprime==false) return 0;
if(val < 3474749660383ULL) return 2;
isprime = IsPrime64_TestSPRP(val, 17);
if(isprime==false) return 0;
if(val < 341550071728321ULL) return 2;
// Finally, apply extra random strong probable-prime
tests
// Follow algorithm from mpz_millerrabin()
for(unsigned r=7; r<reps; r++) {
// Randomly generate a number between 19 and n
unsigned long a = ((unsigned long)rand()) %
(val - 19) + 19;
// Perform the test
isprime = IsPrime64_TestSPRP(val, a);
if(isprime==false) return 0;
}
// Probably prime
return 1;
} // End of handling when the numbers that fit a machine word
// If n is even, it is not a prime
if(mpz_even_p(n) != 0) return 0;
// FIXME Here apply preliminary tests done in GMP: Test
divisibility by small primes
// From here we launch Miller-Rabbin exactly like GMP
return mpz_millerrabin(n, reps);
}
//==================================================
// Main function
//==================================================
void print_help(int argc, char **argv) {
printf("\n");
printf("Syntax: %s [command]\n", argv[0]);
printf("\n");
printf("Options:\n");
printf(" -u, -h, --help Print this help and exit\n");
printf("\n");
printf("Commands:\n");
printf(" isprime32-SPRP <n> <a>\n");
printf(" isprime64-SPRP <n> <a>\n");
printf(" isprimeGMP-prob <n> <reps>\n");
printf(" isprimeGMP-lib <n> <reps>\n");
printf(" isprimeGMP-check <limit>\n");
printf(" isprimeGMP-bench <limit>\n");
printf("\n");
}
int main(int argc, char **argv){
if(argc==1) {
print_help(argc, argv);
exit(EXIT_SUCCESS);
}
int argi = 1;
for( ; argi<argc; argi++) {
char* arg = argv[argi];
if(arg[0]!='-') break;
if(strcmp(arg, "-u")==0 || strcmp(arg, "-h")==0 ||
strcmp(arg, "--help")==0) {
print_help(argc, argv);
exit(EXIT_SUCCESS);
}
else {
printf("Error: Unknown option '%s'\n", arg);
exit(EXIT_FAILURE);
}
}
// If there is no command, exit silently
if(argc - argi < 1) {
exit(EXIT_SUCCESS);
}
srand(time(NULL));
// Get the main command
char* cmd = argv[argi++];
// Utility functions
char* getparam_str() {
if(argi>=argc) {
printf("ERROR : missing parameters for command
'%s'\n", cmd);
exit(EXIT_FAILURE);
}
return argv[argi++];
}
unsigned getparam_ui() {
return strtoul(getparam_str(), NULL, 0);
}
long long getparam_ull() {
return strtoul(getparam_str(), NULL, 0);
}
void isprime_result(bool isprime) {
if(isprime==true) printf("Prime\n");
else printf("Non-prime\n");
exit(EXIT_SUCCESS);
}
void isprime_result_prob(int isprime) {
if(isprime==0) printf("Definitely non-prime\n");
else if(isprime==2) printf("Definitely prime\n");
else printf("Probably prime\n");
exit(EXIT_SUCCESS);
}
// Process the main command
if(strcmp(cmd, "isprime32-SPRP")==0) {
unsigned N = getparam_ui();
unsigned a = getparam_ui();
if(N <= 1 || N % 2 == 0) {
printf("Error: N must be odd and >1\n");
exit(EXIT_FAILURE);
}
bool isprime = IsPrime32_TestSPRP(N, a);
isprime_result(isprime);
}
else if(strcmp(cmd, "isprime64-SPRP")==0) {
uint64_t N = getparam_ull();
uint64_t a = getparam_ui();
if(N <= 1 || N % 2 == 0) {
printf("Error: N must be odd and >1\n");
exit(EXIT_FAILURE);
}
bool isprime = IsPrime64_TestSPRP(N, a);
isprime_result(isprime);
}
else if(strcmp(cmd, "isprimeGMP-prob")==0) {
char* param = getparam_str();
unsigned reps = getparam_ui();
mpz_t N;
mpz_init_set_str(N, param, 10);
int isprime = IsPrimeGMP_Probable(N, reps);
mpz_clear(N);
isprime_result_prob(isprime);
}
else if(strcmp(cmd, "isprimeGMP-lib")==0) {
char* param = getparam_str();
unsigned reps = getparam_ui();
mpz_t N;
mpz_init_set_str(N, param, 10);
int isprime = mpz_probab_prime_p(N, reps);
mpz_clear(N);
isprime_result_prob(isprime);
}
else if(strcmp(cmd, "isprimeGMP-check")==0) {
unsigned limit = getparam_ui();
const unsigned reps = 25;
mpz_t N;
mpz_init(N);
unsigned superior_nb = 0;
unsigned inferior_nb = 0;
for(unsigned long val=1; val<limit; val++) {
mpz_set_ui(N, val);
//#define SETVAL(val) val
//#define SETVAL(val) (100000000 + val)
#define SETVAL(val) (100000000 + rand())
//#define SETVAL(val) rand()
//#define SETVAL(val) ((((uint64_t)rand()) <<
32) + rand())
mpz_set_ui(N, SETVAL(val));
#undef SETVAL
int isprime_loc = IsPrimeGMP_Probable(N, reps);
int isprime_lib = mpz_probab_prime_p(N, reps);
if( (isprime_loc > 0) != (isprime_lib > 0) ) {
printf("Error for %lu: loc = %u, lib =
%u\n", val, isprime_loc, isprime_lib);
}
if(isprime_loc == 2 && isprime_lib == 1)
superior_nb++;
if(isprime_loc == 1 && isprime_lib == 2)
inferior_nb++;
}
mpz_clear(N);
printf("Results: inferior %u, superior %u\n",
inferior_nb, superior_nb);
}
else if(strcmp(cmd, "isprimeGMP-bench")==0) {
unsigned limit = getparam_ui();
const unsigned reps = 25;
unsigned phony_result = 0;
double tbeg;
mpz_t N;
mpz_init(N);
//#define SETVAL(val) val
//#define SETVAL(val) (100000000 + val)
//#define SETVAL(val) (100000000 + rand())
//#define SETVAL(val) rand()
#define SETVAL(val) ((((uint64_t)rand()) << 32) +
rand())
srand(1);
tbeg = Time_GetReal();
for(unsigned long val=1; val<limit; val++) {
mpz_set_ui(N, SETVAL(val));
phony_result += IsPrimeGMP_Probable(N, reps);
}
printf("Local function .. %g s\n", Time_GetReal()-
tbeg);
srand(1);
tbeg = Time_GetReal();
for(unsigned long val=1; val<limit; val++) {
mpz_set_ui(N, SETVAL(val));
phony_result += mpz_probab_prime_p(N, reps);
}
printf("GMP library ..... %g s\n", Time_GetReal()-
tbeg);
#undef SETVAL
// Note: needed to prevent the compiler from removing
function calls
printf("Phony result .... %u\n", phony_result);
mpz_clear(N);
}
else {
printf("ERROR : Unknown command '%s'\n", cmd);
exit(EXIT_FAILURE);
}
exit(EXIT_SUCCESS);
}
===================================================
===================================================
Here is the Makefile to compile the program primes.c
===================================================
# Programs and parameters
CC = gcc
CFLAGS = -Wall -std=gnu99
LD = gcc
LDFLAGS = -lm -lgmp -lmpfr
ifdef DEBUG
CFLAGS += -g -O0
LDFLAGS +=
else
CFLAGS += -march=native -Ofast -DNDEBUG
LDFLAGS += -s
endif
RM = rm
RMFLAGS = -f
# Main target
PROG = primes
SRCS := primes.c
OBJS := $(patsubst %.c, %.o, $(SRCS))
# Link of the main program
$(PROG): $(OBJS)
$(LD) $(LDFLAGS) $(OBJS) -o $(PROG)
# Compilation of the source files
%.o: %.c
$(CC) $(CFLAGS) -c $<
# Delete the object files
clean:
$(RM) $(RMFLAGS) $(PROG) $(OBJS)
===================================================
Regards,
Adrien
On Wed, 2016-10-12 at 22:52 +0200, Adrien Prost-Boucle wrote:
> Hello,
>
> I have recently discovered that there exists a simple and fast way of
> deciding for sure whether some "small" numbers are primes. Small means
> under 341,550,071,728,321, which fits in 64 bits.
> This is according to the following web page and its references:
> http://primes.utm.edu/prove/prove2_3.html
>
> After observing the results of the GMP primality test
> function, mpz_probab_prime_p(), I noticed that the results are
> probabilistic above 1000000 and that the algorithm SPRP descibed in the
> page above are not present.
> So I decided to implement it and do some benchmarking :-)
>
> First, I implemented 32-bit and 64-bit versions of the algorithm so
> tests are done the fastest possible for numbers that fit in 32-bit and
> 64-bit machine registers.
>
> Then, I investigated what was the appropriate new limit under which
> testing divisibility by odd numbers can be done. I lowered the GMP
> limit of 1000000 by 100000. The 32-bit SPRP implementation is fast
> enough to take care of the rest with speedup, while returning proven
> primality result.
>
> And finally I use the 64-bit SPRP function for all numbers that fit
> into an unsigned long on my machine.
>
> The resulting code is attached, in file pprime_p.c.
> It can be used a drop-in replacement for the corresponding file in
> src/mpz/ of GMP (for test only).
> Note that it is _far_ from being ready to put in GMP, for reasons
> discussed below.
>
> I conducted tests to measure average speedup and accuracy of results
> compared to vanilla GMP.
>
> Attached is file primes.c and Makefile for compilation.
> It's for amd64-capable systems.
> Run the following command to compare against GMP:
> ./primes isprimeGMP-check 20000
> Run the following command to have execution time:
> ./primes isprimeGMP-bench 20000
> Modify the macro SETVAL in both commands to run tests on different
> ranges of numbers and with some randomness.
>
> My observations are the following:
>
> In the 32-bit range and up to the last number whose primality is
> proven, 341,550,071,728,321, the speedup varies between 2x and 10x
> (while returning exact primality, not probabilistic).
>
> Above that limit however, there is some slowdown, in the order of 30%-
> 40%. I am unsure where this comes from... I suspect my implementation
> of the powmod64 function is not optimal, the loop can run up to 64
> times due to shift, maybe there is a more optimized solution using some
> macros or ASM instructions?
> Or simply do a few more trial divisions to rapidly find small factors,
> like already done when numbers don't fit an unsigned long?
> Surely the slowdown on this range can be solved. Feedback about that
> would be much appreciated!
>
> My code in pprime_p.c is _far_ from being ready to put in GMP, in fact
> it is more of a preliminary implementation for test and as support for
> discussion about algorithms.
>
> My understanding of GMP internals is quite low.
> In particular, I'm often unsure about how to use utility macros such as umul_ppmm and friends, and about all macros that define on which kind of target we're compiling what should be checked about their capabilities.
> Advice about that is also most welcomed!
>
> Last point of this message, about the tests.
> During implementation of course I occasionally produced bugs which I had a hard time to corner. I have a few suggestions to improve tests related to primality.
>
> tests/mpz/t-pprime_p.c
> => I suggest adding primality test of all values between one known prime and the next prime, and ensure returned primality is correct (because we do know actual primality). Maybe in or or 2 such select ranges. This is linked to other suggestion below.
>
> tests/mpz/t-nextprime.c, function refmpz_nextprime()
> => Make it check if an abnormally high number of increments is done.
> Because this function, as well as the mpz_nextprime() under test, use mpz_probab_prime_p(). So if there is a problem in mpz_probab_prime_p() the current code may run indefinitely, if e.g. all numbers are considered non prime by mpz_probab_prime_p() above a certain number:
>
> void
> refmpz_nextprime (mpz_ptr p, mpz_srcptr t)
> {
> mpz_add_ui (p, t, 1L);
> unsigned c = 0;
> while (! mpz_probab_prime_p (p, 10))
> {
> mpz_add_ui (p, p, 1L);
> c ++;
> if(c > 10000)
> {
> gmp_printf ("Error: too many tested numbers at %Zu\n", p);
> abort();
> }
> }
> }
>
> tests/mpz/t-nextprime.c, function main():
> Make it print the values before abort(), when values are different.
> This eases investigation:
>
> if (mpz_cmp (nxtp, ref_nxtp) != 0)
> {
> gmp_printf ("From %Zu\n", x);
> gmp_printf ("got %Zu\n", nxtp);
> gmp_printf ("want %Zu\n", ref_nxtp);
> abort ();
> }
>
> Best regards,
> Adrien
> _______________________________________________
> gmp-discuss mailing list
> gmp-discuss at gmplib.org
> https://gmplib.org/mailman/listinfo/gmp-discuss
More information about the gmp-discuss
mailing list