Computing a HAP basis
From Polymath Wiki
This is some code I wrote to compute an alternative basis (the HAP basis, or basis of homogeneous arithmetic progressions) for a given sequence. It prints the list of integers n such that {n, 2n, 3n, ...} is a basis vector to stdout.
It hopefully isn't that hard to use, although it's not incredibly well-commented. Replace the 1124-sequence I gave in the declaration of "sequence" with your sequence, and change "#define SEQSIZE foo" to match the length of your sequence. Compile, and run.
/* Computes a HAP basis for a given sequence * Author: Harrison Brown * Date: 1/14/09 * I hereby release this work into the public domain */ // This is a hacky algorithm that could almost certainly be made faster, but right now the sequences are small enough that it doesn't really matter. #include <stdio.h> #include <math.h> #define SEQSIZE 1125 // size of the input sequence plus 1 int main(void){ const char sequence[SEQSIZE] = {0, +1, -1, +1, +1, -1, -1, -1, -1, +1, +1, -1, +1, +1, +1, -1, -1, +1, -1, +1, +1, -1, -1, -1, +1, -1, +1, +1, -1, +1, -1, -1, +1, +1, -1, +1, -1, -1, +1, -1, -1, +1, +1, -1, +1, +1, -1, -1, -1, +1, +1, +1, -1, -1, +1, -1, +1, -1, -1, +1, +1, +1, +1, -1, -1, +1, -1, -1, +1, +1, -1, -1, +1, +1, +1, -1, -1, +1, +1, -1, +1, -1, -1, +1, -1, -1, +1, +1, -1, +1, -1, -1, +1, -1, +1, +1, +1, -1, -1, +1, -1, +1, -1, -1, +1, +1, -1, -1, -1, +1, +1, -1, +1, +1, +1, -1, -1, -1, +1, -1, -1, +1, -1, +1, +1, +1, +1, -1, -1, -1, -1, +1, +1, -1, +1, +1, +1, -1, -1, +1, -1, -1, +1, -1, -1, +1, -1, +1, -1, +1, +1, -1, +1, +1, +1, -1, -1, -1, -1, +1, +1, +1, +1, -1, -1, -1, +1, +1, -1, +1, -1, -1, +1, -1, +1, +1, -1, -1, +1, -1, +1, +1, -1, +1, -1, +1, -1, +1, +1, -1, -1, -1, +1, +1, -1, +1, +1, -1, -1, +1, -1, -1, +1, -1, -1, +1, +1, +1, +1, -1, +1, -1, -1, +1, -1, -1, +1, +1, -1, -1, +1, -1, +1, +1, -1, -1, -1, +1, -1, +1, +1, -1, +1, +1, +1, -1, -1, +1, +1, -1, -1, +1, -1, -1, +1, -1, +1, +1, -1, -1, +1, -1, +1, +1, -1, +1, +1, -1, -1, +1, -1, -1, +1, -1, +1, +1, -1, -1, +1, +1, -1, +1, -1, +1, -1, -1, +1, -1, +1, +1, +1, -1, -1, +1, -1, +1, +1, -1, -1, -1, -1, +1, +1, -1, -1, +1, -1, +1, +1, -1, +1, +1, -1, -1, +1, -1, +1, +1, -1, -1, +1, +1, -1, +1, -1, -1, +1, -1, +1, -1, -1, +1, +1, -1, -1, +1, -1, +1, +1, +1, -1, +1, -1, -1, +1, -1, +1, -1, -1, +1, +1, +1, +1, -1, -1, -1, +1, +1, -1, +1, -1, -1, +1, -1, +1, +1, -1, -1, +1, -1, +1, +1, -1, +1, +1, -1, -1, +1, -1, -1, +1, -1, +1, +1, -1, -1, +1, +1, -1, -1, -1, +1, +1, -1, -1, +1, -1, +1, +1, -1, +1, +1, -1, -1, +1, -1, -1, +1, +1, +1, +1, -1, -1, -1, -1, +1, +1, -1, +1, +1, -1, +1, +1, -1, -1, +1, -1, -1, +1, +1, -1, -1, -1, +1, +1, -1, +1, +1, -1, -1, +1, -1, +1, +1, -1, +1, -1, -1, -1, +1, -1, +1, +1, -1, +1, +1, -1, -1, +1, -1, -1, +1, -1, +1, +1, -1, -1, +1, +1, -1, +1, +1, +1, -1, -1, -1, +1, -1, +1, +1, -1, +1, +1, -1, -1, +1, -1, -1, -1, +1, +1, +1, -1, -1, +1, -1, +1, -1, -1, +1, +1, +1, -1, +1, -1, +1, +1, -1, +1, -1, -1, -1, +1, +1, -1, +1, -1, +1, +1, -1, -1, +1, -1, -1, +1, -1, +1, +1, -1, -1, +1, -1, +1, -1, +1, +1, -1, +1, -1, +1, -1, -1, +1, -1, +1, +1, -1, -1, +1, -1, -1, +1, -1, +1, +1, +1, -1, -1, -1, +1, +1, -1, +1, +1, +1, -1, +1, -1, -1, +1, -1, +1, +1, -1, -1, +1, -1, +1, +1, -1, +1, -1, -1, -1, +1, -1, +1, +1, -1, +1, +1, -1, -1, +1, -1, -1, -1, +1, +1, +1, -1, -1, +1, -1, +1, +1, -1, -1, +1, +1, -1, +1, -1, +1, +1, -1, +1, -1, -1, -1, +1, -1, -1, +1, -1, +1, +1, -1, -1, +1, -1, +1, +1, -1, +1, +1, -1, -1, +1, -1, +1, +1, -1, +1, +1, -1, -1, +1, -1, -1, +1, -1, +1, +1, -1, -1, +1, -1, -1, +1, -1, +1, +1, -1, +1, -1, -1, +1, +1, -1, +1, +1, -1, -1, +1, -1, -1, +1, -1, +1, +1, -1, -1, +1, -1, +1, +1, -1, +1, +1, -1, -1, +1, -1, +1, +1, -1, -1, +1, +1, -1, +1, -1, -1, +1, -1, +1, +1, -1, -1, +1, -1, -1, +1, -1, +1, +1, -1, -1, +1, -1, +1, +1, -1, +1, +1, -1, +1, +1, -1, -1, +1, -1, +1, +1, -1, -1, +1, -1, -1, +1, -1, +1, +1, -1, -1, +1, -1, +1, -1, -1, +1, +1, +1, -1, +1, -1, -1, +1, -1, +1, +1, -1, -1, +1, -1, +1, -1, -1, +1, +1, +1, -1, -1, -1, +1, +1, -1, +1, +1, -1, -1, +1, +1, -1, +1, -1, +1, +1, -1, -1, +1, -1, -1, +1, -1, -1, +1, -1, +1, +1, -1, +1, -1, -1, +1, +1, +1, -1, +1, -1, -1, +1, -1, +1, +1, -1, -1, +1, -1, +1, -1, -1, +1, +1, -1, -1, +1, -1, +1, +1, -1, +1, +1, -1, -1, +1, +1, -1, +1, -1, -1, +1, -1, -1, +1, -1, +1, +1, -1, +1, +1, -1, -1, +1, -1, +1, +1, -1, +1, +1, -1, -1, +1, -1, -1, +1, +1, +1, -1, -1, -1, +1, -1, -1, +1, -1, +1, +1, -1, -1, +1, +1, +1, +1, -1, +1, -1, -1, +1, +1, -1, -1, +1, -1, +1, +1, -1, -1, +1, -1, -1, +1, -1, +1, +1, -1, -1, +1, -1, +1, +1, -1, +1, +1, -1, -1, +1, -1, -1, +1, -1, +1, +1, -1, -1, -1, +1, +1, -1, -1, +1, +1, -1, -1, +1, -1, +1, +1, -1, +1, +1, -1, -1, +1, -1, -1, +1, -1, +1, +1, +1, -1, +1, -1, -1, +1, -1, +1, +1, -1, -1, +1, -1, +1, +1, -1, +1, -1, +1, -1, -1, +1, -1, +1, -1, +1, +1, -1, -1, +1, -1, +1, -1, -1, +1, +1, +1, -1, +1, -1, -1, -1, +1, +1, +1, -1, -1, +1, -1, -1, +1, -1, +1, +1, -1, -1, +1, +1, +1, +1, -1, +1, -1, -1, -1, +1, -1, +1, +1, -1, +1, +1, -1, -1, +1, -1, -1, +1, +1, +1, -1, -1, -1, +1, -1, -1, +1, -1, +1, +1, +1, -1, +1, -1, +1, +1, -1, +1, -1, -1, -1, +1, -1, -1, +1, +1, +1, +1, -1, -1, +1, -1, +1, +1, -1, -1, +1, -1, +1, +1, -1, +1, -1, -1, +1, -1, +1, +1, -1, -1, -1, +1, -1, -1, +1, +1, -1, +1, -1, +1, +1, -1, +1, -1, -1, -1, +1, -1, +1, +1, +1, +1, -1, +1, -1, -1, -1, -1, +1, -1, +1, +1, -1, +1, +1, -1, -1, +1, +1, -1, +1, +1, -1, +1, -1, -1, +1, -1, +1, +1, -1, +1, -1, -1}; char basis[SEQSIZE]; int i, j; int fb; // "factor bound" int resetFB; // will act as a boolean char basisVector, initialize; basis[0] = 0; initialize = sequence[1]; basis[1] = initialize; fb = 1; // fb will be the floor of sqrt(n) throughout resetFB = 0; for(i = 2; i < SEQSIZE; i++){ basisVector = initialize * sequence[i]; if(resetFB){ fb++; resetFB = 0; } if((fb+1) * (fb+1) == i){ basisVector *= basis[fb+1]; //we don't want to double count square roots, so handle them separately resetFB = 1; } for(j = 2; j <= fb; j++){ if(i % j == 0) basisVector *= (basis[j] * basis[i/j]); } basis[i] = basisVector; if(basisVector == -1){ printf("%d ", i); } } printf("\n"); return 0; }