171 lines
2.9 KiB
C
171 lines
2.9 KiB
C
#include <inttypes.h>
|
|
#include <stdio.h>
|
|
#include <stdlib.h>
|
|
#include <string.h>
|
|
#include <assert.h>
|
|
|
|
typedef uint32_t pint;
|
|
typedef uint64_t xint;
|
|
typedef unsigned int uint;
|
|
#define PRIuPINT PRIu32 /* printf macro for pint */
|
|
#define PRIuXINT PRIu64 /* printf macro for xint */
|
|
#define MAX_FACTORS 63 /* because 2^64 is too large for xint */
|
|
|
|
uint8_t *pbits;
|
|
|
|
#define MAX_PRIME (~(pint)0)
|
|
#define MAX_PRIME_SQ 65535U
|
|
#define PBITS (MAX_PRIME / 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()
|
|
{
|
|
FILE *fp;
|
|
pint s, tgt = 4;
|
|
|
|
if (!(pbits = malloc(PBITS))) {
|
|
perror("malloc");
|
|
exit(1);
|
|
}
|
|
|
|
if ((fp = fopen("primebits", "r"))) {
|
|
fread(pbits, 1, PBITS, fp);
|
|
fclose(fp);
|
|
return;
|
|
}
|
|
|
|
memset(pbits, 255, PBITS);
|
|
for (s = 7; s <= MAX_PRIME_SQ; s = next_prime(s)) {
|
|
if (s > tgt) {
|
|
tgt *= 2;
|
|
fprintf(stderr, "sieve %"PRIuPINT"\n", s);
|
|
}
|
|
sieve(s);
|
|
}
|
|
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)
|
|
{
|
|
unsigned char b[8];
|
|
off_t ofs[8];
|
|
int i, q;
|
|
|
|
for (i = 0; i < 8; i++) {
|
|
q = rem_num[i] * p;
|
|
b[i] = ~bit_pos[q % 30];
|
|
ofs[i] = q / 30;
|
|
}
|
|
|
|
for (q = ofs[1], i = 7; i; i--)
|
|
ofs[i] -= ofs[i-1];
|
|
|
|
for (ofs[0] = p, i = 1; i < 8; i++)
|
|
ofs[0] -= ofs[i];
|
|
|
|
for (i = 1; q < PBITS; q += ofs[i = (i + 1) & 7])
|
|
pbits[q] &= b[i];
|
|
}
|
|
|
|
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) {
|
|
if (++addr >= PBITS) return 0;
|
|
bits = 1;
|
|
rem = 0;
|
|
}
|
|
if (addr >= PBITS) return 0;
|
|
while (!(pbits[addr] & bits)) {
|
|
rem++;
|
|
bits <<= 1;
|
|
}
|
|
return p = addr * 30 + rem_num[rem];
|
|
}
|
|
|
|
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;
|
|
|
|
/* check small primes: not strictly necessary */
|
|
if (n <= MAX_PRIME && is_prime(n)) {
|
|
f[0] = n;
|
|
return 1;
|
|
}
|
|
|
|
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[MAX_FACTORS], po;
|
|
|
|
init_primes();
|
|
|
|
for (p = 1; p < 64; p++) {
|
|
po = (1LLU << p) - 1;
|
|
printf("2^%"PRIuPINT" - 1 = %"PRIuXINT, p, po);
|
|
fflush(stdout);
|
|
if ((len = decompose(po, f)) > 1)
|
|
for (i = 0; i < len; i++)
|
|
printf(" %c %"PRIuXINT, i?'x':'=', f[i]);
|
|
putchar('\n');
|
|
}
|
|
|
|
return 0;
|
|
}
|