S&SMC #3: Final Remarks [LONG!]
Message #15 Posted by Ex-PPC member on 14 June 2002, 11:17 p.m.,
in response to message #1 by Ex-PPC member
Thanks to all of you who were interested in my S&SMC #3, the challenge
is over and here are some final notes and remarks:
As has been mentioned in a post, the numbers who are equal to some
mathematical function of their digits are generally referred to as
"Narcissistic Numbers". In Challenge #3, we were looking for numbers
of N digits that are equal to the sum of the N-th powers of their digits.
In the mathematical literature, these are known as "PluPerfect Digital
Invariants", i.e: PPDI. They can be defined for bases other than 10, but
in what follows, base 10 will always be assumed, and we'll call Nth-order PPDI
to a PPDI having exactly N digits. A lot is known about them, for example:
- There is only a finite number of PPDIs, exactly 88. It is very easy
to demonstrate that their number is finite, but finding that there are
exactly 88 of them does require a lot of ingenuity and computer muscle.
- There are no PPDIs of order N when N=2,12,13,15,18,22,26,28,30,36
nor for N>39, so 39 is the largest possible order.
- There is one and only one PPDI when N=6,10,14,20,32,34,37,38
- The largest, 88th PPDI of order N=39 is the 39-digit number:
115132219018763992565095597973971522401
which is equal to the sum of its digits raised to
the 39th power.
- The largest prime PPDi is the 23-digit numer:
35452590104031691935943
which is equal to the sum of its digits raised to the 23rd power.
The solutions for S&SM Challenge #3 are all PPDI of orders N from
1 to 10, i.e:
Order N Solutions
-----------------------------------------------------
1 1, 2, 3, 4, 5, 6, 7, 8, 9
2 none
3 153, 370, 371, 407
4 1634, 8208, 9474
5 54748, 92727, 93084
6 548834
7 1741725, 4210818, 9800817, 9926315
8 24678050, 24678051, 88593477
9 146511208, 472335975, 534494836, 912985153
10 4679307774
Can our beloved HP calculators find these solutions ? Yes, of course.
As I said, almost any model from the HP-25 onwards can find all solutions
up to N=4 with relative ease, even using a straightforward brute-force
approach. Finding solutions for orders N=5 to N=10 in a reasonable amount
of time does require a combination of faster models, such as the HP-71B,
HP-32S/SII, HP-42S or HP-48/49, as well as much more refined, optimized
programming, and even a completely different approach.
Let's see an example of how can you proceed in order to achieve the goal.
In this example, we'll compute all 3 solutions for order N=4, using an HP-71B,
as its BASIC language helps to make the programming easier to understand.
We'll start with a vey simple brute-force approach, then we will refine it
stepwise, creating better and faster versions which will allow us to reach
higher orders. Then, once we reach the limits of our direct approach, we
will hint at a newer, much faster approach which will deliver what we want.
Then, a direct conversion of the HP-71B BASIC code to the HP-41C's native
RPN will be shown.
Caveat emptor: You should bear in mind that, in all cases, the programs shown have been produced strictly for
demonstration purposes and are not meant to be the ultimate
programming solution for this challenge, or even highly
optimized, so I don't claim or intend them to be
state-of-the-art programming at all.
Program P1-71B:
As we are looking for all 4-digit numbers equal to the sum of the 4th
powers of its digits, this means we are looking for numbers which are solutions
of this Diophantine equation:
[ABCD] = 1000*A+100*B+10*C+D = A^4+B^4+C^4+D^4
where A goes from 1 to 9, and B,C,D all go from 0 to 9. A very simple
program, made essentially of four nested FOR-NEXT loops is
quickly produced:
10 DESTROY ALL @ DIM A,B,C,D @ STD @ DELAY 0,0
20 FOR A=1 TO 9 @ FOR B=0 TO 9 @ FOR C=0 TO 9 @ FOR D=0 TO 9
30 IF A^4+B^4+C^4+D^4=1000*A+100*B+10*C+D THEN DISP A;B;C;D
40 NEXT D @ NEXT C @ NEXT B @ NEXT A @ DISP "DONE"
When run, this program displays all three solutions [1634, 8208, 9474], then
terminates displaying "DONE" after verifying there are no more.
Running time: T = 2050 sec.
Program P2-71B:
One of the easiest ways to reduce the running time of a program is to
simplify the computations within the innermost loop. In this case, program
P1 is continually evaluating the expression:
1000*A+100*B+10*C+D
inside the innermost loop. But a bit of thinking or a quick hand simulation
will convince you that this is just the value of the number itself, which
is being incremented by one in each iterarion. So we can simply initialize
this value to 1000 (A=1, B=0, C=0, D=0), and then increment it after each
iteration, thus saving a lot of computation:
10 DESTROY ALL @ DIM A,B,C,D,N @ STD @ DELAY 0,0
20 N=1000 @ FOR A=1 TO 9 @ FOR B=0 TO 9 @ FOR C=0 TO 9 @ FOR D=0 TO 9
30 IF A^4+B^4+C^4+D^4=N THEN DISP N
40 N=N+1 @ NEXT D @ NEXT C @ NEXT B @ NEXT A @ DISP "DONE"
Running time: T = 1854 sec. (1.11 times faster than P1)
Program P3-71B:
A further refinement comes when we realize that we are continually computing
the 4th powers of the digits of N inside the innermost loop. But there are
only 10 distinct values, namely 0^4, 1^4, ..., 9^4, so we can save a lot
of time pre-calculating them outside of all loops and storing them on a
suitably dimensioned array, then simply recalling the values of the powers when
needed, instead of computing them anew each time. As it is much faster to recall an array
element than to raise a number to a power, great time savings are to be expected:
10 DESTROY ALL @ OPTION BASE 0 @ DIM A,B,C,D,N,P(9) @ STD @ DELAY 0,0
15 FOR A=0 TO 9 @ P(A)=A^4 @ NEXT A
20 N=1000 @ FOR A=1 TO 9 @ FOR B=0 TO 9 @ FOR C=0 TO 9 @ FOR D=0 TO 9
30 IF P(A)+P(B)+P(C)+P(D)=N THEN DISP N
40 N=N+1 @ NEXT D @ NEXT C @ NEXT B @ NEXT A @ DISP "DONE"
Running time: T = 460 sec. (4.46 times faster than P1)
Program P4-71B:
Is that all ? Can't we still do more optimization ? Yes. Look at the IF
at line 30. We are constantly recalling and adding up the powers of the digits
to compare them against the value of N. But it's clear that inside the
FOR D=... loop, the value of P(A)+P(B)+P(C) is invariant, as it does not
depend on the loop variable, D. Same applies to P(A)+P(B) inside the FOR C
loop, and P(A) inside the FOR B loop. Keeping this is mind, we avoid
recomputing sums and invariant values repeatedly, but we'll simply compute
them once outside of the respective loops, and store them in intermediate
variables, like this:
10 DESTROY ALL @ OPTION BASE 0 @ DIM A,B,C,D,N,S,T,U,P(9) @ STD @ DELAY 0,0
15 FOR A=0 TO 9 @ P(A)=A^4 @ NEXT A @ N=1000
20 FOR A=1 TO 9 @ S=P(A) @ FOR B=0 TO 9 @ T=S+P(B)
30 FOR C=0 TO 9 @ U=T+P(C) @ FOR D=0 TO 9 @ IF U+P(D)=N THEN DISP N
40 N=N+1 @ NEXT D @ NEXT C @ NEXT B @ NEXT A @ DISP "DONE"
Running time: T = 320 sec. (6.41 times faster than P1)
Program P5-71B:
There's still one trick out of our sleeve. Just notice what happens when, for
instance, we reach the inner FOR C and FOR D loops when A=9 and B=8. At this
point, the sum is already 9^4+8^4 = 10657, which exceeds the maximum possible
value for N, namely N=9999. So, further adding the 4th powers of C and D is no use, since the resulting sum can never equal N. This means we can skip
altogether the full FOR C and FOR D inner loops in this case, so
saving 100 full iterations. Applying this simple idea gives us the program:
10 DESTROY ALL @ OPTION BASE 0 @ DIM A,B,C,D,N,S,T,U,V,P(9) @ STD @ DELAY 0,0
15 FOR A=0 TO 9 @ P(A)=A^4 @ NEXT A @ N=1000
20 FOR A=1 TO 9 @ S=P(A)
25 FOR B=0 TO 9 @ T=S+P(B) @ IF T>N THEN N=N+(10-B)*100 @ GOTO 80
30 FOR C=O TO 9 @ U=T+P(C) @ IF U>N THEN N=N+(10-C)*10 @ GOTO 70
40 FOR D=0 TO 9 @ V=U+P(D) @ IF V>N THEN N=N+10-D @ GOTO 60 ELSE IF V=N THEN PRINT N
50 N=N+1 @ NEXT D
60 NEXT C
70 NEXT B
80 NEXT A @ DISP "DONE"
Running time: T = 242 sec. (8.47 times faster than P1)
[By the way, just in case you are wondering, declaring the variables with INTEGER precision
instead of floating-point DIM is actually slower! HP-71B BASIC always works
with floating-point variables internally, so declaring them INTEGER just adds
a lot of conversions to floating-point precision and back to integer]
So you see, applying just a few, simple common sense ideas to our initial try
has resulted in a program that, while still being very simple and not much larger,
nevertheless runs nearly an order of magnitude faster.
Which is better, the savings are even greater for larger orders of N, so that
this approach allows us to compute all solutions (and to certify there are no
more) in these running times:
Order N P1 P2 P3 P4 P5
-------------------------------------------------
3 02:32 02:21 00:40 00:32 00:25
4 34:10 30:54 07:40 05:20 04:02
5 - - - - 37:20
6 - - - - 5:07:13
Remember, these times are measured from start to until the program stops by
itself, not merely until the last solution is displayed. Times for an HP-32S/SII,
HP-42S or HP-48/49 would be even better, as their CPUs are much faster than the old 640 Khz Saturn CPU of the 71B.
However, for N=7 to N=10 our current approach does not work in reasonable times, and
we find we have two possible options:
- a) stick to our current approach, but change to a faster language, using for
instance FORTH instead or BASIC, or better still, using Saturn
Assembler language (using the 71B Forth/Assembler ROM, for instance).
Converting our program P5 to Saturn assembler is not difficult at all, and
we can take advantage of fast custom-made all-integer routines, much faster
than the floating-point ones we're using in BASIC. The resulting binary
subprogram (not BASIC keyword) is faster than our BASIC version by two full orders of magnitude,
and so allows us to go up to N=7 in under 30 min., N=8 in some 4 hours,
and N=9 in roughly a day and a half continuous running time. The final summit,
N=10, would take more than 2 weeks, so this approach's usefulness ends at N=9.
- b) change our current approach for a completely different one. To that effect,
we notice that (for our example, N=4) we are exhaustively searching all
arrangements from 1000 to 9999, using 4 nested loops, so a total of 9000
cases are tested in all. Our trick in P5 avoids some of these, but only to the point
were we have an exponential increase of roughly 8.5x per additional digit instead of 10x. Useful but insufficient.
However, we are testing far more cases than we actually need !. For instance,
for N=6754, we are computing the sum 6^4+7^4+5^4+4^4 and comparing it
versus 6754. But the sum itself is invariant whatever the order of the
digits may be: we have Sum4(6754) = Sum4(6745) = Sum4(4567) = ...
So, the magic word is lists: you just need to test all different lists,
order doesn't matter !! This tremendously reduces the number of cases we need
to test, to the point where, with clever list-generating programming, we can
succeed even for N=10 ! And, with the help of a fast PC, we can even find
all 88 PPDIs, up to order 39, in a reasonable amount of time.
Of course, I won't give here the solution using this lists approach, that's
for you to take over. But for the sake of all of you who tried this challenge
on your HP-41Cs or your HP-42S, I'm including a straight conversion of P4 to
HP-41C/CV/CX's RPN, for the case N=3:
Program P4-41C:
01 LBL "N3" 16 DSE 16 31 RCL 14
02 9 17 LBL 01 32 +
03 STO 00 18 RCL IND 10 33 RCL 15
04 LBL 00 19 STO 13 34 X=Y?
05 RCL 00 20 RCL 16 35 VIEW X
06 3 21 STO 11 36 ISG 15
07 Y^X 22 LBL 02 37 LBL 04
08 STO IND 00 23 RCL IND 11 38 ISG 12
09 DSE 00 24 RCL 13 39 GTO 03
10 GTO 00 25 + 40 ISG 11
11 100 26 STO 14 41 GTO 02
12 STO 15 27 RCL 16 42 ISG 10
13 1.009 28 STO 12 43 GTO 01
14 STO 10 29 LBL 03 44 "DONE"
15 STO 16 30 RCL IND 12 45 PROMPT
This 45-line program will compute all four solutions for N=3, in these times:
Time to find and display 1st solution (153) .... 24 sec.
2nd solution (370) .... 1 min 37 sec.
3rd solution (371) .... 1 min 37 sec.
4th solution (407) .... 1 min 49 sec.
"DONE" ................ 5 min 6 sec.
This program will also run unchanged on an HP-42S, only
much faster.
That's all. Next S&SMC #4 soon, stay tuned !
|
Go back to the main exhibit hall