Prolog实现Münchausen数求解：如何提速至搜索[0,10^10-1]区间？

下面的代码可求解部分Münchausen数（完美数位不变数），遵循0^0 = 0的约定：

digits(0, P, Q, R) :- !, P = Q, R = P.
digits(N, P, Q, R) :- 
   M is N-1,
   between(0, 9, D),
   S is P*10+D,
   (D = 0 -> T = Q; T is Q+D^D),
   digits(M, S, T, R).

在SWI-Prolog中执行该代码，查询digits(4, 0, 0, R)可得到如下结果：

?- digits(4, 0, 0, R).
R = 0 ;
R = 1 ;
R = 3435 ;
false.

上述代码仅覆盖搜索区间[0,10^4-1]，如何在短时间内完成[0,10^10-1]区间的搜索？

2023年3月13日更新：
根据一篇论文，我们可以彻底解决是否存在其他Münchausen数的问题：

关于3435的一个奇特性质
Daan van Berkel - 2009

这十分令人惊叹，因为Münchausen数的数量是有限的。

内容的提问来源于stack exchange，提问作者user17524790