User:Ledrug/bits

From Rosetta Code

<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>
  3. include <string.h>

typedef uint32_t pint; typedef uint64_t xint; typedef unsigned int uint;

uint8_t *pbits; pint sieved;

  1. define MAX_PRIME (~(uint32_t)1)
  2. define MAX_PRIME_SQ 65535U
  3. define PBITS (MAX_PRIME * 8 / 30 + 1)

pint next_prime(pint); int is_prime(xint); 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); }

int is_prime(xint x) { pint p; if (x > 5) { if (x < MAX_PRIME) return pbits[x/30] & bit_pos[x % 30]; for (p = 2; p; p = next_prime(p)) if (x % p == 0) return 0; return 1; } return x == 2 || x == 3 || x == 5; }

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();

for (p = 1; p < 64; p++) { po = (1LLU << p) - 1; printf("2^%u - 1 = %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>

Retrieved from "https://rosettacode.org/wiki/User:Ledrug/bits?oldid=196658"