Exhausting decimal Munichausen numbers in Prolog under convention 0^0 = 1

Discussion
1

A decimal number 3435 has a special property that

3435 = 3^3+4^4+3^3+5^5,

which is called a Munichausen Number: Munichausen Number (Wikipedia) As 1=1^1 1 is a Munichausen.

Is there any other Munichausen number ? It is not difficult to see that a decimal Munichausen must be less than 10^{10}, provided that convnention 0^0=1 is exploited. Note that
SWI-Prolog exploits it.

To see that 1 and 3435 are all Munichausen, I tried to check it by naive generate-test way for all positive integers < 10^{10.} But 10^{10} seemed to be too big for
that naive method, though I did not try any constraint based technology yet.

Anyway finally building a (irreducible) BDD tree for all 10^{10} decimal digits lists and check
Munichausen for each list of digits. Fortunately it terminates showing the expected result that a decimal Munichausen are only 1 and 3435, after running the query below on m2 mac mini for more than one hour,

However, I am afraid that using BDD is too much brute-force for the problem. I would like to hear about any clever prolog method to exhaust Munichausen numbers less than specified limit (e.g. 10^{10}) . Thanks in advance.

% ?- time( (zdd { K=10, numlist(0, 9, Ns) }, bdd_list_raise(Ns, K, X),
% munichausen(0, 0, X), {fail})).

%@ munichausen_number(3435)
%@ munichausen_number(1)
%@ % 100,000,132,272 inferences, 4716.645 CPU in 4737.086 seconds (100% CPU, 21201541 Lips)
%@ false.

munichausen(X, X, 1, _):- X>0, writeln(munichausen_number(X)).
munichausen(X, Y, U, S):- U>1,
	cofact(U, t(I, L, R), S),
	munichausen(X, Y, I, L, R, S).
%
munichausen(X, Y, _, L, _, S):-	munichausen(X, Y, L, S).
munichausen(0, _, 0, _, R, S):-!, munichausen(0, 0, R, S).
munichausen(X, Y, I, _, R, S):-
	X0 is 10*X + I,
	Y0 is Y + I^I,
	munichausen(X0, Y0, R, S).
3

Initial version of codes of mine did not work for the case N = 9, where N is the length of decimal number, i.e. N=3 for 123.

It is strange for me why BDD version works even for N=10, which is enough for the exhaustive search for the Munichausen numbers. As the BDD version is also based on generate and test method, I doubt that my solution is good.

Does you codes terminate well for N=10 ?

5

Sound great. I will check your codes.

My codes is flexible between the two conventions, and for choice of base (radix). I will play around on these variants.

I have to check this. Thanks.

6

No existence of decimal Munichausen number > 10^{10} came to me from the following simple experimental arithmetic query． Am I missing ?

% ?- between(1, 100, N), 10^N > N*(9^9), writeln(N), fail.
9

Thanks for links. I will go to the links soon.

I ran a script go/3 to enumerate all Munichausen numbers upto radix 10 for both conventions 0^0=0 and 0^0=1. It seems to return to correct answers. It took time about 4 hours, which should be revised.

go(R, C, Nums) enumerates Nums of all Munichausen numbers for radix R and convention flag C.

% ?- between(2, 10, I), member(J, [0,1]),
%	go(I, J, Nums), writeln(go(I,J) = Nums),
%	fail; true.
%@ go(2,0)=[1,0]
%@ go(2,1)=[2,1]
%@ go(3,0)=[8,5,1,0]
%@ go(3,1)=[8,5,1]
%@ go(4,0)=[55,29,28,1,0]
%@ go(4,1)=[55,29,1]
%@ go(5,0)=[264,28,1,0]
%@ go(5,1)=[1]
%@ go(6,0)=[3416,3164,1,0]
%@ go(6,1)=[3416,3164,1]
%@ go(7,0)=[3665,1,0]
%@ go(7,1)=[3665,1]
%@ go(8,0)=[257,256,1,0]
%@ go(8,1)=[1]
%@ go(9,0)=[16871323,923362,917139,96446,28,27,1,0]
%@ go(9,1)=[923362,96446,28,1]
%@ go(10,0)=[438579088,3435,1,0]
%@ go(10,1)=[3435,1]
12

I get a faster method with Prolog first argument indexing:

?- time((canonball, munchhausen(_), fail; true)).
% 4,633,344 inferences, 0.594 CPU in 0.589 seconds (101% CPU, 7803527 Lips)
true.

?- munchhausen(R).
R = 0 ;
R = 1 ;
R = 3435 ;
R = 438579088 ;
false.

The code is based on the observation that the fixpoint
equation for a split number m*10^k + n can be reordered
so that m is on one side and n is on the other side:

m*10^k - F(m) = F(n) - n

The rest is tabulating F(n) - n:

canonball :-
   retractall(cache(_,_)),
   between(0, 99999, N), map(N, Y), C is Y-N,
   assertz(cache(C, N)), fail; true.

munchhausen(R) :-
   between(0, 99999, M), map(M, X), B is 100000*M-X,
   cache(B, N), R is 100000*M+N.

map(0, X) :- !, X = 0.
map(N, X) :-
   M is N//10,
   map(M, Y),
   D is N mod 10,
   (D = 0 -> X=Y; X is Y+D^D).
13

Congratulation !

It seems worth more than to say great.

14

But computing the table up to R=16 is still an open problem for me.
Need more speed-up ideas. Up to R=10 I can verify your table,
It now needs ca. 5 seconds with Dogelog Player on nodeJS:

?- time((between(2,10,R), write(R), write(': '),
(canonball(R), munchhausen(R,S), write(S), write(' '), fail; true),
nl, fail; true)).
2: 0 1
3: 0 1 5 8
4: 0 1 28 29 55
5: 0 1 28 264
6: 0 1 3164 3416
7: 0 1 3665
8: 0 1 256 257
9: 0 1 27 28 96446 917139 923362 16871323
10: 0 1 3435 438579088
% Time 5170 ms, gc 15 ms, 2877753 lips
true.

But if you have a larger radix R, the square root of the range that
needs to be first tabled and later scanned increases more than only
exponential. Its not that we have a constant C and the size is C^R,

its much worse, namely R^(R//2). Quite a challenge to go higher.

15

I have revised the estimation of upper-bound of width of possible Munichausen numbers to the radix - 1. I am fairy sure of correctness of this revision. In fact, I have proved it, which uses interval picturs [r^n, r^(n+1)) (left closed, right open) as you used in a post of of yours in this thread. I’m not sure this revision is relevant to your open problem for case R=16, but it would be good if it were the dase.

16

When R, the radix, is a parameter, I am using this code:

canonball(R) :-
   retractall(cache(_,_)),
   U is R^(R//2)-1,
   between(0, U, N), map(N, R, Y), C is Y-N,
   assertz(cache(C, N)), fail; true.

munchhausen(R, S) :-
   U is R^(R//2),
   V is R^(R-R//2)-1,
   between(0, V, M), map(M, R, X), B is U*M-X,
   cache(B, N), S is U*M+N.

map(0, _, X) :- !, X = 0.
map(N, R, X) :-
   M is N//R,
   map(M, R, Y),
   D is N mod R,
   (D = 0 -> X=Y; X is Y+D^D).

I didn’t try R=16 yet. I guess it takes quite long or it will
crash the computer. The cannonball iteration will use a
U=4294967295, that gives a quite big cache/2 table.

I think its currently eating more and more memory, and
will soon crash. I see:

?- canonball(16).

2b4a1cde63cfa54e76d0978d5bc08c2c76683bda

Yeah after some while I get:

e750366b18aa033abae61d91a11a0516661980f2