freebsd-dev/usr.bin/primes/spsp.c
Colin Percival 6ae1554a5d Final step of eliminating the "games" distribution: Merge src/games
(or what's left of it, at least) into src/usr.bin.

This change will not be MFCed.

Discussed at:   EuroBSDCon 2014
Committed from: EuroBSDCon 2015
2015-10-02 10:08:11 +00:00

182 lines
4.2 KiB
C

/*-
* Copyright (c) 2014 Colin Percival
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions
* are met:
* 1. Redistributions of source code must retain the above copyright
* notice, this list of conditions and the following disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
*
* THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
* ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
* ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
* FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
* DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
* OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
* HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
* OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
* SUCH DAMAGE.
*/
#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");
#include <assert.h>
#include <stddef.h>
#include <stdint.h>
#include "primes.h"
/* Return a * b % n, where 0 <= a, b < 2^63, 0 < n < 2^63. */
static uint64_t
mulmod(uint64_t a, uint64_t b, uint64_t n)
{
uint64_t x = 0;
while (b != 0) {
if (b & 1)
x = (x + a) % n;
a = (a + a) % n;
b >>= 1;
}
return (x);
}
/* Return a^r % n, where 0 <= a < 2^63, 0 < n < 2^63. */
static uint64_t
powmod(uint64_t a, uint64_t r, uint64_t n)
{
uint64_t x = 1;
while (r != 0) {
if (r & 1)
x = mulmod(a, x, n);
a = mulmod(a, a, n);
r >>= 1;
}
return (x);
}
/* Return non-zero if n is a strong pseudoprime to base p. */
static int
spsp(uint64_t n, uint64_t p)
{
uint64_t x;
uint64_t r = n - 1;
int k = 0;
/* Compute n - 1 = 2^k * r. */
while ((r & 1) == 0) {
k++;
r >>= 1;
}
/* Compute x = p^r mod n. If x = 1, n is a p-spsp. */
x = powmod(p, r, n);
if (x == 1)
return (1);
/* Compute x^(2^i) for 0 <= i < n. If any are -1, n is a p-spsp. */
while (k > 0) {
if (x == n - 1)
return (1);
x = powmod(x, 2, n);
k--;
}
/* Not a p-spsp. */
return (0);
}
/* Test for primality using strong pseudoprime tests. */
int
isprime(ubig _n)
{
uint64_t n = _n;
/*
* Values from:
* C. Pomerance, J.L. Selfridge, and S.S. Wagstaff, Jr.,
* The pseudoprimes to 25 * 10^9, Math. Comp. 35(151):1003-1026, 1980.
*/
/* No SPSPs to base 2 less than 2047. */
if (!spsp(n, 2))
return (0);
if (n < 2047ULL)
return (1);
/* No SPSPs to bases 2,3 less than 1373653. */
if (!spsp(n, 3))
return (0);
if (n < 1373653ULL)
return (1);
/* No SPSPs to bases 2,3,5 less than 25326001. */
if (!spsp(n, 5))
return (0);
if (n < 25326001ULL)
return (1);
/* No SPSPs to bases 2,3,5,7 less than 3215031751. */
if (!spsp(n, 7))
return (0);
if (n < 3215031751ULL)
return (1);
/*
* Values from:
* G. Jaeschke, On strong pseudoprimes to several bases,
* Math. Comp. 61(204):915-926, 1993.
*/
/* No SPSPs to bases 2,3,5,7,11 less than 2152302898747. */
if (!spsp(n, 11))
return (0);
if (n < 2152302898747ULL)
return (1);
/* No SPSPs to bases 2,3,5,7,11,13 less than 3474749660383. */
if (!spsp(n, 13))
return (0);
if (n < 3474749660383ULL)
return (1);
/* No SPSPs to bases 2,3,5,7,11,13,17 less than 341550071728321. */
if (!spsp(n, 17))
return (0);
if (n < 341550071728321ULL)
return (1);
/* No SPSPs to bases 2,3,5,7,11,13,17,19 less than 341550071728321. */
if (!spsp(n, 19))
return (0);
if (n < 341550071728321ULL)
return (1);
/*
* Value from:
* Y. Jiang and Y. Deng, Strong pseudoprimes to the first eight prime
* bases, Math. Comp. 83(290):2915-2924, 2014.
*/
/* No SPSPs to bases 2..23 less than 3825123056546413051. */
if (!spsp(n, 23))
return (0);
if (n < 3825123056546413051)
return (1);
/* We can't handle values larger than this. */
assert(n <= SPSPMAX);
/* UNREACHABLE */
return (0);
}