Skip to content
New issue

Have a question about this project? # for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “#”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? # to your account

Fix primality test for very small numbers #102

Open
wants to merge 5 commits into
base: main
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from 1 commit
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions src/pprods.h
Original file line number Diff line number Diff line change
Expand Up @@ -386,4 +386,6 @@ static const uint64_t pprods[] = {
0xaf157a9999de1f, // 14887*14891*14897*14923
0xb1110edec9a5e7, // 14929*14939*14947*14951
};

static const uint32_t pprods_max_prime = 14951;
#endif // PPRODS_H
19 changes: 14 additions & 5 deletions src/primetest.h
Original file line number Diff line number Diff line change
Expand Up @@ -141,9 +141,10 @@ static int is_prime_bpsw(const mpz_t n)
int min_pprods = 5, n_pprods = sizeof(pprods) / sizeof(*pprods);
mpz_t b, d;

/* Discard even numbers. */
if (!mpz_tstbit(n, 0) && mpz_cmp_ui(n, 2))
return 0;
/* Discard even and non-positive numbers. */
if (!mpz_tstbit(n, 0) || n->_mp_size < 1) {
return mpz_cmp_ui(n, 2) ? 0 : 2;
}

/* Adjust the number of GCDs to compute with products of small primes
* if bit length of n is less than 1024. As the computational cost of
Expand All @@ -155,8 +156,16 @@ static int is_prime_bpsw(const mpz_t n)
}

for (i = 0; i < n_pprods; i++) {
if (mpn_gcd_1(n->_mp_d, n->_mp_size, pprods[i]) != 1)
return 0;
if (mpn_gcd_1(n->_mp_d, n->_mp_size, pprods[i]) != 1) {
if (n->_mp_size > 1 || n->_mp_d[0] > pprods_max_prime) {
return 0;
} else if (n->_mp_d[0] < 8) {
/* Special case for primes 3, 5, 7. */
return (1 << n->_mp_d[0]) & 0xac ? 2 : 0;
} else {
break;
}
}
}

mpz_init2(b, n->_mp_size * sizeof(n->_mp_d[0]) * 8);
Expand Down
4 changes: 3 additions & 1 deletion tools/gen_pprods.py
Original file line number Diff line number Diff line change
Expand Up @@ -73,6 +73,8 @@ def print_pprods_h(prods, factors):
static const uint64_t pprods[] = {
%s
};

static const uint32_t pprods_max_prime = %d;
#endif // PPRODS_H
"""
tmpl = tmpl.strip()
Expand All @@ -81,7 +83,7 @@ def print_pprods_h(prods, factors):
factors_str = "*".join(map(str, factors[i]))
extra_spaces = " " * (16 - (pr.bit_length() + 3) // 4)
prods_str += tab + "0x%x, %s// %s\n" % (pr, extra_spaces, factors_str)
print(tmpl % (prods_str.rstrip(),))
print(tmpl % (prods_str.rstrip(), factors[-1][-1]))


n = int(sys.argv[1])
Expand Down