/*
 * factor-sieve.c:  Illustrates the "hopping factor sieve".
 * Time-stamp: <98/06/19 14:57:38 galway>
 * 
 * Compile this with the math library:  cc -lm ...
 * Further information about this program may be found at the web page
 *   http://www.math.uiuc.edu/~galway/SieveStuff/
 * or, contact the author at galway@math.uiuc.edu.
 */

#include <math.h>
#include <stdio.h>
#include <stdlib.h>

/* See code below for explanation of B1.  */
#define B1 0.261497212847643

typedef unsigned int Uns32;

/*
 * Word is ideally the "natural" word size, holding 2^ABITS_PER_WORD bits.
 */
#define ABITS_PER_WORD 5
#define BITS_PER_WORD 32
#define WORD_AMASK    0x1F
typedef Uns32                 Word;

typedef struct DivisorEntry {
  Uns32  Divisor;
  struct DivisorEntry *Next;
} DivisorEntry;

typedef struct FactorSieveEntry {
  DivisorEntry *Divisors;
  DivisorEntry *Leftovers;
} FactorSieveEntry;

typedef struct FactorSieve {
  Uns32             CurrentNumber;
  Uns32             CurrentLoc;
  /* Size in units of entries, SieveData & DivisorTable have same size.  */
  Uns32             Size;
  FactorSieveEntry *SieveData;
  DivisorEntry     *DivisorTable;
} FactorSieve;

/* Set bit n in bitarray.  bitarray is of type *Word.  */
#define SET_BIT(bitarray,n) \
        {bitarray[n>>ABITS_PER_WORD] |= (1<<(n&WORD_AMASK));}

/* Clear bit n in bitarray.  bitarray is of type *Word.  */
#define CLEAR_BIT(bitarray,n) \
        {bitarray[n>>ABITS_PER_WORD] ^= (1<<(n&WORD_AMASK));}

/* TEST_BIT n in bitarray for zero vs non-zero.  */
#define TEST_BIT(bitarray,n) (bitarray[n>>ABITS_PER_WORD]&(1<<(n&WORD_AMASK)))

/************************************************************************/

/*
 * Build an array of bits corresponding to the numbers 0..n.  Bits set
 * to 1 correspond to composite numbers <= n.
 */
Word
*BuildBitSieve(Uns32 n)
{
  Word *result, *destPtr;
  Uns32 N;
  Uns32 m, p;

  N = (n+1+ (1<< ABITS_PER_WORD)) >> ABITS_PER_WORD;

  result = malloc(N*sizeof(Word));
  destPtr = result;
  while (N-- != 0) {
    *destPtr++ = 0;
  }

  /* Cross out multiples of of primes.  */
  for (p=2; p*p <= n; ) {
    for (m=2*p; m <= n; m += p) {
      SET_BIT(result,m);
    }

    do {
      p++;
    }
    while (TEST_BIT(result,p) != 0);
  }

  return result;
}

FactorSieve
*CreateFactorSieve(Uns32 start, Uns32 divisor_bound)
{
  Word *bit_sieve;
  FactorSieve *rslt;
  DivisorEntry *divisor_table, *next;
  FactorSieveEntry *sieve_data;
  Uns32 p, k, p_offset;
  Uns32 prime_count;

  bit_sieve = BuildBitSieve(divisor_bound);
  prime_count = 0;
  for (p = 2; p <= divisor_bound; p++) {
    if(TEST_BIT(bit_sieve, p) == 0) {
      prime_count++;
    }
  }

  divisor_table = malloc(prime_count*sizeof(DivisorEntry));
  /* This is more space than we really nead for "sieve_data"?  */
  k = 0;
  for (p = 2; p <= divisor_bound; p++) {
    if(TEST_BIT(bit_sieve, p) == 0) {
      divisor_table[k].Divisor = p;
      divisor_table[k].Next = NULL;
      k++;
    }
  }

  sieve_data = malloc(prime_count*sizeof(FactorSieveEntry));
  for (k=0; k < prime_count; k++) {
    sieve_data[k].Divisors = NULL;
    sieve_data[k].Leftovers = NULL;
  }

  k = 0;
  for (p = 2; p <= divisor_bound; p++) {
    if(TEST_BIT(bit_sieve, p) == 0) {
      /* Compute offset (distance to next multiple of p).  */
      p_offset = start%p;
      if (p_offset != 0) {
          p_offset = p - p_offset;
      }

      if (p_offset < prime_count) {
        next = sieve_data[p_offset].Divisors;
        sieve_data[p_offset].Divisors = &divisor_table[k];
        divisor_table[k].Next = next;
      }
      else {
        next = sieve_data[0].Leftovers;
        sieve_data[0].Leftovers = &divisor_table[k];
        divisor_table[k].Next = next;
      }
      k++;
    }
  }

  free(bit_sieve);

  rslt = malloc(sizeof(FactorSieve));
  rslt->CurrentNumber = start;
  rslt->CurrentLoc = 0;
  rslt->Size = prime_count;
  rslt->SieveData = sieve_data;
  rslt->DivisorTable = divisor_table;

  return rslt;
}

/* Step "factor sieve" by one.  */
void
AdvanceFactorSieve(FactorSieve *theSieve)
{
  DivisorEntry *divpntr, *next_divpntr, *next, *divisors, *leftovers;
  Uns32 current_loc, new_loc;
  Uns32 p, p_offset;
  FactorSieveEntry *sieve_data;

  current_loc = theSieve->CurrentLoc;
  sieve_data = theSieve->SieveData;

  divisors = sieve_data[current_loc].Divisors;
  sieve_data[current_loc].Divisors = NULL;

  leftovers = sieve_data[current_loc].Leftovers;
  sieve_data[current_loc].Leftovers = NULL;

  /*
   * Have to be careful looping since otherwise we might modify
   * divpntr before following divpntr->Next.
   */
  divpntr = divisors;
  while (divpntr != NULL) {
    next_divpntr = divpntr->Next;
    p_offset = divpntr->Divisor;
    if (p_offset <= theSieve->Size) {
      new_loc = current_loc + p_offset;
      if (new_loc >= theSieve->Size) {
        new_loc -= theSieve->Size;
      }
      next = sieve_data[new_loc].Divisors;
      sieve_data[new_loc].Divisors = divpntr;
      divpntr->Next = next;
    }
    else {
      /* We've "bumped into" ourselves.  */
      next = sieve_data[current_loc].Leftovers;
      sieve_data[current_loc].Leftovers = divpntr;
      divpntr->Next = next;
    }
    divpntr = next_divpntr;
  }

  /* Now do similar loop over leftovers.  */
  divpntr = leftovers;
  while (divpntr != NULL) {
    next_divpntr = divpntr->Next;
    /*
     * Calculation of new_loc is more complicated here.  See
     * CreateFactorSieve for model of what to do.
     */
    p = divpntr->Divisor;
    p_offset = (theSieve->CurrentNumber)%p;
    p_offset = p - p_offset;
    
    if (p_offset <= theSieve->Size) {
      new_loc = current_loc + p_offset;
      if (new_loc >= theSieve->Size) {
        new_loc -= theSieve->Size;
      }
      next = sieve_data[new_loc].Divisors;
      sieve_data[new_loc].Divisors = divpntr;
      divpntr->Next = next;
    }
    else {
      next = sieve_data[current_loc].Leftovers;
      sieve_data[current_loc].Leftovers = divpntr;
      divpntr->Next = next;
    }
    divpntr = next_divpntr;
  }

  current_loc++;
  if (current_loc >= theSieve->Size) {
    current_loc = 0;
  }
  theSieve->CurrentLoc = current_loc;
  theSieve->CurrentNumber += 1;

  return;
}

void
PrintUsage(int argc, char *argv[])
{
  fprintf(stderr, "Usage is %s limit\n",
          argv[0]);

  fprintf(stderr,
          "Factors numbers and collects some statistics up to limit.\n"
          );

  return;
}

int
main(int argc, char *argv[])
{
  unsigned int limit;
  Uns32 bound;
  Uns32 divcount, n, p, pcount;
  double divcount_total;
  FactorSieve *fsieve;
  DivisorEntry *divisors;

  if (argc != 2) {
    PrintUsage(argc,argv);
    exit(EXIT_FAILURE);
  }

  if (sscanf(argv[1], "%d", &limit) != 1) {
    PrintUsage(argc, argv);
    exit(EXIT_FAILURE);
  }

  bound = sqrt(limit);
  if (bound*bound < limit) {
    bound += 1;
  }

  divcount_total = 0.0;
  pcount = 0;
  for (fsieve = CreateFactorSieve(2, bound);
       fsieve->CurrentNumber <= limit;
       AdvanceFactorSieve(fsieve)) {
    divisors = fsieve->SieveData[fsieve->CurrentLoc].Divisors;
    if ((divisors == NULL)
        || (divisors->Divisor == fsieve->CurrentNumber)) {
      pcount++;
    }

    /*
     * Count the number of prime divisors of CurrentNumber (ignoring
     * multiplicity, so 8 has one prime divisor, and 12 has two).
     */
    divcount = 0;
    n = fsieve->CurrentNumber;
    for ( ; divisors != NULL; divisors = divisors->Next) {
      divcount++;
      p = divisors->Divisor;
      while (n%p == 0) {
        n /= p;
      }
    }
    if (n != 1) {
      divcount++;
    }
    divcount_total += divcount;
  }

  printf("%u primes in interval 2..%u\n", pcount, limit);
  printf(" %.8f was average number of prime divisors\n",
         divcount_total/(limit-1));
  /*
   * We predict the average value based on Theorem 430 of
   * AN INTRODUCTION TO THE THEORY OF NUMBERS, 5th Edition, by Hardy
   * and Wright.
   */
  printf(" %.8f predicted.\n", B1 + log(log((double)limit)));

  exit(EXIT_SUCCESS);
}