uptoqdigits(n) = {
	my(res = List([29, 47, 83]));
	forstep(i = 3, n, 2,
		w = withqdigits(i);
		for(j = 1, #w, 
			listput(res, w[j]);
		)	
	);
	Set(res)
}


withqdigits(n) = {
	my(s = sumswithqdigits(n), res = List(), v);
	for(i = 1, #s, 
		v = process(s[i]);
		for(j = 1, #v,
			listput(res, v[j])
		)
	);
	res
}

sumswithqdigits(n) = {
	my(res = List(), ld = [9,11], q = #ld);
	forstep(i = 0, 18, 2,
		forvec(x = vector((n-3)\2, i, [0,18]),
			for(j = 1, q, 
				c = concat([ld[j], Vecrev(x), i, x, ld[j]]);
				fr = fromdigits(c);
				if(issquare(fr, &s),
					if(isprime(s),
						listput(res, c)
					)
				)
			)
		)
	); 
	Set(res)
}

process(v) = {
	my(ld, q = #v, m = v[q\2 + 1]\2, res = List(), qabouthalf = (q-3)\2, part = vector(qabouthalf, i, v[q\2 + 1 + i]));
	if(v[#v] == 11,
		ld = [3,7,9]
	,
		ld = [1,3,7]
	);
	forvec(x = vector(qabouthalf, i, [0, v[q\2+1+i]]),
		for(j = 1, #ld, 
			c = concat([v[#v]-ld[j], Vecrev(part-x), m, x, ld[j]]);
			fr = fromdigits(c);
			if(isprime(fr) && is(fr),
				listput(res, fr)
			)
		)
	);
	Set(res)
	
}

is(n) = {
	my(c = n + fromdigits(Vecrev(digits(n))));
	if(issquare(c, &s),
		return(isprime(s))
	,
		return(0)
	)
}