/* File: lssophie.c */

/* Purpose: Sophie Germain prime generator */

/* Standard headers */
#include <stdio.h>
#include <gmp.h>

#include <sys/mman.h>
#include <sys/resource.h>

static mpz_t two;
static mpz_t m, p;
static mpz_t t, u, v;

static inline int mpz_figurate_p(const mpz_t n, const mpz_t order)
{
	mpz_mul(t, order, order);
	mpz_sub(u, t, order);
	mpz_divexact_ui(t, u, 2);

	mpz_sub_ui(u, n, 1);
	mpz_add(v, u, n);

	mpz_powm(m, two, u, t);
	mpz_powm(p, two, v, t);

	mpz_sub(u, p, m);

	return mpz_congruent_p(u, order, t);
}

int main(void)
{
	static mpz_t sophie[2];

	mpz_init_set_ui(sophie[0], 11);

	mpz_init(sophie[1]);

	mpz_init_set_ui(two, 2);

	mpz_init(m);
	mpz_init(p);

	mpz_init(t);
	mpz_init(u);
	mpz_init(v);

	static int alt = 0;

	mlockall(MCL_CURRENT);

	while (1)
	{
		mpz_set_ui(sophie[!alt], 1);
		mpz_addmul_ui(sophie[!alt], sophie[alt], 2);

		if (mpz_figurate_p(sophie[alt], sophie[!alt]))
		{
			static struct rusage ru;

			getrusage(RUSAGE_SELF, &ru);

			mpz_out_str(stdout, 10, sophie[alt]);
			printf("\t \t%lu.%lus\n", ru.ru_utime.tv_sec, ru.ru_utime.tv_usec / 100000);

			alt = !alt;
			continue;
		}

		while (1)
		{
			mpz_add_ui(sophie[alt], sophie[alt], 12);

			if (mpz_gcd_ui(NULL, sophie[alt], 1078282205) == 1) break;
		}
	}
}

/* End: lssophie.c */