Home | History | Annotate | Line # | Download | only in dist 
      1 /* gmp_nextprime -- generate small primes reasonably efficiently for internal
      2    GMP needs.
      3 
      4    Contributed to the GNU project by Torbjorn Granlund.  Miscellaneous
      5    improvements by Martin Boij.
      6 
      7    THE FUNCTIONS IN THIS FILE ARE INTERNAL WITH MUTABLE INTERFACES.  IT IS ONLY
      8    SAFE TO REACH THEM THROUGH DOCUMENTED INTERFACES.  IN FACT, IT IS ALMOST
      9    GUARANTEED THAT THEY WILL CHANGE OR DISAPPEAR IN A FUTURE GNU MP RELEASE.
     10 
     11 Copyright 2009 Free Software Foundation, Inc.
     12 
     13 This file is part of the GNU MP Library.
     14 
     15 The GNU MP Library is free software; you can redistribute it and/or modify
     16 it under the terms of either:
     17 
     18   * the GNU Lesser General Public License as published by the Free
     19     Software Foundation; either version 3 of the License, or (at your
     20     option) any later version.
     21 
     22 or
     23 
     24   * the GNU General Public License as published by the Free Software
     25     Foundation; either version 2 of the License, or (at your option) any
     26     later version.
     27 
     28 or both in parallel, as here.
     29 
     30 The GNU MP Library is distributed in the hope that it will be useful, but
     31 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
     32 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
     33 for more details.
     34 
     35 You should have received copies of the GNU General Public License and the
     36 GNU Lesser General Public License along with the GNU MP Library.  If not,
     37 see https://www.gnu.org/licenses/.  */
     38 
     39 /*
     40   Optimisation ideas:
     41 
     42   1. Unroll the sieving loops.  Should reach 1 write/cycle.  That would be a 2x
     43      improvement.
     44 
     45   2. Separate sieving with primes p < SIEVESIZE and p >= SIEVESIZE.  The latter
     46      will need at most one write, and thus not need any inner loop.
     47 
     48   3. For primes p >= SIEVESIZE, i.e., typically the majority of primes, we
     49      perform more than one division per sieving write.  That might dominate the
     50      entire run time for the nextprime function.  A incrementally initialised
     51      remainder table of Pi(65536) = 6542 16-bit entries could replace that
     52      division.
     53 */
     54 
     55 #include "gmp-impl.h"
     56 #include <string.h>		/* for memset */
     57 
     58 
     59 unsigned long int
     60 gmp_nextprime (gmp_primesieve_t *ps)
     61 {
     62   unsigned long p, d, pi;
     63   unsigned char *sp;
     64   static unsigned char addtab[] =
     65     { 2,4,2,4,6,2,6,4,2,4,6,6,2,6,4,2,6,4,6,8,4,2,4,2,4,8,6,4,6,2,4,6,2,6,6,4,
     66       2,4,6,2,6,4,2,4,2,10,2,10 };
     67   unsigned char *addp = addtab;
     68   unsigned long ai;
     69 
     70   /* Look for already sieved primes.  A sentinel at the end of the sieving
     71      area allows us to use a very simple loop here.  */
     72   d = ps->d;
     73   sp = ps->s + d;
     74   while (*sp != 0)
     75     sp++;
     76   if (sp != ps->s + SIEVESIZE)
     77     {
     78       d = sp - ps->s;
     79       ps->d = d + 1;
     80       return ps->s0 + 2 * d;
     81     }
     82 
     83   /* Handle the number 2 separately.  */
     84   if (ps->s0 < 3)
     85     {
     86       ps->s0 = 3 - 2 * SIEVESIZE; /* Tricky */
     87       return 2;
     88     }
     89 
     90   /* Exhausted computed primes.  Resieve, then call ourselves recursively.  */
     91 
     92 #if 0
     93   for (sp = ps->s; sp < ps->s + SIEVESIZE; sp++)
     94     *sp = 0;
     95 #else
     96   memset (ps->s, 0, SIEVESIZE);
     97 #endif
     98 
     99   ps->s0 += 2 * SIEVESIZE;
    100 
    101   /* Update sqrt_s0 as needed.  */
    102   while ((ps->sqrt_s0 + 1) * (ps->sqrt_s0 + 1) <= ps->s0 + 2 * SIEVESIZE - 1)
    103     ps->sqrt_s0++;
    104 
    105   pi = ((ps->s0 + 3) / 2) % 3;
    106   if (pi > 0)
    107     pi = 3 - pi;
    108   if (ps->s0 + 2 * pi <= 3)
    109     pi += 3;
    110   sp = ps->s + pi;
    111   while (sp < ps->s + SIEVESIZE)
    112     {
    113       *sp = 1, sp += 3;
    114     }
    115 
    116   pi = ((ps->s0 + 5) / 2) % 5;
    117   if (pi > 0)
    118     pi = 5 - pi;
    119   if (ps->s0 + 2 * pi <= 5)
    120     pi += 5;
    121   sp = ps->s + pi;
    122   while (sp < ps->s + SIEVESIZE)
    123     {
    124       *sp = 1, sp += 5;
    125     }
    126 
    127   pi = ((ps->s0 + 7) / 2) % 7;
    128   if (pi > 0)
    129     pi = 7 - pi;
    130   if (ps->s0 + 2 * pi <= 7)
    131     pi += 7;
    132   sp = ps->s + pi;
    133   while (sp < ps->s + SIEVESIZE)
    134     {
    135       *sp = 1, sp += 7;
    136     }
    137 
    138   p = 11;
    139   ai = 0;
    140   while (p <= ps->sqrt_s0)
    141     {
    142       pi = ((ps->s0 + p) / 2) % p;
    143       if (pi > 0)
    144 	pi = p - pi;
    145       if (ps->s0 + 2 * pi <= p)
    146 	  pi += p;
    147       sp = ps->s + pi;
    148       while (sp < ps->s + SIEVESIZE)
    149 	{
    150 	  *sp = 1, sp += p;
    151 	}
    152       p += addp[ai];
    153       ai = (ai + 1) % 48;
    154     }
    155   ps->d = 0;
    156   return gmp_nextprime (ps);
    157 }
    158 
    159 void
    160 gmp_init_primesieve (gmp_primesieve_t *ps)
    161 {
    162   ps->s0 = 0;
    163   ps->sqrt_s0 = 0;
    164   ps->d = SIEVESIZE;
    165   ps->s[SIEVESIZE] = 0;		/* sentinel */
    166 }
    167