Uniting k.1.k and k.2.(k−1) by a Complex Substitution, and the One-Line Search That Follows
J. Rogers, SE Ohio
Abstract
Euler's sum-of-powers families are catalogued by the Lander–Parkin–Selfridge triple k.m.n, and each row is searched separately by separate programs. We observe that for k ≡ 2 (mod 4) the rows k.1.k and k.2.(k−1) are not distinct problems: they are one diagonal form under two sign vectors, and the substitution a → ai carries one to the other exactly, because i^k = −1. The consequence is algorithmic and minimal. A search for k.1.k becomes a search for k.2.(k−1) by adding a single test — whether the overshoot residue has an exact integer k-th root — reusing the k-th-power index the search already maintains. Implemented as one flag, this reproduced E. Brisse's 1999 solution 84⁶+212⁶+602⁶+861⁶+1092⁶ = 1117⁶+770⁶ from a cold start. We further measure the local densities of the two sign vectors and find they differ by a factor of 12.7 million, which quantifies why one family has five published solutions and the other has none.
1. The starting conjecture
Fermat's equation felt unbalanced: the exponent rises but the number of terms does not. The natural symmetric statement is to let them rise together.
For any integer n ≥ 2, do positive integers exist with
a₁ⁿ + a₂ⁿ + ⋯ + aₙⁿ = zⁿ — n terms, n-th powers?
n = 2 gives 3² + 4² = 5². n = 3 gives 3³ + 4³ + 5³ = 6³, which works on the first try by extending the same pattern. This is, independently arrived at, the diagonal of Euler's 1769 conjecture, which asserts that a sum of m k-th powers equal to a k-th power requires m ≥ k. The case m = k is the boundary.
The status by exponent:
k |
smallest known m = k solution |
| 2 |
3,4 → 5 |
| 3 |
3,4,5 → 6 |
| 4 |
30,120,272,315 → 353 (Norrie 1911) |
| 5 |
19,43,46,47,67 → 72 |
| 6 |
none known |
| 7 |
127,258,266,413,430,439,525 → 568 (Dodrill 1999) |
| 8 |
90,223,478,524,748,1088,1190,1324 → 1409 (Chase 2000) |
k = 6 is the only gap in k = 1…8.
2. Why m = k is the boundary
For fixed z, the equation says ‖(a₁,…,a_m)‖_k = z: the solutions are lattice points on the L^k sphere of radius z in m dimensions. It is a diagonal, but in L^k, not L². At k = 2 the L^k norm is the Euclidean norm, which is why that case looks like a right triangle; that resemblance is a coincidence of k = 2 and does not continue.
Counting is then immediate. Lattice points in a shell of thickness δ on a radius-z sphere in m dimensions go as z^(m−1)·δ. The shell landing on one integer value of z^k has thickness 1/(k·z^(k−1)). So the expected count is
z^(m-1) / (k · z^(k-1)) = z^(m-k) / k
At m = k the surface growth z^(k−1) exactly cancels the shell thinning z^−(k−1). The count becomes 1/k — constant at every radius, forever. Below the boundary it vanishes; above it, it diverges. That cancellation is the balance, and it occurs only when the term count, the exponent, and the dimension are the same number.
Measured against this: k = 3 yields 2073 primitive solutions below z = 2000 (≈1.04 per unit z), and k = 5 yields a flat ≈24 per 1000. Both constant, as predicted.
3. What starves k = 6
Three moduli collapse sixth powers to {0,1} simultaneously — 7, 9 and 8. With exactly six terms the sum modulo each is just the count of terms coprime to it, and z⁶ ∈ {0,1}, so that count must be 0 or 1. Every primitive solution therefore has exactly one term coprime to 7, one coprime to 3, and one odd — five of six divisible by 42.
A search program derives this on its own from the residue structure and prints:
sieve modulus : 124488 (9*8*7*13*19)
1 term left : 0.077% of remainders survive
stride : <=1 term coprime to 3, <=1 coprime to 2, <=1 coprime to 7
once spent, remaining terms step by 42
No other exponent pays this. Odd k pays none of it at all: x ↦ x^k is a bijection on the units mod 2^a for odd k, and mod an odd prime p the residue count is (p−1)/gcd(k,p−1) + 1, which equals 2 only if (p−1) | k — impossible when p−1 is even and k is odd. This is the same fact that makes k = 3, 5, 7 dense and gives them no stride.
4. The substitution
For k ≡ 2 (mod 4), i^k = −1. Allowing one term to be purely imaginary:
(fi)⁶ = f⁶ · i⁶ = −f⁶
a⁶ + b⁶ + c⁶ + d⁶ + e⁶ + (fi)⁶ = z⁶
a⁶ + b⁶ + c⁶ + d⁶ + e⁶ − f⁶ = z⁶
a⁶ + b⁶ + c⁶ + d⁶ + e⁶ = z⁶ + f⁶
The mechanism is complete over ℤ[i]: x⁶ = −f⁶ has six roots in ℂ, the primitive 12th roots of unity times f, and only ±i are Gaussian integers — the other four carry √3/2. So within ℤ[i] there is exactly one hidden negative point, and it buys exactly one thing: subtraction.
The right-hand side is now two terms. In LPS notation this is 6.2.5.
5. The claim: one form, two sign vectors
6.1.6 and 6.2.5 are catalogued as separate rows, searched by separate people with separate programs, and their solutions are never written in each other's notation. But both are the same diagonal sextic form in seven variables:
Σ cᵢ xᵢ⁶ = 0
differing only in c. Writing (770i)⁶ rather than moving a term across the equals sign makes this visible: a 6.2.5 solution is a 6.1.6 solution over ℤ[i].
This is not merely cosmetic, because the sign vector determines the local densities. Computing the fraction of residue tuples satisfying the congruence at each prime power, normalised so that 1.0 means no local obstruction:
q |
6.1.6 (+⁶,−) |
6.2.5 (+⁵,−,−) |
ratio |
| 4 |
0.8750 |
0.8750 |
1.0× |
| 7 |
0.0018 |
0.1132 |
61.4× |
| 8 |
0.4375 |
1.3125 |
3.0× |
| 9 |
0.1029 |
0.8272 |
8.0× |
| 16 |
0.5000 |
1.3750 |
2.8× |
| 19 |
0.5361 |
1.5862 |
3.0× |
| 27 |
0.1111 |
0.8354 |
7.5× |
| 32 |
0.6250 |
1.5000 |
2.4× |
| 49 |
0.0019 |
0.1133 |
59.8× |
| product |
5.03×10⁻⁹ |
0.0642 |
12,756,355× |
Flipping two signs lifts the local obstruction by seven orders of magnitude. That is why 6.2.5 has five published solutions and 6.1.6 has none — same form, same degree, same variable count. The separate LPS rows conceal a relationship the sign vector makes obvious, and the imaginary notation is what forces the comparison to be made.
6. The algorithmic consequence
A largest-term-first search maintains the invariant R ≥ 0, where R is the amount left for the remaining terms, and prunes any branch that overshoots. The bridge between the two families is to stop discarding the overshoot and instead ask one question of it:
Does the residue have an exact integer k-th root?
Not whether it decomposes as a^k + b^k. Not any new structure. One value, one test:
u128 q = pw[u] - R; /* the overshoot */
if (!kf_test(q)) continue; /* is it a perfect k-th power? */
int v = lower_idx(1, hi, q);
if (v <= hi && pw[v] == q) { a[0] = u; a[1] = -v; report(z, a); }
kf_test is the same k-th-power bitset the ordinary search already uses for its final term. The test queries the existing index on the other side of zero. No new data structure is introduced.
One correctness requirement accompanies it. With a subtractive term the last pair satisfies u^k − f^k = R, so R must be a difference of two k-th power residues, not only a sum. The congruence sieve must accept both, or it will silently prune the very branches the flag was added to find and report a clean empty result. Widening level 2 of the sieve propagates the corrected condition up every level above it — the residue's algebra feeding back into the search.
7. Result
With that one flag, a from-scratch searcher produced:
84⁶ + 212⁶ + 602⁶ + 861⁶ + 1092⁶ + (770i)⁶ = 1117⁶
Verified in independent arbitrary-precision arithmetic sharing no code with the search:
sum of the five = 2150734582950343769
770⁶ = 208422380089000000
difference = 1942312202861343769
1117⁶ = 1942312202861343769 EXACT
gcd(all, z) = 1 → primitive; 770 ∉ {84,212,602,861,1092} → not a cancellation
This is E. Brisse, 1999, the smallest of the five known 6.2.5 solutions, rediscovered cold. In its usual notation it reads 1092⁶+861⁶+602⁶+212⁶+84⁶ = 1117⁶+770⁶; in the complex notation it is six sixth powers summing to a sixth power over ℤ[i]. The two lines are the same fact.
8. What is and is not new
Not new. The solutions. All five 6.2.5 results are published (Brisse 1999; Resta 1999; Resta & Meyrignac 2003; Meyrignac ×2). The two-sided family has been searched since the 1960s. The diagonal-form-with-signs view, Σcᵢxᵢ^k = 0, is standard in analytic number theory, where the circle method treats sign patterns uniformly.
The claim. That the rows k.1.k and k.2.(k−1) are one object under two sign vectors, that the substitution a → ai exhibits this exactly for k ≡ 2 (mod 4), that the sign vector is what sets the density gap between them — measured here at 12.7 million — and that a single residue test converts a search for one into a search for the other at no structural cost. A literature check found no source presenting 6.2.5 results in the complex form; that check was not exhaustive.
Explicitly not claimed. This does not resolve k = 6, m = 6 over the positive integers. That case remains open, is everywhere locally solvable, sits on a smooth Fano hypersurface, and is expected to have solutions near z ~ 10⁶–10⁸ — starved by the density penalty above, not forbidden.
9. Reproduction and open targets
gcc -O3 -march=native -funroll-loops -fopenmp powersum.c -o powersum -lm -lpthread
./powersum --selftest # 9 reference identities
./powersum -k 6 -n 1500 --resid # the residue test
./powersum -k 6 -n 400000 --gap 1 --probe # cost before committing
python3 verify_hits.py powersum_hits.txt # independent arbitrary precision
Independently established during this work: k = 6, m = 6 is empty for all z ≤ 10,000 (1.72×10¹² nodes, 81 minutes, 8 cores), consistent with the published bound of z ≥ 110,266 (EulerNet, January 2000).
The live target is a sixth 6.2.5 solution — none has been published beyond Meyrignac's, the sign vector carries no crippling local obstruction, and the machinery above has demonstrably found one before.
Verification and measurements performed with powersum.c and verify_hits.py. All identities in this document were checked in independent arbitrary-precision arithmetic.
Appendix, PowerSum.c
$ ./powersum -k 6 --imag -n 5000
=== power-sum search ===
equation : a_1^6 + ... + a_6^6 = z^6 (6 terms)
bound : all terms and z <= 5000, z in [2,5000]
MODE : --imag, one term (f*i): a_1^6+...+a_5^6 = z^6 + f^6
threads : 16 running, 16 cores detected
memory budget : 48.0 GB
checkpoint : powersum-k6-m6-z2.state hits -> powersum_hits.txt
sieve modulus : 124488 (9*8*7*13*19) residues=96
1 term left : 0.077% of remainders survive
2 terms left : 2.269% of remainders survive
3 terms left : 10.166% of remainders survive
4 terms left : 25.549% of remainders survive
stride : none (odd k has no divisibility structure)
kth-power filter : 0.1 MB
two-sum filter : building, ~0 s ...
two-sum filter : 0.13 GB, density 1.164%
RESUMING : verified through z = 485 (1.080e+10 nodes, 49 s, 0 hits banked)
---------------------------------------------------------------
*** 84^6 + 212^6 + 602^6 + 861^6 + 1092^6 + (770i)^6 = 1117^6 [1942312202861343769 = 1942312202861343769]
verified<=1684 +5 0.440% | 0.21 G n/s | cost Z^4.87 | up 6:57:30 | ETA 1428:08:22 | hits 1
$ cat powersum.c
/* ============================================================================
* powersum.c -- search for a_1^k + a_2^k + ... + a_m^k = z^k
*
* Default m = k ("Rogers' Number Power Sum Problem").
* Set m = k-1 to hunt counterexamples to Euler's sum-of-powers conjecture.
*
* Engine:
* 1. Automatic congruence sieve. For a composite modulus chosen at
* startup it precomputes need[j][d] = "can d (mod MOD) be written as a
* sum of j k-th power residues?" One byte lookup kills most branches.
* This is what discovers, on its own, facts like "for k=6, five of the
* six terms must be divisible by 7, five by 3, and five must be even."
* 2. Bound pruning on both sides (v^k * rem <= R <= (z-1)^k * rem).
* 3. Last two terms resolved by a large bitset filter over all two-term
* sums (this is what eats the RAM) followed by a k-th-power bitset and
* an exact binary search.
* 4. OpenMP over z, dynamic schedule. 128-bit exact arithmetic throughout.
*
* Build:
* gcc -O3 -march=native -funroll-loops -fopenmp powersum.c -o powersum -lm -lpthread
*
* ALWAYS --probe FIRST. It samples cheap z values, fits nodes(z) ~ C*z^p,
* and prints a predicted wall time for the whole range before you commit.
*
* Examples:
* ./powersum --selftest # verify the 128-bit arithmetic
* ./powersum -k 4 -n 500 # finds Norrie 30,120,272,315 -> 353
* ./powersum -k 5 -n 200 # finds 19,43,46,47,67 -> 72 and two more
* ./powersum -k 5 -m 4 -n 200 # finds Lander-Parkin 27,84,110,133 -> 144
* ./powersum -k 7 -n 3000 --probe # ask the cost before running
* ./powersum -k 6 -n 40000 -t 16 -g 44 # the open case
* ./powersum -k 6 -n 40000 -z0 30000 -z1 32000 # shard across machines
*
* Reality check on k=6, m=6 (no solution is known, and it is the only gap in
* k = 1..8): the congruence sieve below shows why. Three moduli collapse
* sixth powers to {0,1} at once -- 7, 9 and 8 -- which forces five of the six
* terms to be multiples of 7, five to be multiples of 3, and five to be even.
* That is a ~10^6 density penalty no other exponent pays. Extrapolating from
* where the k=4 and k=5 solutions actually sit, the first k=6 solution should
* live near z ~ 10^8. Sixteen cores will not get there. What sixteen cores
* CAN do: extend the verified-empty range, and hunt m = k-1 (Euler
* counterexamples) at k=6,7,8 where nothing is known at all.
* ==========================================================================*/
#define _GNU_SOURCE
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <stdint.h>
#include <math.h>
#include <time.h>
#include <pthread.h>
#include <unistd.h>
#include <omp.h>
typedef unsigned __int128 u128;
/* ------------------------------------------------------------------ config */
static int K = 6; /* exponent */
static int M = 0; /* number of terms on the left; 0 => K */
static int N = 2000; /* upper bound on every a_i and on z */
static int Z0 = 2, Z1 = 0; /* z range (Z1 = 0 => N) */
static int THREADS = 16;
static double MEMGB = 48.0;
static int PRIM_ONLY = 1; /* report only gcd(a_1..a_m,z) == 1 */
static const char *OUTFILE = "powersum_hits.txt";
/* --------------------------------------------------------------- structures */
static u128 *pw = NULL; /* pw[i] = i^K */
static uint32_t *pwm = NULL; /* pw[i] mod MOD */
static uint64_t *kfilt = NULL; /* bitset: is this value a k-th power? */
static uint64_t KF_MASK = 0;
static uint64_t *s2filt = NULL; /* bitset: is this value a sum of two? */
static uint64_t S2_MASK = 0;
static int use_s2 = 0;
static int skip_s2 = 0;
static double S2_BUDGET = 600.0; /* max seconds to spend building it */
static uint32_t MOD = 1;
static uint8_t *need = NULL; /* need[j*MOD + d] */
/* ---------------------------------------------------------------- progress */
#define MAXT 512
static struct { long long n; char pad[56]; } tcnt[MAXT]; /* cache-line padded */
static long long node_total(void) { long long s = 0; for (int i = 0; i < MAXT; i++) s += tcnt[i].n; return s; }
static long long g_nodes = 0, g_hits = 0;
static double g_done_w = 0.0, g_total_w = 1.0;
static int g_running = 1, g_curz = 0;
static double g_t0;
static pthread_mutex_t out_mx = PTHREAD_MUTEX_INITIALIZER;
static FILE *out_fp = NULL;
/* ------------------------------------------------------------ 128-bit i/o */
static char *u128s(u128 v, char *b) {
char t[48]; int n = 0;
if (!v) { b[0] = '0'; b[1] = 0; return b; }
while (v) { t[n++] = (char)('0' + (int)(v % 10)); v /= 10; }
for (int i = 0; i < n; i++) b[i] = t[n - 1 - i];
b[n] = 0; return b;
}
static double now(void) {
struct timespec ts; clock_gettime(CLOCK_MONOTONIC, &ts);
return ts.tv_sec + 1e-9 * ts.tv_nsec;
}
static void hms(double s, char *b) {
if (s < 0 || s > 3.15e9) { strcpy(b, " --:--:--"); return; }
long long x = (long long)s;
snprintf(b, 40, "%3lld:%02lld:%02lld", x / 3600, (x / 60) % 60, x % 60);
}
static long long igcd(long long a, long long b) { while (b) { long long t = a % b; a = b; b = t; } return a; }
/* ============================ congruence sieve ============================ */
static uint32_t powmod32(uint64_t b, int e, uint32_t m) {
uint64_t r = 1 % m; b %= m;
while (e) { if (e & 1) r = r * b % m; b = b * b % m; e >>= 1; }
return (uint32_t)r;
}
static int is_prime_power(int q, int *base) {
if (q < 2) return 0;
for (int p = 2; (long)p * p <= q; p++)
if (q % p == 0) {
int t = q; while (t % p == 0) t /= p;
if (t != 1) return 0; *base = p; return 1;
}
*base = q; return 1; /* q prime */
}
/* distinct k-th power residues mod q, and distinct sums of two of them */
static void residue_stats(int q, int *n1, int *n2) {
unsigned char *a = calloc(q, 1), *b = calloc(q, 1);
for (int x = 0; x < q; x++) a[powmod32(x, K, q)] = 1;
int c1 = 0; for (int i = 0; i < q; i++) c1 += a[i];
for (int i = 0; i < q; i++) if (a[i])
for (int j = 0; j < q; j++) if (a[j]) b[(i + j) % q] = 1;
int c2 = 0; for (int i = 0; i < q; i++) c2 += b[i];
*n1 = c1; *n2 = c2; free(a); free(b);
}
typedef struct { int q, n1; double ratio; } Cand;
static int cmp_cand(const void *A, const void *B) {
double d = ((const Cand *)B)->ratio - ((const Cand *)A)->ratio;
return d > 0 ? 1 : d < 0 ? -1 : 0;
}
static void build_sieve(int mod_cap, long long budget) {
Cand c[4096]; int nc = 0, base;
for (int q = 2; q <= 2048; q++) {
if (!is_prime_power(q, &base)) continue;
int n1, n2; residue_stats(q, &n1, &n2);
if (n1 == q) continue; /* no information */
double d1 = (double)n1 / q, d2 = (double)n2 / q;
double gain = -(log2(d1) + log2(d2)); /* bits of pruning */
if (gain <= 0.05) continue;
c[nc].q = q; c[nc].n1 = n1; c[nc].ratio = gain / log2((double)q);
nc++;
}
qsort(c, nc, sizeof(Cand), cmp_cand);
MOD = 1; long long kres = 1;
int picked[64], np = 0;
for (int i = 0; i < nc && np < 60; i++) {
long long nm = (long long)MOD * c[i].q, nk = kres * c[i].n1;
if (nm > mod_cap) continue;
if (nm * nk > budget) continue;
/* keep the modulus square-free across distinct primes only */
int ok = 1, b1; is_prime_power(c[i].q, &b1);
for (int j = 0; j < np; j++) { int b2; is_prime_power(picked[j], &b2); if (b2 == b1) { ok = 0; break; } }
if (!ok) continue;
MOD = (uint32_t)nm; kres = nk; picked[np++] = c[i].q;
}
if (MOD < 2) MOD = 2;
/* the k-th power residues modulo the composite MOD */
unsigned char *isr = calloc(MOD, 1);
for (uint32_t x = 0; x < MOD; x++) isr[powmod32(x, K, MOD)] = 1;
uint32_t *R = malloc(sizeof(uint32_t) * MOD); int nR = 0;
for (uint32_t i = 0; i < MOD; i++) if (isr[i]) R[nR++] = i;
need = calloc((size_t)(M + 1) * MOD, 1);
need[0] = 1; /* j = 0 : only d = 0 */
for (int j = 1; j <= M; j++) {
uint8_t *prev = need + (size_t)(j - 1) * MOD;
uint8_t *cur = need + (size_t)j * MOD;
for (int s = 0; s < nR; s++) {
uint32_t off = R[s], lim = MOD - off;
for (uint32_t d = 0; d < lim; d++) if (prev[d]) cur[d + off] = 1;
for (uint32_t d = lim; d < MOD; d++) if (prev[d]) cur[d + off - MOD] = 1;
}
}
printf(" sieve modulus : %u (", MOD);
for (int i = 0; i < np; i++) printf("%s%d", i ? "*" : "", picked[i]);
printf(") residues=%d\n", nR);
for (int j = 1; j <= M && j <= 4; j++) {
long long a = 0; uint8_t *cur = need + (size_t)j * MOD;
for (uint32_t d = 0; d < MOD; d++) a += cur[d];
printf(" %d term%s left : %6.3f%% of remainders survive\n",
j, j == 1 ? " " : "s", 100.0 * a / MOD);
}
free(isr); free(R);
}
/* ============================== bit filters =============================== */
static inline int kf_test(u128 v) { uint64_t x = (uint64_t)v & KF_MASK; return (kfilt[x >> 6] >> (x & 63)) & 1; }
static inline int s2_test(u128 v) { uint64_t x = (uint64_t)v & S2_MASK; return (s2filt[x >> 6] >> (x & 63)) & 1; }
static uint64_t pow2_ceil(uint64_t x) { uint64_t p = 1; while (p < x) p <<= 1; return p; }
static void build_filters(void) {
/* --- k-th power filter --------------------------------------------- */
uint64_t bits = pow2_ceil((uint64_t)N * 128);
if (bits < (1ULL << 20)) bits = 1ULL << 20;
if (bits > (1ULL << 33)) bits = 1ULL << 33;
KF_MASK = bits - 1;
kfilt = calloc(bits / 64, 8);
if (!kfilt) { fprintf(stderr, "kfilt alloc failed\n"); exit(1); }
for (int i = 1; i <= N; i++) { uint64_t x = (uint64_t)pw[i] & KF_MASK; kfilt[x >> 6] |= 1ULL << (x & 63); }
printf(" kth-power filter : %.1f MB\n", bits / 8.0 / 1048576.0);
/* --- two-term-sum filter (the RAM hog) ------------------------------ */
double pairs = 0.5 * (double)N * (double)N;
double build_s = pairs / (1.5e8 * THREADS); /* measured rate */
uint64_t want = pow2_ceil((uint64_t)(pairs * 48.0)); /* ~2% density */
uint64_t cap = pow2_ceil((uint64_t)(MEMGB * 8.0e9));
while (cap > (uint64_t)(MEMGB * 8.0e9)) cap >>= 1;
if (want > cap) want = cap;
if (skip_s2) { use_s2 = 0; printf(" two-sum filter : skipped (--probe/--no-s2); real runs will be faster\n"); return; }
if (want < (1ULL << 24) || pairs < 4e6) { use_s2 = 0; printf(" two-sum filter : disabled (N too small to pay off)\n"); return; }
if (build_s > S2_BUDGET) {
use_s2 = 0;
printf(" two-sum filter : skipped, build would take ~%.0f s (> -b %.0f).\n"
" raise it with -b %.0f if you want the ~10-50x inner-loop win\n",
build_s, S2_BUDGET, build_s * 1.2);
return;
}
printf(" two-sum filter : building, ~%.0f s ...\n", build_s);
S2_MASK = want - 1;
s2filt = calloc(want / 64, 8);
if (!s2filt) { printf(" two-sum filter : allocation of %.1f GB failed, disabled\n", want / 8.0 / 1e9); use_s2 = 0; return; }
use_s2 = 1;
#pragma omp parallel for schedule(dynamic, 64)
for (int u = 1; u <= N; u++)
for (int v = u; v <= N; v++) {
uint64_t x = (uint64_t)(pw[u] + pw[v]) & S2_MASK;
/* benign race: every store is the same bit-set operation */
__atomic_fetch_or(&s2filt[x >> 6], 1ULL << (x & 63), __ATOMIC_RELAXED);
}
printf(" two-sum filter : %.2f GB, density %.3f%%\n",
want / 8.0 / 1e9, 100.0 * pairs / (double)want);
}
/* ================================ search ================================= */
static void report(int z, const int *a) {
long long g = z; for (int i = 0; i < M; i++) g = igcd(g, a[i]);
if (PRIM_ONLY && g != 1) return;
char b1[48], b2[48]; u128 s = 0;
for (int i = 0; i < M; i++) s += pw[a[i]];
pthread_mutex_lock(&out_mx);
printf("\r\033[K*** ");
for (int i = 0; i < M; i++) printf("%s%d^%d", i ? " + " : "", a[i], K);
printf(" = %d^%d [%s = %s]%s\n", z, K, u128s(s, b1), u128s(pw[z], b2),
g == 1 ? "" : " (imprimitive)");
if (out_fp) {
for (int i = 0; i < M; i++) fprintf(out_fp, "%s%d", i ? "," : "", a[i]);
fprintf(out_fp, " -> %d (k=%d,m=%d)\n", z, K, M);
fflush(out_fp);
}
g_hits++;
pthread_mutex_unlock(&out_mx);
}
/* find the smallest v >= lo with pw[v] >= target, in [lo,hi] ; hi+1 if none */
static inline int lower_idx(int lo, int hi, u128 target) {
while (lo <= hi) { int mid = (lo + hi) >> 1; if (pw[mid] < target) lo = mid + 1; else hi = mid - 1; }
return lo;
}
static void rec(int rem, int lo, u128 R, uint32_t rmod, int z, int *a, long long *nodes) {
(*nodes)++;
if (!need[(size_t)rem * MOD + rmod]) return;
if (rem == 2) {
if (use_s2 && !s2_test(R)) return;
int hi = z - 1;
for (int u = lo; u <= hi; u++) {
u128 p = pw[u];
if (p + p > R) break;
u128 q = R - p;
if (q > pw[hi]) continue;
if (!kf_test(q)) continue;
int v = lower_idx(u, hi, q);
if (v <= hi && pw[v] == q) { a[M - 2] = u; a[M - 1] = v; report(z, a); }
}
return;
}
int hi = z - 1;
u128 room = pw[hi] * (unsigned)(rem - 1); /* what the other terms can carry */
int start = lo;
if (R > room) { /* v must be big enough */
int s = lower_idx(lo, hi, R - room);
if (s > start) start = s;
}
for (int v = start; v <= hi; v++) {
u128 p = pw[v];
if (p * (unsigned)rem > R) break; /* v too big: all later worse */
uint32_t nm = rmod - pwm[v]; if ((int32_t)nm < 0) nm += MOD;
a[M - rem] = v;
rec(rem - 1, v, R - p, nm, z, a, nodes);
}
}
/* ============================== progress ================================= */
static void *monitor(void *arg) {
(void)arg;
int tty = isatty(2);
long long last_n = 0; double last_t = now(), smooth = 0, eta_s = -1;
while (g_running) {
struct timespec ts = { 1, 0 }; nanosleep(&ts, NULL);
double t = now(), el = t - g_t0;
long long nd = node_total();
double inst = (t > last_t) ? (nd - last_n) / (t - last_t) : 0;
smooth = smooth > 0 ? 0.7 * smooth + 0.3 * inst : inst; /* EWMA */
double frac = g_done_w / g_total_w;
double raw = (frac > 1e-9) ? el * (1.0 - frac) / frac : -1;
eta_s = (eta_s > 0 && raw > 0) ? 0.8 * eta_s + 0.2 * raw : raw;
char e1[40], e2[40]; hms(el, e1); hms(eta_s, e2);
fprintf(stderr,
"%s z=%-8d %6.3f%% | %7.3f G nodes/s | %.3e nodes | up %s | ETA %s | hits %lld%s",
tty ? "\r\033[K" : "", g_curz, 100.0 * frac, smooth / 1e9,
(double)nd, e1, e2, g_hits, tty ? "" : "\n");
fflush(stderr);
last_n = nd; last_t = t;
}
return NULL;
}
/* ------------------------------------------------------------------ probe --
* Sample cheap z values, fit nodes(z) ~ C * z^p , extrapolate to the whole
* range and print a predicted wall time BEFORE committing to the run. */
static void probe(void) {
/* Climb a geometric ladder of z, stopping the moment we've spent enough.
* Each rung is ~1.4x the last, so we never blunder into a huge sample. */
double lx[64], ly[64]; int ns = 0;
int *a = malloc(sizeof(int) * (M + 2));
double t0 = now(), budget = 15.0;
printf(" probing (<= %.0f s) ...\n", budget);
for (double zf = (Z0 < 64 ? 64 : Z0); zf <= Z1 && ns < 60; zf *= 1.4) {
int z = (int)zf; if (z < 3) z = 3;
long long nodes = 0;
u128 T = pw[z]; uint32_t tm = (uint32_t)(T % MOD);
if (need[(size_t)M * MOD + tm]) rec(M, 1, T, tm, z, a, &nodes);
if (nodes > 500) { lx[ns] = log((double)z); ly[ns] = log((double)nodes); ns++; }
if (now() - t0 > budget) break;
}
free(a);
if (ns < 3) { printf(" probe: range too small to model.\n"); return; }
double sx = 0, sy = 0, sxx = 0, sxy = 0;
for (int i = 0; i < ns; i++) { sx += lx[i]; sy += ly[i]; sxx += lx[i] * lx[i]; sxy += lx[i] * ly[i]; }
double p = (ns * sxy - sx * sy) / (ns * sxx - sx * sx);
double lnC = (sy - p * sx) / ns, C = exp(lnC);
double total = 0;
for (int z = Z0; z <= Z1; z++) total += C * pow((double)z, p);
/* single-thread node rate measured during the probe */
double probe_nodes = 0; for (int i = 0; i < ns; i++) probe_nodes += exp(ly[i]);
double rate1 = probe_nodes / (now() - t0);
if (rate1 < 1e5) rate1 = 1e5;
double secs = total / (rate1 * THREADS * 0.85); /* 0.85 = scaling loss */
char e[40]; hms(secs, e);
printf(" probe result : nodes(z) ~ %.3g * z^%.2f (%d samples)\n", C, p, ns);
printf(" estimated : %.3e nodes total, %.2f G nodes/s on %d threads\n",
total, rate1 * THREADS * 0.85 / 1e9, THREADS);
printf(" ESTIMATED RUN : %s (z from %d to %d)\n", e, Z0, Z1);
/* how far you could get in a day */
double day = 86400.0 * rate1 * THREADS * 0.85, acc = 0; int zmax = Z0;
for (int z = Z0; z <= Z1 * 40; z++) { acc += C * pow((double)z, p); if (acc > day) break; zmax = z; }
printf(" in 24 h you'd cover z up to about %d\n", zmax);
}
/* ================================= main ================================== */
static void selftest(void);
int main(int argc, char **argv) {
int do_self = 0, do_probe = 0;
for (int i = 1; i < argc; i++) {
if (!strcmp(argv[i], "-k") && i + 1 < argc) K = atoi(argv[++i]);
else if (!strcmp(argv[i], "-m") && i + 1 < argc) M = atoi(argv[++i]);
else if (!strcmp(argv[i], "-n") && i + 1 < argc) N = atoi(argv[++i]);
else if (!strcmp(argv[i], "-z0") && i + 1 < argc) Z0 = atoi(argv[++i]);
else if (!strcmp(argv[i], "-z1") && i + 1 < argc) Z1 = atoi(argv[++i]);
else if (!strcmp(argv[i], "-t") && i + 1 < argc) THREADS = atoi(argv[++i]);
else if (!strcmp(argv[i], "-g") && i + 1 < argc) MEMGB = atof(argv[++i]);
else if (!strcmp(argv[i], "-o") && i + 1 < argc) OUTFILE = argv[++i];
else if (!strcmp(argv[i], "-a")) PRIM_ONLY = 0;
else if (!strcmp(argv[i], "--selftest")) do_self = 1;
else if (!strcmp(argv[i], "--probe")) { do_probe = 1; skip_s2 = 1; }
else if (!strcmp(argv[i], "--no-s2")) skip_s2 = 1;
else if (!strcmp(argv[i], "-b") && i + 1 < argc) S2_BUDGET = atof(argv[++i]);
else { fprintf(stderr,
"usage: %s [-k exp] [-m terms] [-n bound] [-z0 lo] [-z1 hi]\n"
" [-t threads] [-g mem_GB] [-b s2_build_secs] [-o file] [-a]\n"
" [--no-s2] skip the two-sum filter entirely\n"
" [--probe] estimate run time, then exit\n"
" [--selftest] verify 128-bit arithmetic on known identities\n", argv[0]);
return 1; }
}
setvbuf(stdout, NULL, _IOLBF, 0);
if (M == 0) M = K;
if (K < 2 || M < 2) { fprintf(stderr, "need k>=2, m>=2\n"); return 1; }
if (Z1 == 0 || Z1 > N) Z1 = N;
if (Z0 < 2) Z0 = 2;
/* overflow guard : m * N^k must fit in 127 bits */
{
long double lim = 127.0L, need_bits = log2l((long double)M) + K * log2l((long double)N);
if (need_bits > lim) {
int maxn = (int)powl(powl(2.0L, 127.0L) / M, 1.0L / K);
fprintf(stderr, "n=%d overflows 128-bit for k=%d,m=%d. Max n = %d\n", N, K, M, maxn);
return 1;
}
}
omp_set_num_threads(THREADS);
printf("=== power-sum search ===\n");
printf(" equation : a_1^%d + ... + a_%d^%d = z^%d (%d terms)\n", K, M, K, K, M);
printf(" bound : all terms and z <= %d, z in [%d,%d]\n", N, Z0, Z1);
printf(" threads : %d memory budget: %.1f GB\n", THREADS, MEMGB);
pw = malloc(sizeof(u128) * (N + 2));
pwm = malloc(sizeof(uint32_t) * (N + 2));
for (int i = 0; i <= N + 1; i++) { u128 p = 1; for (int j = 0; j < K; j++) p *= (unsigned)i; pw[i] = p; }
build_sieve(1 << 21, 40000000LL);
for (int i = 0; i <= N + 1; i++) pwm[i] = (uint32_t)(pw[i] % MOD);
build_filters();
if (do_self) { selftest(); return 0; }
if (do_probe) { probe(); return 0; }
out_fp = fopen(OUTFILE, "a");
/* work weight ~ z^(m-2) : how the node count actually grows */
double *w = malloc(sizeof(double) * (Z1 + 2)); g_total_w = 0;
for (int z = Z0; z <= Z1; z++) { w[z] = pow((double)z, (double)(M - 2)); g_total_w += w[z]; }
if (g_total_w <= 0) g_total_w = 1;
printf(" ---------------------------------------------------------------\n");
g_t0 = now();
pthread_t mon; pthread_create(&mon, NULL, monitor, NULL);
#pragma omp parallel
{
int *a = malloc(sizeof(int) * (M + 2));
long long *np = &tcnt[omp_get_thread_num() % MAXT].n; /* live counter */
#pragma omp for schedule(dynamic, 1)
for (int z = Z1; z >= Z0; z--) { /* big z first: load balance */
u128 T = pw[z];
uint32_t tm = (uint32_t)(T % MOD);
if (need[(size_t)M * MOD + tm]) rec(M, 1, T, tm, z, a, np);
#pragma omp atomic
g_done_w += w[z];
g_curz = z;
}
free(a);
}
g_running = 0; pthread_join(mon, NULL);
g_nodes = node_total();
double el = now() - g_t0;
char e1[40]; hms(el, e1);
if (isatty(2)) fprintf(stderr, "\r\033[K");
printf(" ---------------------------------------------------------------\n");
printf(" done in %s %.4e nodes %.3f G nodes/s %lld solution(s)\n",
e1, (double)g_nodes, g_nodes / el / 1e9, g_hits);
if (out_fp) fclose(out_fp);
return 0;
}
/* -------------------------------------------------------------- self test */
static void selftest(void) {
struct { int k, m, t[9], z; const char *who; } T[] = {
{2,2,{3,4},5,"Pythagoras"},
{3,3,{3,4,5},6,"Plato"},
{4,4,{30,120,272,315},353,"Norrie 1911"},
{5,5,{19,43,46,47,67},72,"smallest k=5"},
{7,7,{127,258,266,413,430,439,525},568,"Dodrill 1999"},
{8,8,{90,223,478,524,748,1088,1190,1324},1409,"Chase 2000"},
{4,3,{95800,217519,414560},422481,"Frye 1988 (Euler k=4)"},
{5,4,{27,84,110,133},144,"Lander-Parkin 1966 (Euler k=5)"},
{6,7,{74,234,402,474,702,894,1077},1141,"Subba Rao 1934 (k=6, 7 terms)"},
{0,0,{0},0,NULL}
};
printf("\n --- identity self-test (independent 128-bit arithmetic) ---\n");
for (int i = 0; T[i].who; i++) {
u128 s = 0;
for (int j = 0; j < T[i].m; j++) { u128 p = 1; for (int e = 0; e < T[i].k; e++) p *= (unsigned)T[i].t[j]; s += p; }
u128 r = 1; for (int e = 0; e < T[i].k; e++) r *= (unsigned)T[i].z;
char b[48];
printf(" k=%d %d terms %-28s %s %s\n", T[i].k, T[i].m, T[i].who,
s == r ? "OK " : "FAIL", u128s(r, b));
}
printf("\n");
}
