print-perfect-numbers.c

/** Print the _known_ (2013) perfect numbers less than *limit* given on stdin to stdout.

    It uses gmplib (-lgmp) to support arbitrary precision integer arithmetic.
    To compile and print all perfect numbers that are less than 100000000:
      $ gcc print-perfect-numbers.c -lgmp && echo 100000000 | ./a.out
*/
#include <stdio.h>
#include <stdlib.h>
#include <gmp.h>

/** Set *rop* to the perfect number that corresponds to *p* (Mersenne exponent).

    Even perfect = (2**p - 1) * (2**(p - 1))
    iff (2**p - 1) is Mersenne prime -- Euclid–Euler theorem [1].

    Each Mersenne prime generates one even perfect number, and vice versa.

    Odd perfect numbers are larger than 10**1500 if they exist [2].

    [1] http://en.wikipedia.org/wiki/Euclid%E2%80%93Euler_theorem
    [2] http://www.lirmm.fr/~ochem/opn/opn.pdf
*/
void perfect(unsigned long p, mpz_t rop) {
  mpz_set_ui(rop, 1);            /* rop   = 1     */
  mpz_mul_2exp(rop, rop, p);     /* rop <<= p     */
  mpz_sub_ui(rop, rop, 1);       /* rop  -= 1     */
  mpz_mul_2exp(rop, rop, p - 1); /* rop <<= p - 1 */
}

/** Print the _known_ perfect numbers less than *limit* to *fp* file. */
int print_perfect_numbers(FILE* fp, mpz_t limit) {
  static const unsigned long mersenne_exponents[] = {
    /** All known (2014) Mersenne exponents: primes p such that
        (2**p - 1) is prime. Then (2**p - 1) is called a Mersenne prime.

        https://oeis.org/A000043
     */
    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, 859433, 1257787, 1398269, 2976221, 3021377,
    6972593, 13466917, 20996011, 24036583, 25964951, 30402457,
    /* here, other numbers are possible in between

       http://en.wikipedia.org/wiki/Mersenne_prime#List_of_known_Mersenne_primes
     */
    32582657, 37156667, 42643801, 43112609, 57885161,
  };
  size_t n = sizeof(mersenne_exponents) / sizeof(*mersenne_exponents);
  const unsigned long *p = mersenne_exponents, *endptr = &mersenne_exponents[n];
  mpz_t result;
  mpz_init(result);
  for ( ; p != endptr; ++p) {
    perfect(*p, result); /* find the perfect number corresponding to p */
    if (mpz_cmp(result, limit) < 0) { /* result < limit */
      /* print result as a decimal (human readable) */
      if (mpz_out_str(fp, 10, result) == 0) { /* error */
	mpz_clear(result);
	return -1;
      }
      fputs("\n", fp);
    }
    else /* result >= limit */
      break;
  }
  mpz_clear(result);
  return (p - mersenne_exponents);
}

int main(void) {
  /* read the limit from stdin */
  mpz_t limit;
  mpz_init(limit);
  /* accept 0x (hex), 0b (binary), decimal integers */
  fprintf(stderr, "Input limit: ");
  if (mpz_inp_str(limit, stdin, 0) == 0) { /* error */
    mpz_clear(limit);
    fprintf(stderr, "Can't read input integer\n");
    exit(EXIT_FAILURE);
  }
  fprintf(stderr, "\n");

  if (print_perfect_numbers(stdout, limit) < 0) {
    mpz_clear(limit);
    fprintf(stderr, "Can't print result\n");
    exit(EXIT_FAILURE);
  }
  mpz_clear(limit);
  return 0;
}