1/* 2 * mersenne - perform a primality test of 2^p-1, for prime p>1 3 * 4 * Copyright (C) 1999 David I. Bell and Landon Curt Noll 5 * 6 * Primary author: David I. Bell 7 * 8 * Calc is open software; you can redistribute it and/or modify it under 9 * the terms of the version 2.1 of the GNU Lesser General Public License 10 * as published by the Free Software Foundation. 11 * 12 * Calc is distributed in the hope that it will be useful, but WITHOUT 13 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 14 * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General 15 * Public License for more details. 16 * 17 * A copy of version 2.1 of the GNU Lesser General Public License is 18 * distributed with calc under the filename COPYING-LGPL. You should have 19 * received a copy with calc; if not, write to Free Software Foundation, Inc. 20 * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. 21 * 22 * Under source code control: 1991/05/22 21:56:36 23 * File existed as early as: 1991 24 * 25 * Share and enjoy! :-) http://www.isthe.com/chongo/tech/comp/calc/ 26 */ 27 28/* 29 * NOTE: See lucas.cal for a more general routine. 30 */ 31 32 33define mersenne(p) 34{ 35 local u, i, p_mask; 36 37 /* firewall */ 38 if (! isint(p)) 39 quit "p is not an integer"; 40 41 /* two is a special case */ 42 if (p == 2) 43 return 1; 44 45 /* if p is not prime, then 2^p-1 is not prime */ 46 if (! ptest(p,1)) 47 return 0; 48 49 /* lltest: u(i+1) = u(i)^2 - 2 mod 2^p-1 */ 50 u = 4; 51 for (i = 2; i < p; ++i) { 52 u = hnrmod(u^2 - 2, 1, p, -1); 53 } 54 55 /* 2^p-1 is prime iff u(p) = 0 mod 2^p-1 */ 56 return (u == 0); 57} 58