/* Compute and print the prime numbers
   Copyright (C) 2001, 2002 Free Software Foundation, Inc.
   Written by Stephane Carrez (stcarrez@nerim.fr)	

This file is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2, or (at your option) any
later version.

In addition to the permissions in the GNU General Public License, the
Free Software Foundation gives you unlimited permission to link the
compiled version of this file with other programs, and to distribute
those programs without any restriction coming from the use of this
file.  (The General Public License restrictions do apply in other
respects; for example, they cover modification of the file, and
distribution when not linked into another program.)

This file is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program; see the file COPYING.  If not, write to
the Free Software Foundation, 59 Temple Place - Suite 330,
Boston, MA 02111-1307, USA.  */

/*! @page primes Compute the prime numbers using the Eratosthenes algorithm.

   This program computes the prime numbers and print them.
   It is based on the Eratosthenes algorithm.  Bascically, for a given
   number we must check that it can't be divided by the prime numbers
   below it.  In fact we can stop when we checked up to \c sqrt(N).

   The prime numbers are stored in a bitmap because they occupy
   less room in memory.  There are 6542 primes below 65536.  To keep
   them in a list, we need 13084 bytes (if we use unsigned short).
   In our case, the bitmap contains only 8192 bytes (65536 bits).

   To compute the 6542 first primes:

   <pre>
   GCC 2.95.2      293159129 cycles (146.6 s)
   GCC 2.95.3      293100142 cycles (146.6 s) 
   GCC 3.0.4       277176107 cycles (138.6 s)
   </pre>

  @htmlonly
  Source file: <a href="primes_8c-source.html">primes.c</a>
  @endhtmlonly

*/
/*@{*/

#include <stdio.h>
#include <stdarg.h>
#include <sys/param.h>
#include <string.h>
#include <imath.h>

/* The DATA_SIZE depends on your board.  The more memory, the largest
   the bitmap can be.  The last number we can record is 8 times the
   size of that bitmap (1 bit for each number).  */

/* Remove the following test that reduces the bitmap.
   With a 32K memory system, the bitmap can contain up to
   224768 bits and the computation takes... a lot of time... */
#if DATA_SIZE > 8192 + 512
# undef DATA_SIZE
# define DATA_SIZE 8192 + 512
#endif

/* Leave 512 or 128 bytes for stack.  */
#if DATA_SIZE < 512
# define MAX_PRIMES (DATA_SIZE - 128)
#else
# define MAX_PRIMES (DATA_SIZE - 512)
#endif

#define LAST_PRIME (MAX_PRIMES * 8L)

/* Bitmap representing the prime numbers.  */
unsigned char prime_list[MAX_PRIMES];

/* Returns true if 'n' is a prime number recorded in the table.  */
static inline int
is_prime (unsigned long n)
{
  unsigned short bit = (unsigned short) (n) & 0x07;
  
  return prime_list[n >> 3] & (1 << bit);
}

/* Record 'n' as a prime number in the table.  */
static inline void
set_prime (unsigned long n)
{
  unsigned short bit = (unsigned short) (n) & 0x07;

  prime_list[n >> 3] |= (1 << bit);
}

/* Check whether 'n' is a prime number.
   Returns 1 if it's a prime number, 0 otherwise.  */
static int
check_for_prime (unsigned long n)
{
  unsigned long i;
  unsigned char *p;
  unsigned long last_value;
  unsigned char small_n;

  small_n = (n & 0xffff0000) == 0;
  i = 0;

  /* We can stop when we have checked all prime numbers below sqrt(n).  */
  last_value = lsqrt (n);

  /* Scan the bitmap of prime numbers and divide 'n' by the corresponding
     prime to see if it's a multiple of it.  */
  p = prime_list;
  do
    {
      unsigned char val;
      
      val = *p++;
      if (val)
        {
          unsigned short j;
          unsigned long q;

          q = i;
          for (j = 1; val && j <= 0x80; j <<= 1, q++)
            {
              if (val & j)
                {
                  val &= ~j;

                  /* Use 16-bit division if 'n' is small enough.  */
                  if (small_n)
                    {
                      unsigned short r;

                      /* 'n' is a multiple of prime 'q'.  */
                      r = (unsigned short) (n) % (unsigned short) (q);
                      if (r == 0)
                        return 0;
                    }
                  else
                    {
                      unsigned long r;
                      
                      r = n % q;

                      /* 'n' is a multiple of prime 'q'.  */
                      if (r == 0)
                        return 0;
                    }
                }
            }
        }
      i += 8;
    }
  while (i < last_value);
  return 1;
}

#if 0
/* Utility function that can be called from gdb to dump the prime
   numbers.  Do:  `call print_primes()' from gdb.  */
static void
print_primes (void)
{
  long i;
  
  for (i = 0; i < LAST_PRIME; i++)
    if (is_prime (i))
      printf ("%ld\n", i);
}
#endif

/* Set this variable to 1 if you want to print the prime numbers
   while they are found.  */
int verbose = 0;

int
main ()
{
  long i;
  short cnt = 0;
  
  printf ("Computing prime numbers below %ld\n",
          (long) LAST_PRIME);
  memset (prime_list, 0, sizeof (prime_list));
  for (i = 2; i < LAST_PRIME; i++)
    {
      if (check_for_prime (i))
        {
          set_prime (i);
          cnt ++;
          if (verbose)
            printf ("%ld\n", i);
        }
    }
  printf ("Found %ld prime numbers below %ld\n",
          (long) cnt, (long) LAST_PRIME);

  return 0;
}
/*@}*/
