The same OEIS page has:
Catalin Francu, C++ program
The program is made publicly available by OEIS:
From Catalin Francu
C++ source code that generates A005228.
Runs in (approx.) O( n + sqrt(n) + sqrt(sqrt(n)) + sqrt(sqrt(sqrt(n)))...)
to generate the first n elements. Tested up to n = 1,000,000,000.
Input is n, output is a(n). Needs quasi-constant memory.
%o A005228 #include <stdio.h>
typedef unsigned long long uint64;
typedef struct {
uint64 value;
uint64 inc;
} pair;
pair series[100];
void init_all_series() {
pair initial_pair = { 7, 5 }; // pregenerate the first three terms
for (int i = 0; i < 99; series[i++] = initial_pair);
}
void generate(int k);
inline bool test(int k, unsigned number) {
switch (number) {
case 1: case 3: case 7: return 1;
case 2: case 4: case 5: case 6: return 0;
default:
while (series[k].value < number) generate(k);
return (series[k].value == number);
}
}
inline void generate(int k) {
series[k].value += series[k].inc;
if (test(k+1, ++series[k].inc)) series[k].inc++;
}
int main(void) {
int x;
init_all_series();
scanf("%d", &x); // Enter the index
switch (x) {
case 1: puts("1"); break;
case 2: puts("3"); break;
default:
for (int k = 4; k <= x; k++) generate(0);
printf("%llu", series[0].value); // print the x'th term
}
return 0;
}This has the same space consumption as my new version >>36, but linear runtime, so it won't get to R(10^20) even with GMP. The space consumption is achieved with a different method that I still need to have a think on. But this already means that >>36 is algorithmically superior to both the C and C++ versions from OEIS. A bit disappointing.