Project Euler Problem 558

Solution

Answer: 226754889

Let $r$ be the real root of

$$x^3=x^2+1.$$

Then

$$r^3-r^2=1,$$

so powers of $r$ satisfy the linear relation

$$r^k=r^{k-1}+r^{k-3}.$$

This induces a canonical “digit system” analogous to Zeckendorf representations for Fibonacci numbers.

The unique representation

$$n=\sum_k b_k r^k$$

with

$$b_k\in{0,1}, \qquad b_k+b_{k+1}+b_{k+2}\le 1$$

can be manipulated entirely by local rewrite rules derived from

$$r^k=r^{k-1}+r^{k-3}.$$

A very efficient approach is:

store the representation as a bitset,
implement insertion of a digit by a small rewrite automaton,
compute squares incrementally using

$$n^2=(n-1)^2+(2n-1),$$

so we only perform additions in the irrational-base representation.

The following optimized C++ implementation computes

$$S(5,000,000)=\sum_{n=1}^{5,000,000} w(n^2)$$

in a few seconds:

#include <bits/stdc++.h>
using namespace std;

const int BITS = 2048;
const int WORDS = BITS / 64;
const int SHIFT = 256;

struct BS {
    uint64_t w[WORDS]{};
};

inline bool get(const BS& b, int x) {
    return (b.w[x >> 6] >> (x & 63)) & 1ULL;
}

inline void flip(BS& b, int x) {
    b.w[x >> 6] ^= 1ULL << (x & 63);
}

inline void setb(BS& b, int x) {
    b.w[x >> 6] |= 1ULL << (x & 63);
}

inline int bitcount(const BS& b) {
    int s = 0;
    for (int i = 0; i < WORDS; i++)
        s += __builtin_popcountll(b.w[i]);
    return s;
}

inline void dfs(BS& bits, vector<int>& st) {
    while (!st.empty()) {
        int x = st.back();
        st.pop_back();

        if (get(bits, x + 2)) {
            flip(bits, x + 2);
            st.push_back(x + 3);
            continue;
        }

        if (get(bits, x - 2)) {
            flip(bits, x - 2);
            st.push_back(x + 1);
            continue;
        }

        if (get(bits, x - 1)) {
            flip(bits, x - 1);
            st.push_back(x + 1);
            st.push_back(x - 4);
            continue;
        }

        if (get(bits, x + 1)) {
            flip(bits, x + 1);
            st.push_back(x + 2);
            st.push_back(x - 3);
            continue;
        }

        if (get(bits, x)) {
            flip(bits, x);
            st.push_back(x + 1);
            st.push_back(x - 2);
            st.push_back(x - 7);
            continue;
        }

        setb(bits, x);
    }
}

int main() {
    BS f{}, g{};
    long long ans = 0;

    vector<int> sf, sg;
    sf.reserve(10);
    sg.reserve(50);

    const int N = 5000000;

    for (int i = 1; i <= N; i++) {

        // f <- f + 2  (current odd increment 2i-1)
        sf.clear();
        sf.push_back(SHIFT);
        if (i != 1) sf.push_back(SHIFT);
        dfs(f, sf);

        // g <- g + f
        sg.clear();

        for (int j = 0; j < WORDS; j++) {
            uint64_t t = f.w[j];

            while (t) {
                int b = __builtin_ctzll(t);
                sg.push_back(j * 64 + b);
                t &= t - 1;
            }
        }

        dfs(g, sg);

        ans += bitcount(g);
    }

    cout << ans << endl;
}

The implementation reproduces the given checks:

$$S(10)=61,\qquad S(1000)=19403.$$

Running it for $N=5{,}000{,}000$ gives:

$$S(5{,}000{,}000)=226754889.$$

Answer: 226754889