
Munchausen Numbers and How to Find Them - sndean
https://zach.se/munchausen-numbers-and-how-to-find-them/
======
eridius
The initial example seems broken. In Safari, it looks like

    
    
      3³ + 4⁴ + 3³ + 5⁵ = 343533+44+33+55=3435
    

This is apparently because it shows the same example in two different ways in
the HTML, the first using MathML (which renders the superscripts), and the
second using straight HTML (which doesn't handle the superscripts properly).

In Chrome it's broken a bit differently. Here it render as

    
    
      \displaystyle 3^3 + 4^4 + 3^3 + 5^5 = 343533+44+33+55=3435
    

This is similar to the Safari version except Chrome apparently doesn't support
MathML so it renders the annotation instead (\displaystyle 3^3 + …)

~~~
zacharydenton
Thanks, this post has been converted through several different blog engines
over the years and I hadn't noticed the rendering was broken until now. It's
meant to look more like [https://zach.se/project-euler-
solutions/48/](https://zach.se/project-euler-solutions/48/).

------
cperciva
This unintentionally but very clearly demonstrates how a few minutes of
thought can save a much longer period of computation. Rather than testing all
of the integers between 1 and 500,000,000 (as in the article) or between 1 and
10 * 9^9 (for the trivial upper bound), we can simply test all of the
_multisets of 10 digits_ \-- of which there are a mere 92378. (Or multisets of
_up to_ 10 digits if you take 0^0 = 1, in which case you must consider 184755
possibilities.)

------
wingi
There is only one number - look for perfect numbers.

or:
[https://web.archive.org/web/20091027122658/http://geocities....](https://web.archive.org/web/20091027122658/http://geocities.com/~harveyh/narciss.htm#PPDI)
(Armstrong) Numbers

~~~
ezequiel-garzon
Plenty of even perfect numbers. Find some odd ones!

------
zeveb
Here's a version in Lisp:

    
    
        (defun münchhausen-p (i)
          (eq (loop for number = i then quotient
                    for (quotient remainder) = (multiple-value-list (floor number 10))
                    sum (expt remainder remainder)
                    until (zerop quotient))
              i))
        
        (loop for i from 1 to 500000000 when (münchhausen-p i) do (print i))
    

It takes .002 seconds to check the first 5,000 integers, 2.108 for the first
5,000,000 and 268.867 seconds for the first 500,000,000.

I'm sure there could be some clever tricks to make it run faster, but this is
probably good enough, given that there are, I think, only two such numbers.

~~~
lispm
Don't compare numbers with EQ. Use = or EQL.

~~~
zeveb
Yup, you're completely correct. I'm a bit embarrassed.

------
kazinator

      class munchausen_number_finder {
        real_finder_impl_base &m_f;
      public:
        muchausen_number_finder(real_finder_impl_base &b)
        : m_f(b)
        {
        }
        int get_next() { m_f.get_next(); }
      };

~~~
kazinator
Which is to say, you can find them "by proxy".

/me ducks.

------
wolfram74
someone can spot check me, but since 9^9 = 387420489 < 10^9, there can't be
any Munchausen numbers larger than 10^10 as it there's no way for the sum of
10 numbers less than 10^9 is equal to a number larger than 10^10. I feel like
the statement that only 1 and 3435 are Munchausen numbers should have been
proven by exhaustion by now.

~~~
gus_massa
Almost. The numbers can be repeated, if N = 999999999 then Mun(N) = 9 * 9^9 =
3486784401 ~=3.5E9. But the general idea is correct.

[spoiler alert] If we assume that the number N has D digits, then

    
    
      10^(D-1) <= N < 10^D 
    

Now, N has D digits, and each one can be at most 9, so

    
    
      Mun(N) <= 9^9 * D = 387420489 * D < 10^9 * D
    

If we assume that there exist a number N that is a Munchausen number with D
digits, then N = Mun(N) so we can combine the left part of the first equation
with the right part of the second equation

    
    
      10^(D-1) <= N = Mun(N) < 10^9 * D
    

If we only keep the extreme parts of the inequality

    
    
      10^(D-1) < 10^9 * D
    

And now we can take log (in base 10)

    
    
      D-1 < 9 + log(D)
    

And put all the D in the left side:

    
    
      D – log(D) < 10
    

Taking the derivatives, the function f(x)=x-log(x) is essentially increasing
if x is big enough, probably when x>1 or x>2, I'm a bit lazy. (Remember to use
the ln(10) in the derivative because we are using log=log_10 instead of
ln=log_e.)

So we have only to find a value of D that is big enough. For example, if D=12,
then

    
    
      12-log(12) > 12-log(100) = 12-2 = 10
    

So if D is 12 or more, we have a contradiction, so there are no Munchausen
numbers with 12 digits or more.

A more careful analysis can reduce this value a little, but a for between 1
and 10^12-1 to test all the cases should finish in a modern computer in a
reasonable time, specially replacing the pow(n, n) calculation with a cached
value.

~~~
wolfram74
Exploring this in other bases would be amusing, maybe not particularly
enlightening though.

~~~
Someone
[https://oeis.org/A166623](https://oeis.org/A166623)

------
reikonomusha
> Here’s the program written in Python. I used it to find every Munchausen
> number less than 500,000,000. After thirty minutes or so of intense
> computation, it turns out that the only Munchausen numbers are 1 and 3435.

This is an awfully long time to wait for something that is very clearly
optimizable. If we are searching for Munchausen numbers, then we are probably
going to do so in a brute force way. (We could optimize this with various
symmetries, such as digit permutations, but let's ignore this for now.)

If we are going to brute force it, then we are probably not going to exceed
64-bits, which means we can avoid using bignums altogether. Written correctly,
this would be extremely fast in C.

With this code [0], on my machine, with only meager optimization effort (not
thinking too hard about things other than data types), I can go through 3
billion in less than a minute:

    
    
        > (time (find-munchausen 3000000000))
        Computing all Munchausen numbers below three billion.
        Evaluation took:
          57.087 seconds of real time
          57.059879 seconds of total run time (56.840355 user, 0.219524 system)
          99.95% CPU
          159,836,797,604 processor cycles
          401,488 bytes consed
          
        (0 1 3435 438579088)
    

The assembly code for the Munchausen test is relatively tight (though not at
all as good as handwritten, of course), which makes evident why we check so
quickly.

    
    
        ; 8E2:       4C8B05A7FFFFFF   MOV R8, [RIP-89]                ; #(0 1 4 27 ...)
        ; 8E9:       488BD9           MOV RBX, RCX
        ; 8EC:       0F1F4000         NOP
        ; 8F0: L0:   48BA9A99999999999919 MOV RDX, 1844674407370955162
        ; 8FA:       488BC3           MOV RAX, RBX
        ; 8FD:       48F7E2           MUL RAX, RDX
        ; 900:       488BFA           MOV RDI, RDX
        ; 903:       4883E7FE         AND RDI, -2
        ; 907:       B814000000       MOV EAX, 20
        ; 90C:       488BF7           MOV RSI, RDI
        ; 90F:       48D1FE           SAR RSI, 1
        ; 912:       480FAFF0         IMUL RSI, RAX
        ; 916:       4829F3           SUB RBX, RSI
        ; 919:       488BF3           MOV RSI, RBX
        ; 91C:       488BDF           MOV RBX, RDI
        ; 91F:       498B44B001       MOV RAX, [R8+RSI*4+1]
        ; 924:       4839C8           CMP RAX, RCX
        ; 927:       7E0E             JLE L2
        ; 929:       B917001020       MOV ECX, 537919511
        ; 92E: L1:   488BD1           MOV RDX, RCX
        ; 931:       488BE5           MOV RSP, RBP
        ; 934:       F8               CLC
        ; 935:       5D               POP RBP
        ; 936:       C3               RET
        ; 937: L2:   4829C1           SUB RCX, RAX
        ; 93A:       4885DB           TEST RBX, RBX
        ; 93D:       75B1             JNE L0
        ; 93F:       4885C9           TEST RCX, RCX
        ; 942:       B917001020       MOV ECX, 537919511
        ; 947:       41BB4F001020     MOV R11D, 537919567             ; T
        ; 94D:       490F44CB         CMOVEQ RCX, R11
        ; 951:       EBDB             JMP L1
    

[0] [https://bitbucket.org/tarballs_are_good/lisp-
random/raw/mast...](https://bitbucket.org/tarballs_are_good/lisp-
random/raw/master/munchausen.lisp)

~~~
lokedhs
I find it quite remarakable that Python is so popular with maths people, even
though it's so horribly slow. I've seen similar differences between SBCL and
Python quite often, when I have tried reimplmenting code like that in Common
Lisp.

~~~
hcarvalhoalves
> I find it quite remarakable that Python is so popular with maths people,
> even though it's so horribly slow.

It's popular because of Numpy/Scipy, then the speed of native loops doesn't
matter.

