a(n) = {
	res = 0;
	qchecks = 0; 
	pows2 = vector(n, i, 1<<i);
	my(i);
	cans = [1..n];
	bonus = 0; 
	noncans = List(vector(sqrtint(n)-2, i, (i+2)^2));
	forprime(p = ceil(n/2), n, 
		listput(noncans, p);
		bonus++;
	);
	triples = vector(n, i, maketriplesfora(i));
	res-=bonus;
	cans = setminus(Set(cans), Set(noncans)); 
	ulim = n; 
	for(i = 1, #cans, 
		processfora(1<<cans[i], i, 0);
	);
	\\print(qchecks); 
	res + bonus
	
}

processfora(t, n, s) = {
	my(ct, ci, cs = s+1, i, csn = cans[n], hct); 
	qchecks++; 
	ct = t + pows2[csn];
	hct = hammingweight(ct); \\ I'd like to pass this on as variable s to ease calculation but somehow it messes up.
	if(hct + #cans < res + n, 
		return
	);
	for(i = 1, #triples[cans[n]], 
		if(bitand(triples[csn][i], ct) == triples[csn][i],
			return
		);
	);
	res = max(hct, res); 
	for(i = n+1, #cans, 
		if(hct + #cans < res + i,
			return;
		);
		ci = i; 
		processfora(ct, ci, cs);
	)
}

maketriplesfora(n) = { \\ could be faster but good enough for 'small' numbers. 
	my(res = List());
	for(i = 1, #cans,
		for(j = i+1, #cans,
			if(issquare(cans[i]*cans[j]*n) && cans[j] < n,
				listput(res, (1<<cans[i]) + (1<<cans[j]) + (1<<n));
			);
		)
	); res
	
}