From: HJSmith@ix.netcom.com (Harry Smith) Newsgroups: sci.math Subject: Re: nPrimes Date: 3 Feb 1995 21:14:53 GMT In <3gsuen$fvc@kitten.umdc.umu.se> Niklas Frykholm writes: > >schneken@Informatik.TU-Muenchen.DE (Thomas Schnekenburger) wrote: > >> I put the exercise to compute nPrimes(n), >> that is the number of primes smaller than n. >> For example nPrimes(10)=4, nPrimes(100)=25, nPrimes(1000)=168, ... > >> Can someone give me any hints about fast algorithms for computing nPrimes() ? > >How about using the seive of Erastotenes (or however it is spelled). >Start with an array of Boolean the size of n. Then use this algorithm > >1. Take the next "true" element, let's say its x. >2. Set each x'th element in the array to false. > >The numbers remaining when the algorithm is done are primes. >The algorithm is done when n is greater than the square root of n. >For the first 20 numbers it would look like this: (only true shown) > >Step 1: 02 03 04 05 06 07 08 09 10 11 12 13 14 15 16 17 18 19 20 >Step 2: 02 03 05 07 09 11 13 15 17 19 >Step 3: 02 03 05 07 11 13 17 > >5 is greater than the square root of 20, so we are done. > >I think the algorithm is very effective regarding speed, it uses >a lot of memory as a tradeoff. > >-- Niklas > Try my program: /* Start of file Prime8By.h ---------------------------------------------- */ /* Name = 'Prime8By - Generates Primes'; Version = 'Version 1.10, last revised: 1994-09-18 0600 hours'; Author = 'Copyright (c) 1981-1994 by author: Harry J. Smith,'; Address = '19628 Via Monte Dr., Saratoga CA 95070. All rights reserved.'; */ #ifndef Prime8By_H #define Prime8By_H /* Converted from Apple II+ Pascal to Turbo Pascal 5.0 */ /* and changed from SextPrime to Prime8By */ /* 90/07/08 by Harry Smith, Saratoga CA */ /* Convertec to C 1994-09-18 */ /* Algorithm 357 - Collected Algorithms from ACM */ /* An Efficient Prime Number Generator */ /* Driver can generates all the primes < 7,053,000,000.0 = 83987^2 - 816,169 */ /* 83987 is the 8190th prime */ /*****************************************************************************/ /* Developed in TURBO Pascal 5.5, converted to C */ #define JLim (8000) /* Size of working array */ #define MaxT (7053000000.0) /* Maximum integer tested */ typedef double PrimeType; extern PrimeType far *IP; /* Working storage for primes */ /* The following PROCEDUREs can be called externally: */ int NPrime(int M); /* M = Number of new primes desired this call, must be >= 2 on 1st call */ /* The current length of the tables in array IQ and JQ is returned */ /* This is approx. the number of primes < Sqrt(current prime) */ /* It must not reach QS = 8190 */ /* May set IP[1] < 0 to restart if IP != NULL */ #endif /* ifndef Prime8By_H */ /* End of file Prime8By.h ------------------------------------------------ */ /* Start of file Prime8By.C ---------------------------------------------- */ /* Name = 'Prime8By - Generates Primes'; Version = 'Version 1.10, last revised: 1994-09-18 0600 hours'; Author = 'Copyright (c) 1981-1994 by author: Harry J. Smith,'; Address = '19628 Via Monte Dr., Saratoga CA 95070. All rights reserved.'; */ #include #include #include #include #include "Prime8By.h" /* Converted from Apple II+ Pascal to Turbo Pascal 5.0 */ /* and changed from SextPrime to Prime8By */ /* 90/07/08 by Harry Smith, Saratoga CA */ /* Convertec to C 1994-09-18 */ /* Algorithm 357 - Collected Algorithms from ACM */ /* An Efficient Prime Number Generator */ /* Driver can generates all the primes < 7,053,000,000.0 = 83987^2 - 816,169 */ /* 83987 is the 8190th prime */ /*****************************************************************************/ /* Developed in TURBO Pascal 5.5, converted to C */ /* const int JLim = 8000; . Size of working array */ /* const double MaxT = 7053000000.0; . Maximum integer tested */ PrimeType *IP = NULL; /* Working storage for primes */ /**************************************************************/ /* implementation */ const int QS = 8190; /* Size of IQ and JQ arrays, must be large enough to hold */ /* all primes <= Sqrt(MaxT) */ PrimeType far *IQ; /* Secondary table of primes squared */ long far *JQ; /* Secondary table of primes < Sqrt(current prime) */ /* 2, 3, 5, 7, ... 83987 */ /*--------------------------------------*/ int NPrime(int M) /* M = Number of new primes desired this call, must be >= 2 on 1st call */ /* The current length of the tables in array IQ and JQ is returned */ /* This is approx. the number of primes < Sqrt(current prime) */ /* It must not reach QS = 8190 */ { /* Typed constants are Static variables */ static long IJ; /* Current length of the tables in array IQ and JQ */ static long IK; /* Max index into JQ to test NJ */ static long INC; /* = 2, 4, 2, 4, ... */ static long J; /* Array index */ static long NJ; /* Test number for secondary table of primes, MUST BE LONG */ static PrimeType NI; /* Test number - 2J for primes */ static PrimeType NJP; /* NJ in PrimeType format */ long I; /* Array index */ long JQI = 0; /* JQ[I] to store in IP[J], MUST BE LONG*/ long K; /* Array index */ if (IP == NULL) { if ((IP = (PrimeType *) farcalloc(JLim+1, sizeof( PrimeType))) == NULL) { printf("Not enough memory\n"); exit(1); } if ((IQ = (PrimeType *) farcalloc(QS+1, sizeof( PrimeType))) == NULL) { printf("Not enough memory\n"); exit(1); } if ((JQ = (long *) farcalloc(QS+1, sizeof( long))) == NULL) { printf("Not enough memory\n"); exit(1); } IP[1] = -1; } K = 0; if (IP[1] >= 0) goto L60; /* Set initial conditions */ for (I = 1; I <= QS; I++) { IQ[(int)I] = 0; JQ[(int)I] = 0; /* Not actually needed */ } IP[1] = 2; JQ[1] = 2; IK = 2; INC = 2; IQ[2] = 9; JQ[2] = 3; IQ[1] = 3; IJ = 3; IQ[3] = 25; JQ[3] = 5; NJ = 5; K = 1; /* Prepare to delete a sequence of composite numbers */ L10: J = K + 1; NI = IQ[1] - J - J; IQ[1] = NI + JLim + JLim; /* Diagnostic printf("IQ[1]=%.0f NI=%.0f\n", IQ[1] , NI); */ for (I = J; I <= JLim; I++) IP[(int)I] = 0; L20: I = IJ; if (IQ[(int)IJ] >= IQ[1]) goto L50; /* Extend the list of primes in array JQ counting so as to omit multiples of 2 and 3 */ L30: NJ += INC; INC = 6 - INC; if ((JQ[(int)IK+1] * JQ[(int)IK+1]) <= NJ) IK++; for (J = 3; J <= IK; J++) if (!(NJ % JQ[(int)J])) goto L30; IJ++; JQ[(int)IJ] = NJ; NJP = NJ; IQ[(int)IJ] = NJP * NJP; /* Diagnostic printf("IQ[%ld]=%.0f\n", IJ, IQ[(int)IJ]); */ goto L20; /* If J + J + NI is composite, enter its smallest prime factor in IP[J]. If J + J + NI is prime, then set IP[J] to 0 */ L40: IP[(int)J] = JQI; J += JQI; /* Put prime in all slots divisable by it */ if (J < JLim) goto L40; IQ[(int)I] = NI + J + J; /* Diagnostic printf("IQ[%ld]=%.0f J=%ld NI=%.0f\n", I, IQ[(int)I], J, NI); */ L50: I--; JQI = JQ[(int)I]; J = ((long)(IQ[(int)I] - NI)) / 2; if (J < JLim) goto L40; if (I != 1) goto L50; J = K; /* Pack the next M primes in IP[1], ..., IP[M] */ L60: J++; if (IP[(int)J] != 0) goto L60; if (J == JLim) goto L10; K++; IP[(int)K] = NI + J + J; if (K != M) goto L60; return (int)IJ; } /* NPrime */ /* End of file Prime8By.C ------------------------------------------------- */ /* Start of file PRTst8By.C ---------------------------------------------- */ char *Name = "PRTst8By - Displays one Prime at a time to test Prime8By Unit"; char *Version = "Version 1.10, last revised: 1994/09/21, 0600 hours"; char *Author = "Copyright (c) 1981-1994 by author: Harry J. Smith,"; char *Address = "19628 Via Monte Dr., Saratoga, CA 95070. All rights reserved."; #include #include #include #include "Prime8By.h" /* Algorithm 357 - Collected Algorithms from ACM */ /* An Efficient Prime Number Generator */ /* Developed in Turbo Pascal 5.0, converted to C */ /* PrimeType IP[]; . Working storage for primes, defined in Prime8By */ PrimeType II; /* Count of the number of primes generated */ int M; /* Number of primes each call */ int N; /* The current length of the internal tables returned */ char Ch; /* Character input by ReadKey */ const char ESC = 27; /*--------------------------------------*/ void TestIt( void) { while (kbhit()) Ch = getch(); if (Ch == ESC) { printf("\nDo you really want to terminate this program (Y/N)? "); Ch = getch(); if (Ch == 'Y' || Ch == 'y') Ch = 'Y'; else Ch = 'N'; printf("%c\n", Ch); if (Ch == 'Y') exit(0); } } /* TestIt */ /*--------------------------------------*/ void Init( void) /* Initialize program */ { printf(" [1;33;44m [2J"); /* Ansi.Sys ESC seq. for YELLOW on BLUE, clrscr */ printf("\n%s\n%s\n%s\n%s\n\n", Name, Version, Author, Address); printf("Type any key to continue... (or Esc to quit)\n"); Ch = getch(); if (Ch == ESC) exit(0); if (IP != NULL) IP[1] = -1; /* start or restart NPrimes */ M = 2500; N = NPrime(M); /* Cannot get just 1 the first time */ if (IP != NULL) IP[1] = -1; /* start or restart NPrimes */ N = NPrime(M); /* Cannot get just 1 the first time */ II = 0; } /* Init */ /*--------------------------------------*/ int main( void) { /* PRTst8By */ Init(); do { N = NPrime(M); N = NPrime(M); N = NPrime(M); II = II + 4 * M; printf("%d ", N); printf("%.0f-%.0f ", II, IP[M]); if (kbhit()) TestIt(); if (IP[M] < MaxT) N = NPrime(M); } while (IP[M] < MaxT); return 0; } /* PRTst8By */ /* End of file PRTst8By.C ------------------------------------------------ */ -- | Harry J. Smith | 19628 Via Monte Dr., Saratoga, CA 95070-4522, USA | Home Phone: 408 741-0406, Work Phone: 408 235-5088 (Voice Mail) | EMail: HJSmith@ix.netcom.com on the Internet via Netcom NetCruiser --