Skip to content

Instantly share code, notes, and snippets.

@jrodatus jrodatus/prob_psq.c
Last active Aug 18, 2017

Embed
What would you like to do?
prob_psq
/*
* Demo of a close to log-time probabilistic test for squaredness, per:
* http://math.stackexchange.com/a/1360053/253804
*
* Requires libgmp, the GNU Multiple Precision Arithmetic Library:
* https://gmplib.org/
*
* To build, run:
* gcc -O2 prob_psq.c -lgmp -o prob_psq
*/
#include <gmp.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#define HR "========================================"\
"========================================"
/* generic POSIX errno's */
#define ENONE 0
#define ENOMEM 12 /* Cannot allocate memory */
#define EINVAL 22 /* Invalid argument */
#define EDOM 33 /* Numerical argument out of domain */
#define ERANGE 34 /* Result too large */
#define doerr(e) (printf("Error: %s\n", #e), e)
#define UL_MAX ((unsigned long)(-1)) /* usually 2^64-1 */
#define COL_SPC 3
#define DEFAULT_NP 30 /* default # primes, if not set by user */
typedef unsigned long ul;
static gmp_randstate_t rst;
static ul *primes = NULL;
static ul np = 0;
static ul nt = 0;
static mpz_t B;
int gen_primes() {
int err = ENONE;
mpz_t P;
ul i;
if (!np) {
/* set to something reasonable if uninitialized */
printf("# primes not set! Defaulting to %d\n", DEFAULT_NP);
np = DEFAULT_NP;
}
if (!(primes = reallocf(primes, np * sizeof(ul)))) {
err = doerr(ENOMEM);
goto out_err;
}
mpz_init_set_ui(P, 2); /* start at 3 */
printf("Generating %lu primes...\n", np);
for (i = 0; i < np; ++i) {
mpz_nextprime(P, P);
if (mpz_cmp_ui(P, UL_MAX) > 0) {
gmp_printf("\nprime too large for unsigned long: %Zd\n",
P);
err = doerr(ERANGE);
goto out_err;
}
primes[i] = mpz_get_ui(P);
printf(" %lu/%lu\r", i, np);
fflush(stdout);
}
printf("\nDone: P_max=%lu\n\n", primes[np-1]);
out_err:
return err;
}
int probab_perfect_square(mpz_t N, int *errp) {
int err = ENONE;
mpz_t P;
int i;
int prob_psq;
mpz_init(P);
if (!primes || !np) {
if ((err = gen_primes())) {
prob_psq = -1;
goto out_clear;
}
}
if (mpz_cmp_ui(N, primes[np-1]) <= 0) {
gmp_printf("N=%Zd is <= largest prime (%lu)\n", N, primes[np-1]);
err = doerr(EDOM);
prob_psq = -1;
goto out_clear;
}
prob_psq = 1; /* start by supposing N is a perfect square */
for (i = 0; i < np; ++i) {
mpz_set_ui(P, primes[i]);
if (mpz_legendre(N, P) < 0) {
/* if N is a quadratic non-residue (mod P),
* then it cannot possibly be a perfect square */
prob_psq = 0;
goto out_clear;
}
}
/* N passed all tests */
out_clear:
mpz_clear(P);
if (errp) {
*errp = err;
}
return prob_psq;
}
int set_conditions() {
int err = ENONE;
printf("Number of primes: ");
scanf("%lu", &np);
printf("Number of tests (T): ");
scanf("%lu", &nt);
printf("Upper bound for N: ");
gmp_scanf("%Zd", B);
putchar('\n');
if ((err = gen_primes())) {
goto out_err;
}
if (mpz_cmp_ui(B, primes[np-1]) <= 0) {
gmp_printf("Upper test bound (%Zd) is <= "
"largest prime (%lu)\n", B, primes[np-1]);
err = doerr(EINVAL);
goto out_err;
}
out_err:
return err;
}
int test_residue_stats() {
static const char *col[] = {"P", "fraction of N's that were residues"};
int err = ENONE;
mpz_t N, P;
ul *nqr = NULL;
ul i, j;
int col0_w;
int row_w;
mpz_init(N);
mpz_init(P);
if ((err = set_conditions())) {
goto out_free;
}
if (!(nqr = calloc(np, sizeof(ul)))) {
err = doerr(ENOMEM);
goto out_free;
}
puts(HR);
gmp_printf("We will test T=%lu random integers N, where:\n\n"
"\tP_max = %lu < N <= %Zd\n\n"
"Each row P shows the number of these N's that were "
"quadratic residues (mod P),\nwritten as a fraction of T.\n\n"
"If ~50%% of randomly-chosen N's were residues (mod P),\n"
"that suggests a probability of 50%% that a given N > P "
"will be a residue.\n\nConsistent with the Law of Large "
"Numbers,\na large T (e.g. >1000) is needed to converge to "
"this result.\n", nt, primes[np-1], B);
puts(HR);
for (j = 0; j < nt; ++j) {
/* require N > P */
do {
mpz_urandomm(N, rst, B);
} while (mpz_cmp_ui(N, primes[np-1]) <= 0);
for (i = 0; i < np; ++i) {
mpz_set_ui(P, primes[i]);
nqr[i] += (mpz_legendre(N, P) == 1);
}
printf(" %lu/%lu\r", j, nt);
fflush(stdout);
}
col0_w = snprintf(NULL, 0, "%lu", primes[np-1]);
printf("\n\n%*s%*s%s\n", col0_w, col[0], COL_SPC, " ", col[1]);
row_w = col0_w + COL_SPC + strlen(col[1]);
for (i = 0; i < row_w; ++i) {
putchar('-');
}
putchar('\n');
for (i = 0; i < np; ++i) {
printf("%*lu%*s%0.9f (%lu/%lu)\n", col0_w, primes[i], COL_SPC,
" ", (double)nqr[i]/nt, nqr[i], nt);
}
out_free:
free(nqr);
mpz_clear(N);
mpz_clear(P);
return err;
}
int test_probab_prime() {
int err = ENONE;
mpz_t N, P;
int prob_psq;
int is_psq;
ul i, fp;
mpz_init(N);
mpz_init(P);
if ((err = set_conditions())) {
goto out_clear;
}
gmp_printf("Testing %lu random values of N, %lu < N <= %Zd...\n", nt,
primes[np-1], B);
for (i = fp = 0; i < nt; ++i) {
/* require N > P */
do {
mpz_urandomm(N, rst, B);
} while (mpz_cmp_ui(N, primes[np-1]) <= 0);
if ((prob_psq = probab_perfect_square(N, &err)) < 0) {
putchar('\n');
goto out_clear;
}
is_psq = mpz_perfect_square_p(N);
fp += (prob_psq != is_psq);
printf(" %lu/%lu\r", i, nt);
fflush(stdout);
}
printf("\n\n# primes : %lu\n", np);
printf("# N's tested : %lu\n", nt);
printf("# false positives : %lu\n", fp);
printf("success rate : %0.9f\n", (1.0 - ((double)fp / nt)));
out_clear:
mpz_clear(N);
mpz_clear(P);
return err;
}
int main(int argc, char **argv) {
int err = ENONE;
int tt;
gmp_randinit_mt(rst); /* mersenne twister */
gmp_randseed_ui(rst, time(NULL));
mpz_init_set_ui(B, 0);
printf("Test Mode\n\n"
" 0. Statistics for N > P being a quadratic residue mod P\n"
" 1. Probable perfect square algorithm\n\n"
"Enter test [0|1] ");
scanf("%d", &tt);
if (tt == 0) {
err = test_residue_stats();
} else if (tt == 1) {
err = test_probab_prime();
} else {
printf("Invalid test mode (%d)\n", tt);
err = doerr(EINVAL);
}
mpz_clear(B);
return err;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.