first(n) = {my(res = vector(n)); res[1] = 1; m = Map(); 
	mapput(m, 1, 1); 
	for(i = 4, 6, mapput(m, i, i-2));
	mapput(m, 9, 5); 
	lastWithFirstdigit = vector(5);
	lastWithFirstdigit[1] = 1;
	for(i = 2, n, 
		lt = (res[i-1]/10^valuation(res[i-1], 10))%10; 
		ind = mapget(m, lt);
		nt = nextsquare(lastWithFirstdigit[ind], lt); 
		res[i] = nt;
		lastWithFirstdigit[ind] = nt
	); res
} 

nextsquare(sq, lt) = {
	my(s = sqrtint(sq) + 1); c = s^2;
	ld = leadingdigit(s^2); 
	if(ld == lt,
		if(c % 10 == 0, 
			return(nextsquare(c, lt))
		,
			return(c)
		)
	,
		if(ld > lt, 
			c*=10;
		);
		cs = ceil(sqrt(lt*10^logint(c, 10)))^2; 
		while(leadingdigit(cs) != lt,
			c*=10; 
			cs = ceil(sqrt(lt*10^logint(c, 10)))^2;
		);
		
		if(cs % 10 == 0, 
			return(nextsquare(cs, lt))
		,
			return(cs)
		)
	)
}

leadingdigit(n) = n\10^logint(n, 10);