home *** CD-ROM | disk | FTP | other *** search
/ NetNews Usenet Archive 1992 #27 / NN_1992_27.iso / spool / alt / slack / 3123 < prev    next >
Encoding:
Text File  |  1992-11-16  |  4.2 KB  |  114 lines
  1. Newsgroups: alt.slack
  2. Path: sparky!uunet!snorkelwacker.mit.edu!paperboy.osf.org!rosencrantz!rogers
  3. From: rogers@rosencrantz.osf.org (Andrew Rogers)
  4. Subject: Re: Ways to annoy the entire Net
  5. Message-ID: <1992Nov16.161816.2610@osf.org>
  6. Sender: news@osf.org (USENET News System)
  7. Organization: HICAM/OSD
  8. References: <17249@mindlink.bc.ca> <1992Nov8.193137.13509@netcom.com> <1992Nov13.203344.27430@osf.org>
  9. Date: Mon, 16 Nov 1992 16:18:16 GMT
  10. Lines: 102
  11.  
  12. In article <1992Nov13.203344.27430@osf.org> rogers@rosencrantz.osf.org (Andrew Rogers) writes:
  13. >I once posted the largest Mersenne prime to sci.math, which pissed a lot of
  14. >people off (and they've discovered two or three larger ones since then)...
  15.  
  16. A couple of people have asked about this, so here's the program I used to
  17. generate the aforementioned net.annoyer:
  18.  
  19. /* 
  20.  * Mersenne prime program
  21.  *
  22.  * Mersenne primes are primes of the form (2^p)-1, where p itself is prime; as
  23.  * of 5/92, 32 such primes are known to exist.  This program prints the value
  24.  * of any Mersenne prime in decimal; for efficiency's sake, calculations are
  25.  * performed using multiple-precision arithmetic in base 10000, multiplying
  26.  * the intermediate result by 2^13 (8192) on each iteration.  (These values
  27.  * were chosen for a 32-bit machine; see "machine dependencies" below.)
  28.  *
  29.  * Usage: a.out <n>     (1 <= n <= 32)
  30.  *
  31.  * Compile-time flags: -DWIDTH=n sets width of output line to n (default = 72);
  32.  * -DFULL prints additional information before the result.
  33.  *
  34.  * Author: A.W. Rogers, based on an idea in "Computer Language" circa 1985.
  35.  */
  36.  
  37. #include <stdio.h>
  38.  
  39. /* machine dependencies - values below assume long = 32 bits */
  40. #define BASE    10000L  /* largest power of 10 < sqrt(ULONG_MAX) */
  41. #define EXP     4       /* LOG10(BASE) */
  42. #define POWER   13      /* largest power of 2 < BASE */
  43.  
  44. /* other dependencies */
  45. #define MAX_N   (sizeof(power)/sizeof(power[0]))        /* largest index */
  46. #define MAX_P   756839                  /* largest value of p (cf. power[]) */
  47. #define NDIG    ((MAX_P * 302)/1000)    /* approx. # of digits in 2^MAX_P */
  48.  
  49. #ifndef WIDTH
  50. #define WIDTH   72              /* output width (columns per line) */
  51. #endif
  52.  
  53. /* depending on your hardware, other types may be more efficient here */
  54. typedef long word_t;            /* values 0..BASE */
  55. typedef long tmp_t;             /* values 0..BASE^2 */
  56.  
  57. main(argc, argv)
  58.     int argc;
  59.     char **argv;
  60. {
  61.     /* powers of two corresponding to known Mersenne primes (as of 5/92) */
  62.     static long power[] = {
  63.         2,      3,      5,      7,      13,     17,     19,     31,
  64.         61,     89,     107,    127,    521,    607,    1279,   2203,
  65.         2281,   3217,   4253,   4423,   9689,   9941,   11213,  19937,
  66.         21701,  23209,  44497,  86243,  110503, 132049, 216091, 756839
  67.        };
  68.  
  69.     word_t val[(NDIG / EXP) + 1], carry;
  70.     tmp_t tmp;
  71.     long p, i, j, nwds;
  72.     int n, atoi();
  73.     char fmt1[10], fmt2[15], line[WIDTH + EXP + 1], *s, *t;
  74.  
  75.     /* validate argument - must be valid Mersenne prime index */
  76.     if (argc != 2 || (n = atoi(argv[1])) < 1 || n > MAX_N) {
  77.         fprintf(stderr, "usage: %s 1..%d\n", argv[0], MAX_N);
  78.         exit(1);
  79.     }
  80.  
  81.     /* calculate (2^p)-1 using multiple-precision arithmetic */
  82.     p = power[n-1];
  83.     val[0] = 1L << (p % POWER);                     /* starting value */
  84.     for (nwds = 1, i = p % POWER; i < p; i += POWER) {
  85.         for (carry = j = 0; j < nwds; j++) {
  86.             tmp = ((tmp_t)val[j] << POWER) + carry;
  87.             carry = tmp / BASE;
  88.             val[j] = tmp % BASE;
  89.         }
  90.         if (carry)                              /* new word? */
  91.             val[nwds++] = carry;
  92.     }
  93.     val[0]--;               /* subtract 1 from LS word (can't carry) */
  94.  
  95.     /* set up output formats and print the result */
  96.     sprintf(fmt1, "%%0%dd", EXP);                   /* e.g. "%04d" */
  97.     sprintf(fmt2, "%%%d.%ds\n", WIDTH, WIDTH);      /* e.g. "%72.72s" */
  98.  
  99. #ifdef FULL             /* assumes WIDTH is large enough to hold this */
  100.     sprintf(line, "M[%d] = (2 ^ %d) - 1 = %d", n, p, val[nwds-1]);
  101. #else
  102.     sprintf(line, "%d", val[nwds-1]);
  103. #endif
  104.     s = line + strlen(line);
  105.     for (t = line + WIDTH, i = nwds - 2; i >= 0; i--)
  106.         if ((s += sprintf(s, fmt1, val[i])) > t) {
  107.             printf(fmt2, line);     /* print previous line */
  108.             strcpy(line, t);        /* copy excess to next */
  109.             s = line + strlen(line);
  110.         }
  111.  
  112.     printf("%s\n", line);                   /* final partial line */
  113. }
  114.