Computing a HAP basis
From Polymath Wiki
Jump to navigationJump to search
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;
}