Computing a HAP basis

From Polymath1Wiki
Revision as of 14:44, 14 January 2010 by Harrison (Talk | contribs)

(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)
Jump to: navigation, 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;
}