mersenne.cal

早期freebsd实现

/*
 * Copyright (c) 1993 David I. Bell
 * Permission is granted to use, distribute, or modify this source,
 * provided that this copyright notice remains intact.
 *
 * Perform a primality test of 2^p-1, for prime p>1.
 */

define mersenne(p)
{
	local u, i, p_mask;

	/* firewall */
	if (! isint(p))
		quit "p is not an integer";

	/* two is a special case */
	if (p == 2)
		return 1;

	/* if p is not prime, then 2^p-1 is not prime */
	if (! ptest(p,10))
		return 0;

	/* calculate 2^p-1 for later mods */
	p_mask = 2^p - 1;

	/* lltest: u(i+1) = u(i)^2 - 2 mod 2^p-1 */
	u = 4;
	for (i = 2; i < p; ++i) {
		u = u^2 - 2;
		u = u&p_mask + u>>p;
		if (u > p_mask)
			u = u&p_mask + 1;
	}

	/* 2^p-1 is prime iff u(p) = 0 mod 2^p-1 */
	return (u == p_mask);
}

global lib_debug;
if (lib_debug >= 0) {
    print "mersenne(p) defined";
}