mirror of
https://github.com/rn10950/RetroZilla.git
synced 2024-09-19 19:50:18 +02:00
269 lines
6.5 KiB
C
269 lines
6.5 KiB
C
/*
|
|
* 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;
|
|
}
|