2004-09-17

primes in e

Here's how to find the first 10-digit prime in e. If you get hired, you're welcome.

The following could probably be simplified to use stdio, and you could use any of a dozen different arbitrary precision mathematics libraries. I use GMP. This program prints successive 10-digit values of e until it finds a possible prime, then quits. To skip the composites, try "./eprime | tail -1". This code runs almost instantly on my iBook G4. Of course, your actual mileage may vary.

#include "buffer.h"
#include "stralloc.h"
#include "strerr.h"
#include <gmp.h>

#define FATAL "eprime: fatal: "
void die_nomem(){strerr_die2x(111,FATAL,"out of memory");}

static stralloc strnum = {0};

int main() {

  mpf_t sum,k,f,tmp;
  mpz_t n;
  long n_digit = 1000;
  const long iterations = 500;
  long l;
  int i;
  char buf[10];
  char *ch;

  if (!stralloc_ready(&strnum,n_digit)) die_nomem();

  mpf_set_default_prec(4096);
  mpf_init_set_si(sum,1);
  mpf_init(k);
  mpf_init_set_si(f,1);
  mpf_init(tmp);
  mpz_init(n);

  for(l=1;l<=iterations;++l) {
    mpf_set_si(k,l);
    mpf_mul(f,f,k);
    mpf_ui_div(tmp,1,f);
    mpf_add(sum,sum,tmp);
  }

  if (!mpf_get_str(strnum.s,(mp_exp_t *)buf,10,n_digit,sum))
    strerr_die2x(112,FATAL,"can't convert to string");
  if (!strnum.s) 
    strerr_die2x(112,FATAL,"string empty");

  mpf_clear(sum);
  mpf_clear(k);
  mpf_clear(f);
  mpf_clear(tmp);

  strnum.len = str_len(strnum.s);

  for(l=0;10+l<strnum.len;++l) {
    ch = strnum.s+l;
    for(i = 0;i<10;++i) buf[i] = ch[i];

    buffer_put(buffer_1,buf,str_len(buf));
    buffer_putsflush(buffer_1,"\n");
    mpz_set_str(n,buf,10);

    if (0<mpz_probab_prime_p(n,50)) break;
  }

  mpz_clear(n);
  return 0;
}

No comments: