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