GMP function to generate prime number bigger than X

Hans Petter Selasky hps at selasky.org
Wed Jan 26 14:41:54 CET 2022


On 1/26/22 14:17, Torbjörn Granlund wrote:
> Hans Petter Selasky <hps at selasky.org> writes:
> 
>    I haven't tried. Do you know if it is fast?
> 
> Why don't you try it?
> 
>    Basically I wanted to avoid the Miller Rabin test, because any prime
>    will do, not just the next one.
> 
> That statement makes no sense.

Hi Torbjörn,

Let me restate that. I need a prime, which also has the property, that 
two in the power of "N", does not become equal to one before "N" reaches 
"prime" minus one. Not all primes have the property. Not sure what this 
is called in formal language.

My program uses 6 seconds roughly for the same primes and sizes.

mpz_nextprime() uses/adds 2 seconds roughly for the same primes and sizes.

I didn't look into any optimizations, nor how this scales.

clang++ -O2 -I/usr/local/include -L/usr/local/lib -lgmp -lgmpxx -lm 
test.cpp && time ./a.out

--HPS

cat << EOF > test.cpp
#include <stdint.h>
#include <string.h>
#include <stdbool.h>
#include <math.h>
#include <iostream>
#include <gmpxx.h>

typedef mpz_class v_t;

static uint64_t
mbin_factor_slow_64(uint64_t x)
{
	uint64_t limit;
	uint64_t y;

	if (x <= 1)
		return (0);
	if (!(x & 1))
		return (2);

	limit = ((uint64_t)sqrt(x) + 4) | 1;

	if (limit > x)
		limit = x;
	if (limit < 3)
		limit = 3;

	for (y = 3; y != limit; y += 2) {
		if ((x % y) == 0)
			return (y);
	}
	return (0);
}

int
main()
{
	v_t k = 2;

	for (int a = 3; a < 16000; a++) {
		if (mbin_factor_slow_64(a) == 0)
			continue;
#if 1
		v_t r;
		v_t b = 2;
		v_t e = a * k;
		v_t m = a * k + 1;

		mpz_powm(r.get_mpz_t(), b.get_mpz_t(), e.get_mpz_t(), m.get_mpz_t());
		if (r == 1) {
			e /= 2;
			mpz_powm(r.get_mpz_t(), b.get_mpz_t(), e.get_mpz_t(), m.get_mpz_t());
			if (r == 1)
				continue;
			k *= a;
			std::cout << "factor " << (k + 1) << "\n";
#if 1
			v_t t = k + 199;
			v_t tr;

			mpz_nextprime(tr.get_mpz_t(), t.get_mpz_t());

			std::cout << "factor " << tr << "\n";
#endif
		}
#endif
	}
	return (0);
}
EOF


More information about the gmp-discuss mailing list