first(n) = { my(prods = vector(n^2), sums = vector(2*n^2), res = vector(n), news, t, leastunseensum = 3); res[1] = 3; prods[1] = 1; sums[1] = sums[2] = 1; for(i = 2, n, news = List(); if(isprime(i), news = i*[1..i] , if(i % 2 == 0, news = i*[(i/2)..i]; , for(j = 1, i, d = divisors(i*j); if(d[#d \ 2 + 1] >= i, listput(news, i*j); ) ) ); ); for(j = 1, #news, prods[news[j]] = 1; if(sums[news[j]] == 0, t++; sums[news[j]] = 1; ) ); i2 = i^2; for(k = 1, #news, for(j = max(leastunseensum, news[k]+1), i2 + news[k], if(sums[j] == 0 && prods[j - news[k]] == 1, sums[j] = 1; t++; ) ); ); for(j = leastunseensum, #sums, if(sums[j] == 0, leastunseensum = j; break ) ); res[i] = res[i-1] + t; t = 0; ); kill(prods); kill(sums); return(res); }