/* * 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";