primes.c

00001 /* Compute and print the prime numbers
00002    Copyright (C) 2001, 2002, 2003 Free Software Foundation, Inc.
00003    Written by Stephane Carrez (stcarrez@nerim.fr)       
00004 
00005 This file is free software; you can redistribute it and/or modify it
00006 under the terms of the GNU General Public License as published by the
00007 Free Software Foundation; either version 2, or (at your option) any
00008 later version.
00009 
00010 In addition to the permissions in the GNU General Public License, the
00011 Free Software Foundation gives you unlimited permission to link the
00012 compiled version of this file with other programs, and to distribute
00013 those programs without any restriction coming from the use of this
00014 file.  (The General Public License restrictions do apply in other
00015 respects; for example, they cover modification of the file, and
00016 distribution when not linked into another program.)
00017 
00018 This file is distributed in the hope that it will be useful, but
00019 WITHOUT ANY WARRANTY; without even the implied warranty of
00020 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00021 General Public License for more details.
00022 
00023 You should have received a copy of the GNU General Public License
00024 along with this program; see the file COPYING.  If not, write to
00025 the Free Software Foundation, 59 Temple Place - Suite 330,
00026 Boston, MA 02111-1307, USA.  */
00027 
00054 
00055 #include <stdio.h>
00056 #include <stdarg.h>
00057 #include <sys/param.h>
00058 #include <sys/sio.h>
00059 #include <string.h>
00060 #include <imath.h>
00061 
00062 /* The DATA_SIZE depends on your board.  The more memory, the largest
00063    the bitmap can be.  The last number we can record is 8 times the
00064    size of that bitmap (1 bit for each number).  */
00065 
00066 /* Remove the following test that reduces the bitmap.
00067    With a 32K memory system, the bitmap can contain up to
00068    224768 bits and the computation takes... a lot of time... */
00069 #if DATA_SIZE > 8192 + 512
00070 # undef DATA_SIZE
00071 # define DATA_SIZE 8192 + 512
00072 #endif
00073 
00074 /* Leave 512 or 128 bytes for stack.  */
00075 #if DATA_SIZE < 512
00076 # define MAX_PRIMES (DATA_SIZE - 128)
00077 #else
00078 # define MAX_PRIMES (DATA_SIZE - 512)
00079 #endif
00080 
00081 #define LAST_PRIME (MAX_PRIMES * 8L)
00082 
00083 /* Bitmap representing the prime numbers.  */
00084 unsigned char prime_list[MAX_PRIMES];
00085 
00086 /* Returns true if 'n' is a prime number recorded in the table.  */
00087 static inline int
00088 is_prime (unsigned long n)
00089 {
00090   unsigned short bit = (unsigned short) (n) & 0x07;
00091   
00092   return prime_list[n >> 3] & (1 << bit);
00093 }
00094 
00095 /* Record 'n' as a prime number in the table.  */
00096 static inline void
00097 set_prime (unsigned long n)
00098 {
00099   unsigned short bit = (unsigned short) (n) & 0x07;
00100 
00101   prime_list[n >> 3] |= (1 << bit);
00102 }
00103 
00104 /* Check whether 'n' is a prime number.
00105    Returns 1 if it's a prime number, 0 otherwise.  */
00106 static int
00107 check_for_prime (unsigned long n)
00108 {
00109   unsigned long i;
00110   unsigned char *p;
00111   unsigned long last_value;
00112   unsigned char small_n;
00113 
00114   small_n = (n & 0xffff0000) == 0;
00115   i = 0;
00116 
00117   /* We can stop when we have checked all prime numbers below sqrt(n).  */
00118   last_value = lsqrt (n);
00119 
00120   /* Scan the bitmap of prime numbers and divide 'n' by the corresponding
00121      prime to see if it's a multiple of it.  */
00122   p = prime_list;
00123   do
00124     {
00125       unsigned char val;
00126       
00127       val = *p++;
00128       if (val)
00129         {
00130           unsigned short j;
00131           unsigned long q;
00132 
00133           q = i;
00134           for (j = 1; val && j <= 0x80; j <<= 1, q++)
00135             {
00136               if (val & j)
00137                 {
00138                   val &= ~j;
00139 
00140                   /* Use 16-bit division if 'n' is small enough.  */
00141                   if (small_n)
00142                     {
00143                       unsigned short r;
00144 
00145                       /* 'n' is a multiple of prime 'q'.  */
00146                       r = (unsigned short) (n) % (unsigned short) (q);
00147                       if (r == 0)
00148                         return 0;
00149                     }
00150                   else
00151                     {
00152                       unsigned long r;
00153                       
00154                       r = n % q;
00155 
00156                       /* 'n' is a multiple of prime 'q'.  */
00157                       if (r == 0)
00158                         return 0;
00159                     }
00160                 }
00161             }
00162         }
00163       i += 8;
00164     }
00165   while (i < last_value);
00166   return 1;
00167 }
00168 
00169 #if 0
00170 /* Utility function that can be called from gdb to dump the prime
00171    numbers.  Do:  `call print_primes()' from gdb.  */
00172 static void
00173 print_primes (void)
00174 {
00175   long i;
00176   
00177   for (i = 0; i < LAST_PRIME; i++)
00178     if (is_prime (i))
00179       printf ("%ld\n", i);
00180 }
00181 #endif
00182 
00183 /* Set this variable to 1 if you want to print the prime numbers
00184    while they are found.  */
00185 int verbose = 0;
00186 
00187 int
00188 main ()
00189 {
00190   long i;
00191   short cnt = 0;
00192 
00193   serial_init ();
00194   printf ("Computing prime numbers below %ld\n",
00195           (long) LAST_PRIME);
00196   memset (prime_list, 0, sizeof (prime_list));
00197   for (i = 2; i < LAST_PRIME; i++)
00198     {
00199       if (check_for_prime (i))
00200         {
00201           set_prime (i);
00202           cnt ++;
00203           if (verbose)
00204             printf ("%ld\n", i);
00205         }
00206     }
00207   printf ("Found %ld prime numbers below %ld\n",
00208           (long) cnt, (long) LAST_PRIME);
00209 
00210   return 0;
00211 }