P:=proc(q, h) local a, b, d, n, t, x; x:=h+1; a:=1/h; b:=ilog10(h)+1;
d:=h; print(d); t:=1/a; for n from x to q do if
trunc(evalf(a+t/(n-t), 100)*10^(b+ilog10(n)+1))=d*10^(ilog10(n)+1)+n
then b:=b+ilog10(n)+1; d:=d*10^(ilog10(n)+1)+n; a:=a+t/(n-t); t:=n;
x:=n+1; print(n); fi; od; end: P(10^20, 10);
