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