4.4BSD/usr/src/contrib/calc-1.26.4/lib/mersenne.cal

/*
 * 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 (!isnum(lib_debug) || lib_debug>0) print "mersenne(p) defined";