prime_sieve

[Marco Bodrato]

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
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 3 of the License, or (at your
option) any later version.

see https://www.gnu.org/licenses/.  */

#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					\
	}							\
      __mask = __mask << 1 | __mask >> (GMP_LIMB_BITS-1);	\
      __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/