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 }