Comment from Eric Rains, Sep 26 2009:

      I think the attached file computes the first 2^24 terms in the EKG 
sequence A064413 (the code's about 8 years old, though, so I'm not really sure 
exactly what it's doing).

/* compile this with gcc -O3 -Wall */
#include <stdio.h>

#if 1
/* stop when a(n)>AMAX */
#define AMAX (0x4000000)
/* square root of AMAX */
#define SAMAX (0x2000)
#else
#define AMAX (0x10000)
#define SAMAX (0x100)
#endif

/* this is done as a structure to minimize the memory access overhead; */
/* one memory access should suffice for all                            */
typedef struct {
    /* has n been used? */
  unsigned int hit;
    /* next available multiple of n, if n prime */
  unsigned int gap;
    /* smallest prime factor of n */
  unsigned int small;
    /* largest factor of n not divisible by small(n) */
  unsigned int quot;} info;
  
info scratch[AMAX];

#define hit(n)   (scratch[(n)-1].hit)
#define gap(n)   (scratch[(n)-1].gap)
#define small(n) (scratch[(n)-1].small)
#define quot(n)  (scratch[(n)-1].quot)

int main(void)
{
int i,j,p;
int n=1;
int anow;
int old=0;

for (i=1;i<=AMAX;i++) {
  hit(i)=0;
  gap(i)=i;
  small(i)=i;
  quot(i)=1;}

setbuf(stdout,NULL);

/* we're sieving anyway, so why bother with a list of primes? */
for (p=2;p<=SAMAX;p++) {
  if (small(p)==p) {
    for (i=p+p,j=2;i<=AMAX;i+=p,j++) {
      if (small(i)>=p) {
        small(i)=p;
        if (small(j)==p) quot(i)=quot(j);
        else quot(i)=j;
        gap(i)=0;}}}}

anow=2;
hit(anow)=1;
/*printf("a(1)=2\n");*/
/*printf("1\n");*/

while (anow<=AMAX) {
  int t1,t2,t3,anext;
  n++;
  t1=anow;
  anext=0;
  do {
    t2=small(t1);
    t1=quot(t1);
    t3=gap(t2);
    if (t3==anow) {
      while ((t3<=AMAX)&&hit(t3)) t3+=t2;
      gap(t2)=t3;}
    if ((anext==0)||(t3<anext)) anext=t3;} while (t1!=1);
  if (anext<=AMAX) hit(anext)=1;

  old=anow;
  anow=anext;
  /* if (anext==n+1) printf("%d\n",n+1);*/
  printf("a(%d)=%d\n",n,anext);}
printf("%d %d\n",n,anow);
return(0);
}