User:Ledrug/bits: Difference between revisions
No edit summary |
No edit summary |
||
Line 75: | Line 75: | ||
#include <stdlib.h> |
#include <stdlib.h> |
||
#include <stdint.h> |
#include <stdint.h> |
||
#include <string.h> |
|||
typedef |
typedef uint32_t pint; |
||
typedef uint64_t xint; |
typedef uint64_t xint; |
||
typedef unsigned int uint; |
|||
uint8_t *pbits; |
|||
pint sieved; |
|||
#define |
#define MAX_PRIME (~(uint32_t)1) |
||
#define |
#define MAX_PRIME_SQ 65535U |
||
#define PBITS (MAX_PRIME * 8 / 30 + 1) |
|||
pint next_prime(pint); |
|||
int is_prime(pint); |
|||
void sieve(pint); |
|||
uint8_t bit_pos[30] = { |
|||
#define P_OFFSET(x) (x / UBITS) |
|||
0, 1<<0, 0, 0, 0, 0, |
|||
0, 1<<1, 0, 0, 0, 1<<2, |
|||
0, 1<<3, 0, 0, 0, 1<<4, |
|||
0, 1<<5, 0, 0, 0, 1<<6, |
|||
0, 0, 0, 0, 0, 1<<7, |
|||
}; |
|||
uint8_t rem_num[] = { 1, 7, 11, 13, 17, 19, 23, 29 }; |
|||
#define MAX_PRIME (1ULL << 32) |
|||
#define MAX_PRIME_SQ (1ULL << 16) |
|||
uint *pbits; |
|||
void |
void init_primes() |
||
{ |
{ |
||
pint tgt = 4; |
|||
pbits = malloc(PBITS); |
|||
off_t ofs; |
|||
pbits = calloc(sizeof(uint), (MAX_PRIME / UBITS)); |
|||
FILE * |
FILE *fp = fopen("primebits", "r"); |
||
if ( |
if (fp) { |
||
fread(pbits, |
fread(pbits, 1, PBITS, fp); |
||
sieved = MAX_PRIME_SQ; |
|||
fclose(pf); |
|||
fclose(fp); |
|||
return; |
return; |
||
} |
} |
||
memset(pbits, 255, PBITS); |
|||
sieved = 5; |
|||
for (p = 3; p <= 7; p += 2) { |
|||
while (sieved < MAX_PRIME_SQ) { |
|||
for (q = p * p; q < 105 * UBITS; q += p * 2) { |
|||
if (sieved > tgt) { |
|||
ofs = P_OFFSET(q); |
|||
tgt *= 2; |
|||
fprintf(stderr, "sieve %u\n", sieved); |
|||
} |
} |
||
sieve(sieved); |
|||
} |
} |
||
fp = fopen("primebits", "w"); |
|||
for (p = 105; p < MAX_PRIME / UBITS; p ++) |
|||
fwrite(pbits, 1, PBITS, fp); |
|||
fclose(fp); |
|||
} |
|||
void sieve(pint p) |
|||
for (p = 11; p <= MAX_PRIME_SQ;) { |
|||
{ |
|||
for (q = p * p; q < MAX_PRIME; q += p * 2) { |
|||
xint p1 = p * p; |
|||
ofs = P_OFFSET(q); |
|||
pint o = sieved; |
|||
pbits[ofs] |= P_BIT(q); |
|||
unsigned char bits; |
|||
} |
|||
uint32_t addr; |
|||
for (p += 2; p <= MAX_PRIME_SQ; p += 2) { |
|||
off_t shift[8]; |
|||
if (!(pbits[ P_OFFSET(p) ] & P_BIT(p))) break; |
|||
int n_shift = 0, i; |
|||
if (sieved > 5) |
|||
while (sieved < p + 30) { |
|||
addr = p1 / 30; |
|||
bits = bit_pos[p1 % 30]; |
|||
if (!n_shift || addr > shift[n_shift - 1]) |
|||
shift[n_shift++] = addr; |
|||
pbits[addr] &= ~bits; |
|||
sieved = next_prime(sieved); |
|||
p1 = p * sieved; |
|||
} |
|||
sieved = next_prime(o); |
|||
i = 0; |
|||
addr = o = shift[0] + p; |
|||
while (addr < PBITS) { |
|||
pbits[addr] = pbits[addr - p]; |
|||
addr = o + shift[i++]; |
|||
if (i == n_shift) { |
|||
i = 0; |
|||
o += p; |
|||
} |
} |
||
} |
} |
||
pf = fopen("primebits", "w"); |
|||
fwrite(pbits, sizeof(uint), MAX_PRIME / UBITS, pf); |
|||
fclose(pf); |
|||
} |
} |
||
pint next_prime(pint p) |
|||
int is_prime(xint x) |
|||
{ |
{ |
||
off_t addr; |
|||
uint8_t bits, rem; |
|||
if (x == 2) return 1; |
|||
if ((x & 1) == 0) return 0; |
|||
if ( |
if (p > 5) { |
||
addr = p / 30; |
|||
return !(pbits[ P_OFFSET(x) ] & P_BIT(x)); |
|||
bits = bit_pos[ p % 30 ] << 1; |
|||
for (rem = 0; (1 << rem) < bits; rem++); |
|||
while (pbits[addr] < bits || !bits) { |
|||
return 0; |
|||
addr++; |
|||
return 1; |
|||
bits = 1; |
|||
rem = 0; |
|||
} |
|||
if (addr >= PBITS) return 0; |
|||
while (!(pbits[addr] & bits)) { |
|||
rem++; |
|||
bits <<= 1; |
|||
} |
|||
p = addr * 30 + rem_num[rem]; |
|||
return p; |
|||
} |
|||
switch(p) { |
|||
case 2: return 3; |
|||
case 3: return 5; |
|||
case 5: return 7; |
|||
} |
|||
return 2; |
|||
} |
} |
||
int decompose(xint |
int decompose(xint n, xint *f) |
||
{ |
{ |
||
pint p = 0; |
|||
int i = 0; |
|||
while (n >= (xint)p * p) { |
|||
if (!(p = next_prime(p))) break; |
|||
while (n % p == 0) { |
|||
n /= p; |
|||
if (!is_prime(p)) continue; |
|||
while (x % p == 0) { |
|||
x /= p; |
|||
f[i++] = p; |
f[i++] = p; |
||
} |
} |
||
} |
} |
||
if ( |
if (n > 1) f[i++] = n; |
||
return i; |
return i; |
||
} |
} |
||
Line 157: | Line 210: | ||
int main() |
int main() |
||
{ |
{ |
||
int i, |
int i, len; |
||
pint p = 0; |
|||
xint fac[100], x; |
|||
xint f[100], po; |
|||
sieve(); |
|||
init_primes(); |
|||
while ((p = next_prime(p)) < 63) { |
|||
po = (1LLU << p) - 1; |
|||
printf("2^%u = %llu = ", p, po); |
|||
printf("2^%d - 1 = %llu = ", i, x); |
|||
fflush(stdout); |
fflush(stdout); |
||
n = decompose(x, fac); |
|||
len = decompose(po, f); |
|||
while (n--) printf("%llu%s", fac[n], n ? " x " : "\n"); |
|||
for (i = 0; i < len; i++) |
|||
printf("%llu%s", f[i], i == len - 1 ? "\n" : " x "); |
|||
} |
} |
||
return 0; |
return 0; |
Revision as of 00:44, 5 August 2011
<lang c>#include <stdio.h>
- include <stdlib.h>
- include <stdint.h>
typedef uint32_t pint; typedef uint64_t xint; typedef unsigned int uint;
int is_prime(xint);
inline int next_prime(pint p) { if (p == 2) return 3; for (p += 2; p > 1 && !is_prime(p); p += 2); if (p == 1) return 0; return p; }
int is_prime(xint n) {
- define NCACHE 8192
- define S (sizeof(uint) * 2)
static uint cache[NCACHE] = {0};
pint p = 2; int ofs, bit = -1;
if (n < NCACHE * S) { ofs = n / S; bit = 1 << ((n & (S - 1)) >> 1); if (cache[ofs] & bit) return 1; }
do { if (n % p == 0) return 0; if (p * p > n) break; } while ((p = next_prime(p)));
if (bit != -1) cache[ofs] |= bit; return 1; }
int decompose(xint n, pint *out) { int i = 0; pint p = 2; while (n > p * p) { while (n % p == 0) { out[i++] = p; n /= p; } if (!(p = next_prime(p))) break; } if (n > 1) out[i++] = n; return i; }
int main() { int i, j, len; xint z; pint out[100]; for (i = 2; i < 64; i = next_prime(i)) { z = (1ULL << i) - 1; printf("2^%d - 1 = %llu = ", i, z); fflush(stdout); len = decompose(z, out); for (j = 0; j < len; j++) printf("%u%s", out[j], j < len - 1 ? " x " : "\n"); }
return 0; }</lang> <lang c>#include <stdio.h>
- include <stdlib.h>
- include <stdint.h>
- include <string.h>
typedef uint32_t pint; typedef uint64_t xint; typedef unsigned int uint;
uint8_t *pbits; pint sieved;
- define MAX_PRIME (~(uint32_t)1)
- define MAX_PRIME_SQ 65535U
- define PBITS (MAX_PRIME * 8 / 30 + 1)
pint next_prime(pint); int is_prime(pint); void sieve(pint);
uint8_t bit_pos[30] = { 0, 1<<0, 0, 0, 0, 0, 0, 1<<1, 0, 0, 0, 1<<2, 0, 1<<3, 0, 0, 0, 1<<4, 0, 1<<5, 0, 0, 0, 1<<6, 0, 0, 0, 0, 0, 1<<7, };
uint8_t rem_num[] = { 1, 7, 11, 13, 17, 19, 23, 29 };
void init_primes() { pint tgt = 4; pbits = malloc(PBITS);
FILE *fp = fopen("primebits", "r"); if (fp) { fread(pbits, 1, PBITS, fp); sieved = MAX_PRIME_SQ; fclose(fp); return; }
memset(pbits, 255, PBITS); sieved = 5; while (sieved < MAX_PRIME_SQ) { if (sieved > tgt) { tgt *= 2; fprintf(stderr, "sieve %u\n", sieved); } sieve(sieved); } fp = fopen("primebits", "w"); fwrite(pbits, 1, PBITS, fp); fclose(fp); }
void sieve(pint p) { xint p1 = p * p; pint o = sieved; unsigned char bits; uint32_t addr; off_t shift[8]; int n_shift = 0, i;
if (sieved > 5) while (sieved < p + 30) { addr = p1 / 30; bits = bit_pos[p1 % 30]; if (!n_shift || addr > shift[n_shift - 1]) shift[n_shift++] = addr;
pbits[addr] &= ~bits; sieved = next_prime(sieved); p1 = p * sieved; } sieved = next_prime(o);
i = 0; addr = o = shift[0] + p; while (addr < PBITS) { pbits[addr] = pbits[addr - p]; addr = o + shift[i++]; if (i == n_shift) { i = 0; o += p; } }
}
pint next_prime(pint p) { off_t addr; uint8_t bits, rem;
if (p > 5) { addr = p / 30; bits = bit_pos[ p % 30 ] << 1; for (rem = 0; (1 << rem) < bits; rem++); while (pbits[addr] < bits || !bits) { addr++; bits = 1; rem = 0; } if (addr >= PBITS) return 0; while (!(pbits[addr] & bits)) { rem++; bits <<= 1; } p = addr * 30 + rem_num[rem]; return p; } switch(p) { case 2: return 3; case 3: return 5; case 5: return 7; } return 2; }
int decompose(xint n, xint *f) { pint p = 0; int i = 0; while (n >= (xint)p * p) { if (!(p = next_prime(p))) break; while (n % p == 0) { n /= p; f[i++] = p; } } if (n > 1) f[i++] = n; return i; }
int main() { int i, len; pint p = 0; xint f[100], po; init_primes();
while ((p = next_prime(p)) < 63) { po = (1LLU << p) - 1; printf("2^%u = %llu = ", p, po); fflush(stdout);
len = decompose(po, f); for (i = 0; i < len; i++) printf("%llu%s", f[i], i == len - 1 ? "\n" : " x "); } return 0; }</lang>