Changeset b15787 in git
 Timestamp:
 Jul 1, 2010, 12:07:25 PM (13 years ago)
 Branches:
 (u'jengelhdatetime', 'ceac47cbc86fe4a15902392bdbb9bd2ae0ea02c6')(u'spielwiese', 'f875bbaccd0831e36aaed09ff6adeb3eb45aeb94')
 Children:
 a15000c3a4796528e68ed866058804cdf1e3aa8d
 Parents:
 ea947e637d9452966fb62c164fc2c965f8954886
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

Singular/misc_ip.cc
rea947e rb15787 106 106 i is assumed to be small enough to address a valid index 107 107 in the given array */ 108 bool getValue(const int i, const unsigned int* ii) 109 { 110 int index = i / 32; 111 int offset = i % 32; 108 bool getValue(const unsigned i, const unsigned int* ii) 109 { 110 if (i==2) return true; 111 if ((i & 1)==0) return false; 112 unsigned I= i/2; 113 unsigned index = I / 32; 114 unsigned offset = I % 32; 112 115 unsigned int v = 1 << offset; 113 116 return ((ii[index] & v) != 0); … … 119 122 i is assumed to be small enough to address a valid index 120 123 in the given array */ 121 void setValue(const int i, bool value, unsigned int* ii) 122 { 123 int index = i / 32; 124 int offset = i % 32; 124 void setValue(const unsigned i, bool value, unsigned int* ii) 125 { 126 if ((i&1)==0) return; // ignore odd numbers 127 unsigned I=i/2; 128 unsigned index = I / 32; 129 unsigned offset = I % 32; 125 130 unsigned int v = 1 << offset; 126 if (value && ((ii[index] & v) != 0)) return; 127 if ((!value) && ((ii[index] & v) == 0)) return; 128 if (value && ((ii[index] & v) == 0)) { ii[index] += v; return; } 129 if ((!value) && ((ii[index] & v) != 0)) { ii[index] = v; return; } 131 if (value) { ii[index] = v; } 132 else { ii[index] &= ~v; } 130 133 } 131 134 … … 161 164 (This may result in an earlier termination.) */ 162 165 163 int s = 67108864; /* = 2^26*/166 int s = 1<<25; /* = 2^25 */ 164 167 int maxP = 2147483647; /* = 2^31  1, by the way a Mersenne prime */ 165 if ((pBound != 0) && (pBound < maxP)) maxP = pBound; 168 if ((pBound != 0) && (pBound < maxP)) 169 { 170 maxP = pBound; 171 } 172 if (maxP< (214748364763) ) s=(maxP+63)/64; 173 else s=2147483647/64+1; 166 174 unsigned int* isPrime = new unsigned int[s]; 167 /* the lowest bit of isPrime[0] stores whether 0 is a prime, 168 next bit is for 1, next for 2, etc. i.e. 169 intended usage is: isPrime[0] = 10100000100010100010100010101100. 170 This goes on up to isPrime[67108863]; the highest bit of this 171 unsigned int corresponds to 67108863*32+31 = 2^311. 175 /* the lowest bit of isPrime[0] stores whether 1 is a prime, 176 next bit is for 3, next for 5, etc. i.e. 177 intended usage is: isPrime[0] = ... 172 178 We shall make use only of bits which correspond to numbers = 173 179 1 or +1 mod 6. */ 174 for (i = 0; i < s; i++) isPrime[i] = 4294967295; /* all 32 bits set */ 175 int p = 5; bool add2 = true; 180 for (i = 0; i < s; i++) isPrime[i] = ~0;/*4294967295*/; /* all 32 bits set */ 181 int p=9; while((p<maxP) && (p>0)) { setValue(p,false,isPrime); p+=6; } 182 p = 5; bool add2 = true; 176 183 /* due to possible overflows, we need to check whether p > 0, and 177 184 likewise i > 0 below */ … … 183 190 temp = divTimes(nn, p, &e); 184 191 nlDelete(&nn, NULL); nn = temp; 185 if (e > 0) { primes[pCounter] = p; multiplicities[pCounter++] = e; } 192 if (e > 0) 193 { primes[pCounter] = p; multiplicities[pCounter++] = e; } 186 194 /* invalidate all multiples of p, starting with 2*p */ 187 195 i = 2 * p; 188 while ((0 < i) && (i < = s)) { setValue(i, false, isPrime); i += p; }196 while ((0 < i) && (i < maxP)) { setValue(i, false, isPrime); i += p; } 189 197 /* move on to the next prime in the sieve; we either add 2 or 4 190 198 in order to visit just the numbers equal to 1/+1 mod 6 */
Note: See TracChangeset
for help on using the changeset viewer.