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 }