User:Ledrug/bits

<lang c>#include <stdio.h>

1. include <stdlib.h>
2. 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) {

1. define NCACHE 8192
2. 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>

1. include <stdlib.h>
2. include <stdint.h>

typedef unsigned int uint; typedef uint64_t xint;

1. define UBYTES sizeof(uint)
2. define UBITS (UBYTES * 16)

1. define P_OFFSET(x) (x / UBITS)
2. define P_BIT(x) (1 << ((x & (UBITS - 1)) >> 1))

1. define MAX_PRIME (1ULL << 32)
2. define MAX_PRIME_SQ (1ULL << 16)

uint *pbits;

void sieve() { uint64_t p, q; off_t ofs; pbits = calloc(sizeof(uint), (MAX_PRIME / UBITS)); FILE *pf = fopen("primebits", "r"); if (pf) { fread(pbits, sizeof(uint), MAX_PRIME / UBITS, pf); fclose(pf); return; }

*pbits = 1; for (p = 3; p <= 7; p += 2) { for (q = p * p; q < 105 * UBITS; q += p * 2) { ofs = P_OFFSET(q); pbits[ofs] |= P_BIT(q); } } for (p = 105; p < MAX_PRIME / UBITS; p ++) pbits[p] = pbits[p - 105];

for (p = 11; p <= MAX_PRIME_SQ;) { for (q = p * p; q < MAX_PRIME; q += p * 2) { ofs = P_OFFSET(q); pbits[ofs] |= P_BIT(q); } for (p += 2; p <= MAX_PRIME_SQ; p += 2) { if (!(pbits[ P_OFFSET(p) ] & P_BIT(p))) break; } } pf = fopen("primebits", "w"); fwrite(pbits, sizeof(uint), MAX_PRIME / UBITS, pf); fclose(pf); }

int is_prime(xint x) { xint i; if (x == 2) return 1; if ((x & 1) == 0) return 0; if (x < MAX_PRIME) return !(pbits[ P_OFFSET(x) ] & P_BIT(x)); for (i = 3; i * i < x; i++) if (!(pbits[ P_OFFSET(x) ] & P_BIT(x)) && x % i == 0) return 0; return 1; }

int decompose(xint x, xint *f) { int i; xint p; for (i = 0; !(x & 1); f[i++] = 2, x /= 2);

for (p = 3; p * p < x; p += 2) { if (!is_prime(p)) continue; while (x % p == 0) { x /= p; f[i++] = p; } } if (x > 1) f[i++] = x; return i; }

int main() { int i, n; xint fac[100], x; sieve();

for (i = 2; i < 64; i++) { if (!is_prime(i)) continue; x = (1LLU << i) - 1; printf("2^%d - 1 = %llu = ", i, x); fflush(stdout); n = decompose(x, fac); while (n--) printf("%llu%s", fac[n], n ? " x " : "\n"); } return 0; }</lang>