addvec(vlist, v) = {
	for (i=1, #v,
		listput(vlist, v[i]);
	);
	vlist;
}

fill(vc) = {
	my(vlist = List());
	for (i=1, #vc,
		vlist = addvec(vlist, vc[i]);
	);
	vlist;
}

kin(vlist, k) = vecsearch(vecsort(Vec(vlist), , 8), k);
precp(p) = precprime(p-1);
supf(v) = if (#v > 1, vector(#v-1, k, v[k+1]), []);
supl(v) = if (#v > 1, vector(#v-1, k, v[k]), []);

G(x, p) = {
	if (p==2, return ([1, 2]));
	if (p==3, if (x < 18, return ([1, 2]), return ([1, 3, 6, 2])));
	my(px = precprime(sqrtint(x\2)));
	if (px == 0, return);
	if (p > px, return (G(x, px)));
	my(retg = List());
	my(p1 = precp(p));
	my(p2 = precp(p1));
	retg = addvec(retg, [1]);
	retg = addvec(retg, p*G(x/p, p2));
	retg = addvec(retg, [2*p*p1]);
	retg = addvec(retg, p*p1*supl(G(x/(p*p1), p2))); 
	retg = addvec(retg, supf(G(x, p1)));
	Vec(retg);
}

g(p) = G(2*p^2, p);

d(p) = {
	my(v = g(p));
	nv = concat(vector(#v-1, k, v[k+1]), 1);
	nv * nextprime(p+1);
}

lista(nn) = {
	vd = vector(nn+1, k, if (k==1, [1], d(prime(k-1))));
	vc = vector(nn+1, k, vd[k]);
	vq = vector(nn+1);
	vq[1] = 1;
	for (k=1, nn,
		vlist = fill(vc);
		if (kin(vlist, k),
			vq[k+1] = vq[k];
		,
			vq[k+1] = max(nextprime(vq[k]+1), nextprime(k^2));
			p = vq[k+1];
			if (primepi(p) > nn, break);
			vc[primepi(p)+1] = concat([p*k^2, k, nextprime(p+1)*p*k^2], d(p));
		);
	);
	Vec(vlist);
}