User:Ledrug/bits: Difference between revisions

From Rosetta Code
Content added Content deleted
mNo edit summary
m (scrub)
Line 1: Line 1:
<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(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 && (xint)p * p <= x; 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>

Revision as of 02:29, 5 August 2011

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