User:Ledrug/bits: Difference between revisions

From Rosetta Code
Content added Content deleted
No edit summary
No edit summary
Line 75:
#include <stdlib.h>
#include <stdint.h>
#include <string.h>
 
typedef unsigneduint32_t int uintpint;
typedef uint64_t xint;
typedef unsigned int uint;
uint8_t *pbits;
pint sieved;
 
#define UBYTESMAX_PRIME sizeof(uint~(uint32_t)1)
#define UBITSMAX_PRIME_SQ (UBYTES * 16)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)
#define 0, P_BIT(x) (1 <<0, ((x0, &0, (UBITS0, - 1)) >> 1))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 sieveinit_primes()
{
uint64_tpint p,tgt q= 4;
pbits = malloc(PBITS);
off_t ofs;
 
pbits = calloc(sizeof(uint), (MAX_PRIME / UBITS));
FILE *pffp = fopen("primebits", "r");
if (pffp) {
fread(pbits, sizeof(uint)1, MAX_PRIME / UBITSPBITS, pffp);
sieved = MAX_PRIME_SQ;
fclose(pf);
fclose(fp);
return;
}
 
*memset(pbits, =255, 1PBITS);
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);
pbits[ofs]tgt |*= P_BIT(q)2;
fprintf(stderr, "sieve %u\n", sieved);
}
sieve(sieved);
}
fp = fopen("primebits", "w");
for (p = 105; p < MAX_PRIME / UBITS; p ++)
fwrite(pbits[p], =1, pbits[p -PBITS, 105]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)
{
xintoff_t iaddr;
uint8_t bits, rem;
if (x == 2) return 1;
 
if ((x & 1) == 0) return 0;
if (xp <> MAX_PRIME5) {
addr = p / 30;
return !(pbits[ P_OFFSET(x) ] & P_BIT(x));
for (i bits = 3;bit_pos[ ip *% i30 ] << x1; i++)
iffor (!(pbits[rem P_OFFSET(x)= ]0; & P_BIT(x))1 &&<< xrem) %< ibits; == 0rem++);
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 xn, xint *f)
{
intpint ip = 0;
xintint pi = 0;
forwhile (in >= 0; !(xxint)p &* 1p); f[i++] = 2, x /= 2);{
if (!(p = next_prime(p))) break;
 
for while (pn = 3;% p * p < x; p +== 20) {
n /= p;
if (!is_prime(p)) continue;
while (x % p == 0) {
x /= p;
f[i++] = p;
}
}
if (xn > 1) f[i++] = xn;
return i;
}
Line 157 ⟶ 210:
int main()
{
int i, nlen;
pint p = 0;
xint fac[100], x;
xint f[100], po;
sieve();
init_primes();
 
forwhile (i(p = 2; inext_prime(p)) < 64; i++63) {
ifpo = (!is_prime(i)1LLU << p) continue- 1;
xprintf("2^%u = (1LLU%llu <<= i)", -p, 1po);
printf("2^%d - 1 = %llu = ", i, x);
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;

Revision as of 00:44, 5 August 2011

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

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