From: rogers@rosencrantz.osf.org (Andrew Rogers)
Subject: Re: Ways to annoy the entire Net
Organization: HICAM/OSD
Date: Mon, 16 Nov 1992 16:18:16 GMT

In article <1992Nov13.203344.27430@osf.org> rogers@rosencrantz.osf.org (Andrew Rogers) writes:
>I once posted the largest Mersenne prime to sci.math, which pissed a lot of
>people off (and they've discovered two or three larger ones since then)...

A couple of people have asked about this, so here's the program I used to
generate the aforementioned net.annoyer:

/* 
 * Mersenne prime program
 *
 * Mersenne primes are primes of the form (2^p)-1, where p itself is prime; as
 * of 5/92, 32 such primes are known to exist.  This program prints the value
 * of any Mersenne prime in decimal; for efficiency's sake, calculations are
 * performed using multiple-precision arithmetic in base 10000, multiplying
 * the intermediate result by 2^13 (8192) on each iteration.  (These values
 * were chosen for a 32-bit machine; see "machine dependencies" below.)
 *
 * Usage: a.out <n>     (1 <= n <= 32)
 *
 * Compile-time flags: -DWIDTH=n sets width of output line to n (default = 72);
 * -DFULL prints additional information before the result.
 *
 * Author: A.W. Rogers, based on an idea in "Computer Language" circa 1985.
 */

#include <stdio.h>

/* machine dependencies - values below assume long = 32 bits */
#define BASE    10000L  /* largest power of 10 < sqrt(ULONG_MAX) */
#define EXP     4       /* LOG10(BASE) */
#define POWER   13      /* largest power of 2 < BASE */

/* other dependencies */
#define MAX_N   (sizeof(power)/sizeof(power[0]))        /* largest index */
#define MAX_P   756839                  /* largest value of p (cf. power[]) */
#define NDIG    ((MAX_P * 302)/1000)    /* approx. # of digits in 2^MAX_P */

#ifndef WIDTH
#define WIDTH   72              /* output width (columns per line) */
#endif

/* depending on your hardware, other types may be more efficient here */
typedef long word_t;            /* values 0..BASE */
typedef long tmp_t;             /* values 0..BASE^2 */

main(argc, argv)
	int argc;
	char **argv;
{
	/* powers of two corresponding to known Mersenne primes (as of 5/92) */
	static long power[] = {
		2,      3,      5,      7,      13,     17,     19,     31,
		61,     89,     107,    127,    521,    607,    1279,   2203,
		2281,   3217,   4253,   4423,   9689,   9941,   11213,  19937,
		21701,  23209,  44497,  86243,  110503, 132049, 216091, 756839
       };

	word_t val[(NDIG / EXP) + 1], carry;
	tmp_t tmp;
	long p, i, j, nwds;
	int n, atoi();
	char fmt1[10], fmt2[15], line[WIDTH + EXP + 1], *s, *t;

	/* validate argument - must be valid Mersenne prime index */
	if (argc != 2 || (n = atoi(argv[1])) < 1 || n > MAX_N) {
		fprintf(stderr, "usage: %s 1..%d\n", argv[0], MAX_N);
		exit(1);
	}

	/* calculate (2^p)-1 using multiple-precision arithmetic */
	p = power[n-1];
	val[0] = 1L << (p % POWER);                     /* starting value */
	for (nwds = 1, i = p % POWER; i < p; i += POWER) {
		for (carry = j = 0; j < nwds; j++) {
			tmp = ((tmp_t)val[j] << POWER) + carry;
			carry = tmp / BASE;
			val[j] = tmp % BASE;
		}
		if (carry)                              /* new word? */
			val[nwds++] = carry;
	}
	val[0]--;               /* subtract 1 from LS word (can't carry) */

	/* set up output formats and print the result */
	sprintf(fmt1, "%%0%dd", EXP);                   /* e.g. "%04d" */
	sprintf(fmt2, "%%%d.%ds\n", WIDTH, WIDTH);      /* e.g. "%72.72s" */

#ifdef FULL             /* assumes WIDTH is large enough to hold this */
	sprintf(line, "M[%d] = (2 ^ %d) - 1 = %d", n, p, val[nwds-1]);
#else
	sprintf(line, "%d", val[nwds-1]);
#endif
	s = line + strlen(line);
	for (t = line + WIDTH, i = nwds - 2; i >= 0; i--)
		if ((s += sprintf(s, fmt1, val[i])) > t) {
			printf(fmt2, line);     /* print previous line */
			strcpy(line, t);        /* copy excess to next */
			s = line + strlen(line);
		}

	printf("%s\n", line);                   /* final partial line */
}
Newsgroups: alt.slack