Hacker News new | past | comments | ask | show | jobs | submit login
Munchausen Numbers and How to Find Them (zach.se)
36 points by sndean on Sept 21, 2016 | hide | past | favorite | 20 comments

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 + …)

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/.

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.)

There is only one number - look for perfect numbers.

or: https://web.archive.org/web/20091027122658/http://geocities.... (Armstrong) Numbers

Plenty of even perfect numbers. Find some odd ones!

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))
    (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.

Well, 9^9 = 387,420,489. So a ten-digit Munchausen number that had enough nines in it could be as big as 3,874,204,890. So going to 500,000,000 isn't quite enough. You can stop at 4 billion, though...

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

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

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

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

/me ducks.

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.

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.

For an 11-digit number, Mun(N) <= 11 * 9^9 = 4,261,625,379 - which only has 10 digits. Therefore, no Munchausen number can have 11 (or more) digits.

For a 10-digit number, Mun(N) <= 10 * 9^9 = 3,874,204,890. So you can stop there.

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

https://oeis.org/A046253 agrees, but adds 0 and 438579088, using 0^0 = 0.

Also, http://mathworld.wolfram.com/MuenchhausenNumber.html gives an explanation why they are called Münchhausen Numbers. It states:

these numbers "raise themselves" analogously to way in which Baron Hieronymus von Münchhausen allegedly raised himself by riding a cannonball.

I would think it refers to bootstrapping (the baron pulled himself out of a morass by pulling on his bootstraps)

And for those playing with this in their favorite language: I would precompute an array of ten d^d values and hoist the d^d computation out of the loop for a very nice speed up.

Finally: int may not be sufficient in C if you iterate up to 9 digits. 9^9 needs 29 bits => 999,999,999 would need over 32 => unsigned int already is living very dangerously (probably too dangerously)

> 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...

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.

> 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.

Guidelines | FAQ | Lists | API | Security | Legal | Apply to YC | Contact