/*
 * sieve.c
 *
 * Finds prime numbers using the Sieve of Eratosthenes
 *
 * This implementation uses a bitmap to represent all odd integers in a
 * given range.  We iterate over this bitmap, crossing off the
 * multiples of each prime we find.  At the end, all the remaining set
 * bits correspond to prime integers.
 *
 * Here, we make two passes -- once we have generated a sieve-ful of
 * primes, we copy them out, reset the sieve using the highest
 * generated prime from the first pass as a base.  Then we cross out
 * all the multiples of all the primes we found the first time through,
 * and re-sieve.  In this way, we get double use of the memory we
 * allocated for the sieve the first time though.  Since we also
 * implicitly ignore multiples of 2, this amounts to 4 times the
 * values.
 *
 * This could (and probably will) be generalized to re-use the sieve a
 * few more times.
 *
 * ***** BEGIN LICENSE BLOCK *****
 * Version: MPL 1.1/GPL 2.0/LGPL 2.1
 *
 * The contents of this file are subject to the Mozilla Public License Version
 * 1.1 (the "License"); you may not use this file except in compliance with
 * the License. You may obtain a copy of the License at
 * http://www.mozilla.org/MPL/
 *
 * Software distributed under the License is distributed on an "AS IS" basis,
 * WITHOUT WARRANTY OF ANY KIND, either express or implied. See the License
 * for the specific language governing rights and limitations under the
 * License.
 *
 * The Original Code is the MPI Arbitrary Precision Integer Arithmetic library.
 *
 * The Initial Developer of the Original Code is
 * Michael J. Fromberger.
 * Portions created by the Initial Developer are Copyright (C) 1998
 * the Initial Developer. All Rights Reserved.
 *
 * Contributor(s):
 *
 * Alternatively, the contents of this file may be used under the terms of
 * either the GNU General Public License Version 2 or later (the "GPL"), or
 * the GNU Lesser General Public License Version 2.1 or later (the "LGPL"),
 * in which case the provisions of the GPL or the LGPL are applicable instead
 * of those above. If you wish to allow use of your version of this file only
 * under the terms of either the GPL or the LGPL, and not to allow others to
 * use your version of this file under the terms of the MPL, indicate your
 * decision by deleting the provisions above and replace them with the notice
 * and other provisions required by the GPL or the LGPL. If you do not delete
 * the provisions above, a recipient may use your version of this file under
 * the terms of any one of the MPL, the GPL or the LGPL.
 *
 * ***** END LICENSE BLOCK ***** */
/* $Id: sieve.c,v 1.3 2004/04/27 23:04:37 gerv%gerv.net Exp $ */

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

typedef unsigned char  byte;

typedef struct {
  int   size;
  byte *bits;
  long  base;
  int   next;
  int   nbits;
} sieve;

void sieve_init(sieve *sp, long base, int nbits);
void sieve_grow(sieve *sp, int nbits);
long sieve_next(sieve *sp);
void sieve_reset(sieve *sp, long base);
void sieve_cross(sieve *sp, long val);
void sieve_clear(sieve *sp);

#define S_ISSET(S, B)  (((S)->bits[(B)/CHAR_BIT]>>((B)%CHAR_BIT))&1)
#define S_SET(S, B)    ((S)->bits[(B)/CHAR_BIT]|=(1<<((B)%CHAR_BIT)))
#define S_CLR(S, B)    ((S)->bits[(B)/CHAR_BIT]&=~(1<<((B)%CHAR_BIT)))
#define S_VAL(S, B)    ((S)->base+(2*(B)))
#define S_BIT(S, V)    (((V)-((S)->base))/2)

int main(int argc, char *argv[])
{
  sieve   s;
  long    pr, *p;
  int     c, ix, cur = 0;

  if(argc < 2) {
    fprintf(stderr, "Usage: %s <width>\n", argv[0]);
    return 1;
  }

  c = atoi(argv[1]);
  if(c < 0) c = -c;

  fprintf(stderr, "%s: sieving to %d positions\n", argv[0], c);

  sieve_init(&s, 3, c);

  c = 0;
  while((pr = sieve_next(&s)) > 0) {
    ++c;
  }

  p = calloc(c, sizeof(long));
  if(!p) {
    fprintf(stderr, "%s: out of memory after first half\n", argv[0]);
    sieve_clear(&s);
    exit(1);
  }

  fprintf(stderr, "%s: half done ... \n", argv[0]);

  for(ix = 0; ix < s.nbits; ix++) {
    if(S_ISSET(&s, ix)) {
      p[cur] = S_VAL(&s, ix);
      printf("%ld\n", p[cur]);
      ++cur;
    }
  }

  sieve_reset(&s, p[cur - 1]);
  fprintf(stderr, "%s: crossing off %d found primes ... \n", argv[0], cur);
  for(ix = 0; ix < cur; ix++) {
    sieve_cross(&s, p[ix]);
    if(!(ix % 1000))
      fputc('.', stderr);
  }
  fputc('\n', stderr);

  free(p);

  fprintf(stderr, "%s: sieving again from %ld ... \n", argv[0], p[cur - 1]);
  c = 0;
  while((pr = sieve_next(&s)) > 0) {
    ++c;
  }
  
  fprintf(stderr, "%s: done!\n", argv[0]);
  for(ix = 0; ix < s.nbits; ix++) {
    if(S_ISSET(&s, ix)) {
      printf("%ld\n", S_VAL(&s, ix));
    }
  }

  sieve_clear(&s);

  return 0;
}

void sieve_init(sieve *sp, long base, int nbits)
{
  sp->size = (nbits / CHAR_BIT);

  if(nbits % CHAR_BIT)
    ++sp->size;

  sp->bits = calloc(sp->size, sizeof(byte));
  memset(sp->bits, UCHAR_MAX, sp->size);
  if(!(base & 1))
    ++base;
  sp->base = base;
  
  sp->next = 0;
  sp->nbits = sp->size * CHAR_BIT;
}

void sieve_grow(sieve *sp, int nbits)
{
  int  ns = (nbits / CHAR_BIT);

  if(nbits % CHAR_BIT)
    ++ns;

  if(ns > sp->size) {
    byte   *tmp;
    int     ix;

    tmp = calloc(ns, sizeof(byte));
    if(tmp == NULL) {
      fprintf(stderr, "Error: out of memory in sieve_grow\n");
      return;
    }

    memcpy(tmp, sp->bits, sp->size);
    for(ix = sp->size; ix < ns; ix++) {
      tmp[ix] = UCHAR_MAX;
    }

    free(sp->bits);
    sp->bits = tmp;
    sp->size = ns;

    sp->nbits = sp->size * CHAR_BIT;
  }
}

long sieve_next(sieve *sp)
{
  long out;
  int  ix = 0;
  long val;

  if(sp->next > sp->nbits)
    return -1;

  out = S_VAL(sp, sp->next);
#ifdef DEBUG
  fprintf(stderr, "Sieving %ld\n", out);
#endif

  /* Sieve out all multiples of the current prime */
  val = out;
  while(ix < sp->nbits) {
    val += out;
    ix = S_BIT(sp, val);
    if((val & 1) && ix < sp->nbits) { /* && S_ISSET(sp, ix)) { */
      S_CLR(sp, ix);
#ifdef DEBUG
      fprintf(stderr, "Crossing out %ld (bit %d)\n", val, ix);
#endif
    }
  }

  /* Scan ahead to the next prime */
  ++sp->next;
  while(sp->next < sp->nbits && !S_ISSET(sp, sp->next))
    ++sp->next;

  return out;
}

void sieve_cross(sieve *sp, long val)
{
  int  ix = 0;
  long cur = val;

  while(cur < sp->base)
    cur += val;

  ix = S_BIT(sp, cur);
  while(ix < sp->nbits) {
    if(cur & 1) 
      S_CLR(sp, ix);
    cur += val;
    ix = S_BIT(sp, cur);
  }
}

void sieve_reset(sieve *sp, long base)
{
  memset(sp->bits, UCHAR_MAX, sp->size);
  sp->base = base;
  sp->next = 0;
}

void sieve_clear(sieve *sp)
{
  if(sp->bits) 
    free(sp->bits);

  sp->bits = NULL;
}