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);
}