Line data    Source code

       1             : /* Dynamic Eratosthenes sieve.
       2             :  
       3             :   Copyright 2001-2016 Paul Zimmermann and Alexander Kruppa.
       4             :   Imported from CADO-NFS, which imported it from GMP-ECM.
       5             :   (use 'unsigned long' instead of 'double'; thread-safe)
       6             : 
       7             :   This file is part of the ECM Library.
       8             : 
       9             :   The ECM Library is free software; you can redistribute it and/or modify
      10             :   it under the terms of the GNU Lesser General Public License as published by
      11             :   the Free Software Foundation; either version 2.1 of the License, or (at your
      12             :   option) any later version.
      13             : 
      14             :   The ECM Library is distributed in the hope that it will be useful, but
      15             :   WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
      16             :   or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
      17             :   License for more details.
      18             : 
      19             :   You should have received a copy of the GNU Lesser General Public License
      20             :   along with the ECM Library; see the file COPYING.LIB.  If not, write to
      21             :   the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston,
      22             :   MA 02110-1301, USA.
      23             : */
      24             : 
      25             : /* compile with -DMAIN to use as a standalone program */
      26             : 
      27             : #include <stdio.h>
      28             : #include <stdlib.h>
      29             : #include <string.h>
      30             : #include "getprime_r.h"
      31             : 
      32             : /* provided for in cado.h, but we want getprime.c to be standalone */
      33             : #ifndef ASSERT
      34             : #define ASSERT(x)
      35             : #endif
      36             : 
      37             : /* This function returns successive odd primes, starting with 3.
      38             :    To perform a loop over all primes <= B1, do the following
      39             :    (compile this file with -DMAIN to count primes):
      40             : 
      41             :       prime_info_t pi;
      42             :       prime_info_init (pi);
      43             :       for (p = 2; p <= B1; p = getprime_mt (pi))
      44             :          {
      45             :             ...
      46             :          }
      47             : 
      48             :       prime_info_clear (pi);
      49             : */
      50             : 
      51             : void
      52        2033 : prime_info_init (prime_info_t i)
      53             : {
      54        2033 :   i->offset = 0;
      55        2033 :   i->current = -1;
      56        2033 :   i->primes = NULL;
      57        2033 :   i->nprimes = 0;
      58        2033 :   i->sieve = NULL;
      59        2033 :   i->len = 0;
      60        2033 :   i->moduli = NULL;
      61        2033 : }
      62             : 
      63             : void
      64        2033 : prime_info_clear (prime_info_t i)
      65             : {
      66        2033 :   free (i->primes);
      67        2033 :   free (i->sieve);
      68        2033 :   free (i->moduli);
      69        2033 : }
      70             : 
      71             : /* this function is thread-safe */
      72             : ecm_uint
      73  2199833337 : getprime_mt (prime_info_t i)
      74             : {
      75  2199833337 :   if (i->len)
      76             :     {
      77  2199831310 :       unsigned char *ptr = i->sieve + i->current;
      78 22934861484 :       while (!*++ptr);
      79  2199831310 :       i->current = ptr - i->sieve;
      80             :     }
      81             :   else
      82        2027 :     i->current = 0;
      83             : 
      84  2199833337 :   if (i->current < i->len) /* most calls will end here */
      85  2198923583 :     return i->offset + 2 * i->current;
      86             : 
      87             :   /* otherwise we have to sieve */
      88      909754 :   i->offset += 2 * i->len;
      89             : 
      90             :   /* first enlarge sieving table if too small */
      91      909754 :   if ((ecm_uint) i->len * i->len < i->offset)
      92             :     {
      93       12391 :       free (i->sieve);
      94       12391 :       i->len *= 2;
      95       12391 :       i->sieve = (unsigned char *) malloc ((i->len + 1 )
      96             :                                            * sizeof (unsigned char));
      97             :       /* assume this "small" malloc will not fail in normal usage */
      98             :       ASSERT(i->sieve != NULL);
      99       12391 :       i->sieve[i->len] = 1; /* End mark */
     100             :     }
     101             : 
     102             :   /* now enlarge small prime table if too small */
     103      909754 :   if ((i->nprimes == 0) ||
     104      907727 :       ((ecm_uint) i->primes[i->nprimes - 1] * (ecm_uint)
     105      907727 :        i->primes[i->nprimes - 1] < i->offset + i->len))
     106             :       {
     107       11003 :         if (i->nprimes == 0) /* initialization */
     108             :           {
     109        2027 :             i->nprimes = 1;
     110        2027 :             i->primes = (ecm_uint*) malloc (i->nprimes
     111        2027 :                                                 * sizeof(ecm_uint));
     112             :             /* assume this "small" malloc will not fail in normal usage */
     113             :             ASSERT(i->primes != NULL);
     114        2027 :             i->moduli = (ecm_uint*) malloc (i->nprimes
     115        2027 :                                                 * sizeof(ecm_uint));
     116             :             /* assume this "small" malloc will not fail in normal usage */
     117             :             ASSERT(i->moduli != NULL);
     118        2027 :             i->len = 1;
     119        2027 :             i->sieve = (unsigned char *) malloc((i->len + 1) *
     120             :                                        sizeof(unsigned char)); /* len=1 here */
     121             :             /* assume this "small" malloc will not fail in normal usage */
     122             :             ASSERT(i->sieve != NULL);
     123        2027 :             i->sieve[i->len] = 1; /* End mark */
     124        2027 :             i->offset = 5;
     125        2027 :             i->sieve[0] = 1; /* corresponding to 5 */
     126        2027 :             i->primes[0] = 3;
     127        2027 :             i->moduli[0] = 1; /* next odd multiple of 3 is 7, i.e. next to 5 */
     128        2027 :             i->current = -1;
     129        2027 :             return 3;
     130             :           }
     131             :         else
     132             :           {
     133             :             ecm_uint k, p, j, ok;
     134             : 
     135        8976 :             k = i->nprimes;
     136        8976 :             i->nprimes *= 2;
     137        8976 :             i->primes = (ecm_uint*) realloc (i->primes, i->nprimes *
     138             :                                                  sizeof(ecm_uint));
     139        8976 :             i->moduli = (ecm_uint*) realloc (i->moduli, i->nprimes *
     140             :                                                  sizeof(ecm_uint));
     141             :             /* assume those "small" realloc's will not fail in normal usage */
     142             :             ASSERT(i->primes != NULL && i->moduli != NULL);
     143      236977 :             for (p = i->primes[k-1]; k < i->nprimes; k++)
     144             :               {
     145             :                 /* find next (odd) prime > p */
     146             :                 do
     147             :                   {
     148   389661164 :                     for (p += 2, ok = 1, j = 0; (ok != 0) && (j < k); j++)
     149   388772020 :                       ok = p % i->primes[j];
     150             :                   }
     151      889144 :                 while (ok == 0);
     152      228001 :                 i->primes[k] = p;
     153             :                 /* moduli[k] is the smallest m such that offset + 2*m = k*p */
     154      228001 :                 j = i->offset % p;
     155      228001 :                 j = (j == 0) ? j : p - j; /* -offset mod p */
     156      228001 :                 if ((j % 2) != 0)
     157      116146 :                   j += p; /* ensure j is even */
     158      228001 :                 i->moduli[k] = j / 2;
     159             :               }
     160             :           }
     161             :       }
     162             : 
     163             :   /* now sieve for new primes */
     164             :   {
     165             :     ecm_int k;
     166             :     ecm_uint j, p;
     167             :     
     168      907727 :     memset (i->sieve, 1, sizeof(unsigned char) * (i->len + 1));
     169  2732708148 :     for (j = 0; j < i->nprimes; j++)
     170             :       {
     171  2731800421 :         p = i->primes[j];
     172 51996210013 :         for (k = i->moduli[j]; k < i->len; k += p)
     173 49264409592 :           i->sieve[k] = 0;
     174  2731800421 :         i->moduli[j] = k - i->len; /* for next sieving array */
     175             :       }
     176             :   }
     177             : 
     178             :   {
     179      907727 :     unsigned char *ptr = i->sieve - 1;
     180     6968258 :     while (!*++ptr);
     181      907727 :     i->current = ptr - i->sieve;
     182             :   }
     183             : 
     184             :   ASSERT(i->current < i->len); /* otherwise we found a prime gap >= sqrt(x)
     185             :                                   around x */
     186      907727 :   return i->offset + 2 * i->current;
     187             : }
     188             : 
     189             : #ifdef MAIN
     190             : int
     191             : main (int argc, char *argv[])
     192             : {
     193             :   unsigned long p, B;
     194             :   unsigned long pi = 0;
     195             :   prime_info_t i;
     196             : 
     197             :   if (argc != 2)
     198             :     {
     199             :       fprintf (stderr, "Usage: getprime <bound>\n");
     200             :       exit (EXIT_FAILURE);
     201             :     }
     202             : 
     203             :   B = strtoul (argv[1], NULL, 0);
     204             :   
     205             :   prime_info_init (i);
     206             : 
     207             :   for (pi = 0, p = 2; p <= B; p = getprime_mt (i), pi++);
     208             :   printf ("pi(%lu)=%lu\n", B, pi);
     209             : 
     210             :   prime_info_clear (i); /* free the tables */
     211             : 
     212             :   return 0;
     213             : }
     214             : #endif