# prime_sieve

Marco Bodrato bodrato at mail.dm.unipi.it
Sun Jan 3 14:58:07 UTC 2016

```Ciao,

Il Dom, 3 Gennaio 2016 3:42 pm, Marco Bodrato ha scritto:
> I attach the code.

The list does not like my attachments...

-----8<-----------8<----------
/* gmp_numberofprimesupto_ui (N) -- Returns the number the primes in
the range [0..N].
gmp_prime_sieve (S, E, *c, F) -- Calls F(c,p) for each prime p in
the range [S..E]; when the range is over or F returns a non-zero
value, it returns the number of used primes.

THE FUNCTION IN THIS FILE USE INTERNAL FUNCTIONS OF GMP, WITH A
MUTABLE INTERFACE.  IT WAS ONLY TESTED WITH GMP-6.1; IT IS ALMOST
GUARANTEED THAT IT WILL NOT WORK WITH A FUTURE GNU MP RELEASE.

Copyright 2012, 2015, 2016 Marco Bodrato.

This file is free software; you can redistribute it and/or modify it
Free Software Foundation; either version 3 of the License, or (at your
option) any later version.

#include "gmp.h"
#include "gmp-impl.h"
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

/*********************************************************/
/* Section sieve: sieving functions and tools for primes */
/*********************************************************/

static mp_limb_t
n_to_bit (mp_limb_t n) { return ((n-5)|1)/3U; }

/* id_to_n (x) = bit_to_n (x-1) = (id*3+1)|1*/
static mp_limb_t
id_to_n  (mp_limb_t id)  { return id*3+1+(id&1); }

static mp_size_t
primesieve_size (mp_limb_t n) { return n_to_bit(n) / GMP_LIMB_BITS + 1; }

/*************************************************************/
/* Section macros: common macros, for swing/fac/bin (&sieve) */
/*************************************************************/

#define LOOP_ON_SIEVE_CONTINUE(prime,end,sieve)			\
__max_i = (end);						\
\
do {							\
++__i;							\
if (((sieve)[__index] & __mask) == 0)			\
{							\
mp_limb_t prime;					\
prime = id_to_n(__i)

#define LOOP_ON_SIEVE_BEGIN(prime,start,end,off,sieve)		\
do {								\
mp_limb_t __mask, __index, __max_i, __i;			\
\
__i = (start)-(off);					\
__index = __i / GMP_LIMB_BITS;				\
__mask = CNST_LIMB(1) << (__i % GMP_LIMB_BITS);		\
__i += (off);						\
\
LOOP_ON_SIEVE_CONTINUE(prime,end,sieve)

#define LOOP_ON_SIEVE_STOP					\
}							\
__index += __mask & 1;					\
}  while (__i <= __max_i)

#define LOOP_ON_SIEVE_END					\
LOOP_ON_SIEVE_STOP;						\
} while (0)

typedef int prime_callback_func (void *ctx, unsigned long prime);

unsigned long
gmp_prime_sieve (unsigned long start, unsigned long end,
void *ctx, prime_callback_func *f)
{
unsigned long ret;

ret = (start <= 2) & (2 <= end);
if ((ret) && (*f) (ctx, 2))
return 1;
if ((start <= 3) & (3 <= end))
{
++ret;
if ((*f) (ctx, 3))
return ret;
}
if (end > 4)
{
mp_limb_t *sieve;
mp_size_t size;

size = primesieve_size (end);

sieve = __GMP_ALLOCATE_FUNC_LIMBS (size);
gmp_primesieve (sieve, end);
start = MAX (start, 5) | 1;
LOOP_ON_SIEVE_BEGIN (prime, n_to_bit(start) + (start % 3 == 0),
n_to_bit (end), 0, sieve);
++ret;
if ((*f) (ctx, prime)) break;
LOOP_ON_SIEVE_END;

__GMP_FREE_FUNC_LIMBS (sieve, size);
}
return ret;
}

unsigned long
gmp_numberofprimesupto_ui (unsigned long n)
{
static const unsigned table[] = { 0, 0, 1, 2, 2 };

if (n < numberof (table))
return table[n];
else
{
mp_limb_t *sieve;
mp_size_t size;
unsigned int ret;

size = primesieve_size (n);

sieve = __GMP_ALLOCATE_FUNC_LIMBS (size);
ret = gmp_primesieve (sieve, n);
__GMP_FREE_FUNC_LIMBS (sieve, size);

return table[numberof (table)-1] + ret;
}
}

#ifdef MAIN
int
count3mod4 (void *ctx, unsigned long prime)
{
*(unsigned long*) ctx += (prime > 2) & (prime >> 1);
return 0;
}

int
main (int argc, char *argv[])
{
unsigned long s, N, c;

s = 0;
N = strtoul (argv[1], NULL, 0);
if (argc > 2)
{
s = N;
N = strtoul (argv[2], NULL, 0);
}

c = 0;
printf ("pi(%lu .. %lu)=%lu\n", s, N, gmp_prime_sieve (s, N, &c,
&count3mod4));
printf ("pi(%lu .. %lu,=3 mod 4)=%lu\n", s, N, c);

return 0;
}
#endif

#ifdef MAINPI
int
main (int argc, char *argv[])
{
unsigned long N;

N = strtoul (argv[1], NULL, 0);

printf ("pi(0 .. %lu)=%lu\n", N, gmp_numberofprimesupto_ui (N));

return 0;
}
#endif

-----8<-----------8<----------

--
http://bodrato.it/

```