The Museum of HP Calculators

HP Forum Archive 21

[ Return to Index | Top of Index ]

Programming exercise
Message #1 Posted by Gerson W. Barbosa on 22 Aug 2012, 10:00 p.m.

What is my licence plate number?

Hint #1: The four-digit number is the smallest number expressible as the sum of two cubes in two different ways, that is, n = a^3 + b^3 = c^3 + d^3; a, b, c and d are positive integers;

Hint #2: n < 2000;

Hint #3: the three-letter pefix is CAB backwards :-)

My straightforward approach (not to call it 'brute-force' :-) could be used on the HP-41CV, but I have preferred to do it on my HP 50g. My RPL program is somewhat clumsy because I wrote a QBASIC program first and then converted it to RPL, but it does find n, a, b, c and d in 70 seconds on the real calculator when the input is 1999. When the input is 9999, another larger solution is found, but it takes about 10 minutes before the solutions are shown. The ten solutions less than 100,000 are found in about 9 minutes on the emulator. There is much room for improvement as no mathematical approaches have been tried-- neither could I. (I would have chosen API-3142, but it wasn't available anymore. On second thought I should have chosen something as simple as ABC-1234 :-)

Have fun!

      
Re: Programming exercise
Message #2 Posted by Don Shepherd on 22 Aug 2012, 10:45 p.m.,
in response to message #1 by Gerson W. Barbosa

So how long have you been driving a taxicab in a suburb of London?

            
Re: Programming exercise
Message #3 Posted by Gerson W. Barbosa on 22 Aug 2012, 11:15 p.m.,
in response to message #2 by Don Shepherd

Ah, you know the story! Well, most everyone here does, I think :-)

                  
Re: Programming exercise
Message #4 Posted by Valentin Albillo on 23 Aug 2012, 4:27 a.m.,
in response to message #3 by Gerson W. Barbosa

Quote:
Ah, you know the story! Well, most everyone here does, I think :-)

Ramanujan, perchance ?

Best regards from V.
                        
Re: Programming exercise
Message #5 Posted by LHH on 23 Aug 2012, 4:46 a.m.,
in response to message #4 by Valentin Albillo

This is just the kind of problem I like to try to solve in my head when I go to bed. Usually it's calculating the time to a star at various speeds and distances (or something along those lines). I may have mentioned before my (young) kids and I once tried to estimate the number of atoms in the universe and came up with 10^75. Recently I saw an estimate of 10^80 elemental particles so we actually came pretty close! This would have been more difficult just remembering all the possible cubes but I'm sure it would have put me to sleep quickly (which is the desired effect)! Anyway, just doodling I wrote down all the possible cubes and have an answer (BAC-1729) but I did cheat by using a biological computer (I will accept my disqualification in a manly fashion).

                              
Re: Programming exercise
Message #6 Posted by Valentin Albillo on 23 Aug 2012, 5:59 a.m.,
in response to message #5 by LHH

. Hi, LHH: .

Quote:
This is just the kind of problem I like to try to solve in my head when I go to bed.

Hehe, just in case you're too sleepy to give it a try, this short HP-71B code will find the numeric part of the answer in negligible time:

>CAT
workfile   BASIC   150 07/26/04 15:00

>LIST

10 DESTROY ALL @ INPUT "Max=";M @ R=INT(M^(1/3)) @ DIM S(2*R*R*R)
20 FOR I=1 TO R @ FOR J=I TO R @ N=I*I*I+J*J*J
30 IF S(N) THEN DISP N;"is the minimum" @ END ELSE S(N)=1
40 NEXT J @ NEXT I @ DISP "None found up to";M

>RUN
Max=1000
None found up to 1000

>RUN
Max=2000
 1729 is the minimum

It's a trivial brute-force approach and it certainly could be made faster by simply pre-computing a cubes table up to the required limit so that N=I*I*I+J*J*J would be computed instead as N=C(I)+C(J) and other such tricks but it's pretty pointless as the running time is already completely negligible.

Other trivial optimizations would include:

- taking out of the J loop the invariant computation of I*I*I or C(I) if using the cubes table

- using a string S$ of the correct length instead of a vector filling it up initially to all "0" characters, afterwards storing a "1" to mark each achieved sum in the correct position within the string after previously checking if the position is already marked so that a minimum solution has been found. This would save lots of memory as the string would take much less space than the array (which could have been declared INTEGER to save about 50% memory, by the way).

- if your HP model has the equivalent of flag or bit arrays, this would be ideal, as the S array or the S$ string are just mimicking them

etc, etc, etc, but as it's so trivial and runs instantly I didn't bother, too lazy today ... XD

Best regards from V.

Edited: 23 Aug 2012, 6:11 a.m.

                                    
Re: Programming exercise
Message #7 Posted by Gilles Carpentier on 23 Aug 2012, 10:53 a.m.,
in response to message #6 by Valentin Albillo

HI

10 DESTROY ALL @ INPUT "Max=";M @ R=INT(M^(1/3)) @ DIM S(2*R*R*R)
20 FOR I=1 TO R @ FOR J=I TO R @ N=I*I*I+J*J*J
30 IF S(N) THEN DISP N;"is the minimum" @ END ELSE S(N)=1
40 NEXT J @ NEXT I @ DISP "None found up to";M

Interesting. But how are sure that N is the minimum ?

or:

20 FOR I=1 TO R-1 @ FOR J=I+1 TO R @ N=I*I*I+J*J*J ?

                                    
Re: Programming exercise
Message #8 Posted by C.Ret on 23 Aug 2012, 11:34 a.m.,
in response to message #6 by Valentin Albillo

Hi.

I used the same algorithm as you on my Commodore C128D. In your code, you use the aray S() to "flag" with bi-cubic positions have already bee encounted. That was my first approach idea. But, using a whole integer array to only record '0' and '1' appears to me as a too binary. Instead of the '1' flag, I put the a value (a>0), that makes possible a more explicit display of the results n = a,b = c ,d :

LIST
10 A=0:B%=0:N=0:INPUT "Max:";M%:R%=(M%/2)^(1/3)
20 DIM P%(M%),C(2*R%)
30 FOR N=1 TO 2*R%:C(N)=N*N*N:NEXT N 
40 FOR A=1 TO R:B%=A:DO:B%=B%+1:N=C(A)+C(B%):IF N>M% THEN EXIT
50 IF P%(N) THEN PRINT N,"="A;B%"="P%(N);INT((N-P%(N))^(1/3)):ELSE P%(N)=A
60 LOOP:NEXT A:END

READY.
RUN
Max:? 32000
 1729 = 9  10 = 1  12 
 4104 = 9  15 = 2  16 
 13832 = 18  20 = 2  24 
 20663 = 19  24 = 10  27 
READY.
Note that it takes a few seconds on a Commodore C128/C128D so that pre-computing the cubes (line 30) makes no significant changes. Except perhaps better arithmetic precision, the Commodore 8bits are well known for really bad math! Max limit is 32000 due to conventional memory limit (128 ko). Max limit may be extend using memory add on such as graphic memory or unconventional RAM banks (but also need PEEK/POKE bad programming practice current practices at C64 times!).

Edited: 23 Aug 2012, 11:40 a.m.

                                          
Re: Programming exercise
Message #9 Posted by Gerson W. Barbosa on 23 Aug 2012, 12:30 p.m.,
in response to message #8 by C.Ret

When I see this and other approaches here, mine looks quite silly! It could be done within two loops only once, instead I used two nested loops twice. Anyway, here it is, just for the records: 

%%HP: T(3)A(R)F(,);
\<< 3, XROOT IP \-> r1
  \<< { } 1, r1
    FOR a a r1
      FOR b a SQ a * b SQ b * + +
      NEXT
    NEXT OBJ\-> \->ARRY r1 SQ r1 + 2, / 1, r1 1, - \-> n e1 s1
    \<< 1, 1, n 1, -
      FOR i GETI UNROT i r1 1, + ==
        IF
        THEN s1 'r1' STO+ 1, 'e1' STO+ 's1' 1, STO-
        END e1 r1 s1 \-> e2 r2 s2
        \<< OVER i 1, + DUP n
          FOR j GETI j r2 1, + ==
            IF
            THEN s2 'r2' STO+ 1, 'e2' STO+ 's2' 1, STO-
            END 6, PICK ==
            IF
            THEN 5, PICK 1, IQUOT \->STR "=" + e1 1, IQUOT \->STR + "^3+" + 6, PICK e1 3, ^ - 3, XROOT 1, IQUOT \->STR
              + "^3=" + e2 1,IQUOT \->STR + "^3+" + 6, PICK e2 3, ^ - 3, XROOT 1, IQUOT \->STR + "^3" + 6, ROLLD
            END
          NEXT
        \>> DROP2 ROT DROP
      NEXT DROP
    \>>
  \>> DROP
\>>

40000 TN -->  7: "1729=1^3+12^3=9^3+10^3"
              6: "4104=2^3+16^3=9^3+15^3"
              5: "13832=2^3+24^3=18^3+20^3"
              4: "39312=2^3+34^3=15^3+33^3"
              3: "32832=4^3+32^3=18^3+30^3"
              2: "40033=9^3+34^3=16^3+33^3"
              1: "20683=10^3+27^3=19^3+24^3" 

(after whopping 154 seconds on the emulator!)
                              
Re: Programming exercise
Message #10 Posted by Gerson W. Barbosa on 23 Aug 2012, 12:44 p.m.,
in response to message #5 by LHH

Quote:
Anyway, just doodling I wrote down all the possible cubes and have an answer (BAC-1729) but I did cheat by using a biological computer

Your superior carbon-based computer is absolutely right! I ought to have tried using mine as well :-)

                        
Re: Programming exercise
Message #11 Posted by Gerson W. Barbosa on 23 Aug 2012, 12:49 p.m.,
in response to message #4 by Valentin Albillo

Hello Valentin,

Quote:
Ramanujan, perchance ?

Yes, as we are aware of:

http://en.wikipedia.org/wiki/Taxicab_number

http://mathworld.wolfram.com/TaxicabNumber.html

Best regards,

Gerson.

      
Re: Programming exercise
Message #12 Posted by Thomas Klemm on 23 Aug 2012, 6:04 a.m.,
in response to message #1 by Gerson W. Barbosa

The program is for the HP-42S but can be adapted easily to other models:

STO 00
LBL 00
RCL 00
3
Y^X
STO 02
1
STO 01
LBL 01
3
Y^X
RCL+ 02
STO 03
3
1/X
Y^X
IP
STO 04
RCL 00
1
+
STO 05
LBL 02
RCL 04
RCL 05
X>Y?
GTO 06
3
Y^X
STO 07
1
STO 06
LBL 03
RCL 06
3
Y^X
RCL+ 07
RCL 03
X<=Y?
GTO 04
1
STO+ 06
GTO 03
LBL 04
X#Y?
GTO 05
CLA
AIP
|-"LF"
ARCL 00
|-","
ARCL 01
|-" "
ARCL 05
|-","
ARCL 06
AVIEW
LBL 05
1
STO+ 05
GTO 02
LBL 06
RCL 00
RCL 01
1
+
STO 01
X<Y?
GTO 01
1
STO+ 00
GTO 00

As I don't have the real calculater at hand I can't tell you how fast it is to find the smallest solution. However the result is displayed immediately on my iPhone. Only for the next solutions it takes some time to consider.

Usage: 

1 R/S
R/S
...

Thanks for the nice challenge: I had some fun.

Kind regards
Thomas

Edit: corrected usage

Edited: 23 Aug 2012, 4:13 p.m. after one or more responses were posted

            
Re: Programming exercise
Message #13 Posted by Gerson W. Barbosa on 23 Aug 2012, 1:07 p.m.,
in response to message #12 by Thomas Klemm

Hello Thomas,

Quote:
Thanks for the nice challenge: I had some fun.

I'm glad you have liked it. This is not a challenge however, this is just an exercise. I would pose a challenge only if I had a hard to beat solution, which would be very unlikely here :-)

I tried your program in my HP-42S but I can't get it to work. It's a 107-byte program. I'll double-check the listing later.

Cheers,

Gerson.

                  
Re: Programming exercise
Message #14 Posted by Thomas Klemm on 23 Aug 2012, 3:31 p.m.,
in response to message #13 by Gerson W. Barbosa

Hi Gerson

After adding an END to the program I get the following listing of the labels:

00 { 111-Byte Prgm }
02 LBL 00
09 LBL 01
23 LBL 02
33 LBL 03
44 LBL 04
58 LBL 05
62 LBL 06
73 END

This might help to cross-check with what you have already entered.

It appears the usage I gave was wrong. Instead of 0 start with 1. Apologies for that lapsus. However you should get:

1729
10,9 12,1

4104
15,9 16,2

13832
20,18 24,2

...

Though I've double checked the listing I might still have made some typos. Please forgive me as I'm typing all this on an iPhone which is a pain.

Then once more thanks for the exercise and for taking the time to test my solution. For me it was still a challenge. Part of it was the limitation of the available memory.

Cheers
Thomas

                        
Re: Programming exercise
Message #15 Posted by Gerson W. Barbosa on 23 Aug 2012, 7:40 p.m.,
in response to message #14 by Thomas Klemm

Hello Thomas,

Quote:
It appears the usage I gave was wrong. Instead of 0 start with 1. Apologies for that lapsus. However you should get:

1729
10,9 12,1

4104
15,9 16,2

13832
20,18 24,2

...

The lapsus was mine, sorry! I had forgotten to insert the concatenation characters (shame on me!). Starting with 0 also works. The first solution is found after 32 seconds. The second one is shown about one minute later. Very nice!

Quote:
For me it was still a challenge. Part of it was the limitation of the available memory.

Doing it fast and with little memory is really challenging. When I said it was not a challenge I meant I didn't intend to challenge no one here, rather invite you do do a programming exercise and find better solutions than the clumsy one I had found. My naive method consisted of finding all the possible sums of two cubes and storing them in a list, a second list would store a. When this was done I would search in the first list for repeated sums and take the respective a and c in the second list; b and d would be simply calculated. Because I intended to port the program to the HP-41CV, I replaced the second array with a scheme to keep track of a and b. The QBASIC program below is easier to follow than the RPL version elsewhere, but it still looks awkward :-)

10 CLS
12 INPUT MX
15 R1 = INT(MX ^ (1 / 3))
17 DIM M((R1 * (R1 + 1)) / 2)
20 N = 0
25 FOR A = 1 TO R1
30   FOR B = A TO R1
32     N = N + 1
35     M(N) = A ^ 3 + B ^ 3
50   NEXT B
55 NEXT A
57 E1 = 1: S1 = R1 - 1
60 FOR I = 1 TO N - 1
65   T = M(I)
67   IF I = R1 + 1 THEN R1 = R1 + S1: E1 = E1 + 1: S1 = S1 - 1
68   E2 = E1: R2 = R1: S2 = S1
70   FOR J = I + 1 TO N
71   IF J = R2 + 1 THEN R2 = R2 + S2: E2 = E2 + 1: S2 = S2 - 1
72   IF M(J) <> T THEN 94
74   PRINT USING "#######"; T;
76   PRINT " = ";
78   PRINT USING "##"; E1;
80   PRINT "^3 + ";
82   PRINT USING "##"; (T - E1 ^ 3) ^ (1 / 3);
84   PRINT "^3 = ";
86   PRINT USING "##"; E2;
88   PRINT "^3 + ";
90   PRINT USING "##"; (T - E2 ^ 3) ^ (1 / 3);
92   PRINT "^3"
94   NEXT J
96 NEXT I

? 100000
   1729 =  1^3 + 12^3 =  9^3 + 10^3
   4104 =  2^3 + 16^3 =  9^3 + 15^3
  13832 =  2^3 + 24^3 = 18^3 + 20^3
  39312 =  2^3 + 34^3 = 15^3 + 33^3
  46683 =  3^3 + 36^3 = 27^3 + 30^3
  32832 =  4^3 + 32^3 = 18^3 + 30^3
  40033 =  9^3 + 34^3 = 16^3 + 33^3
  20683 = 10^3 + 27^3 = 19^3 + 24^3
  65728 = 12^3 + 40^3 = 31^3 + 33^3
  64232 = 17^3 + 39^3 = 26^3 + 36^3

The problem in using QBASIC is it much faster than our calculators, real or emulated ones. The ten solutions above are found in about three seconds in my notebook, which disguises a bad algorithm.

Cheers,

Gerson.

            
Re: Programming exercise
Message #16 Posted by Gerson W. Barbosa on 27 Aug 2012, 7:21 p.m.,
in response to message #12 by Thomas Klemm

Thomas,

I tried another approach on Emu71 using no matrices at all. Is it similar to what you have used on Free42? The timing is not accurate (sometimes I get about twice the actual seconds). Thanks!

Cheers,

Gerson.

>LIST
10 DESTROY ALL @ T=TIME @ INPUT N @ R=INT(N^(1/3))
15 FOR A=1 TO R @ A3=A*A*A
20 FOR B=A+1 TO R @ B3=B*B*B
25 FOR C=A+1 TO R @ C3=C*C*C
30 FOR D=C+1 TO R @ D3=D*D*D
35 IF B3-C3=D3-A3 THEN PRINT A3+B3;A;B;C;D
40 NEXT D @ NEXT C @ NEXT B @ NEXT A @ PRINT TIME-T


>RUN
? 1000000
 1729  1  12  9  10
 4104  2  16  9  15
 13832  2  24  18  20
 39312  2  34  15  33
 704977  2  89  41  86
 46683  3  36  27  30
 216027  3  60  22  59
 32832  4  32  18  30
 110656  4  48  36  40
 314496  4  68  30  66
 216125  5  60  45  50
 439101  5  76  48  69
 110808  6  48  27  45
 373464  6  72  54  60
 593047  7  84  63  70
 149389  8  53  29  50
 262656  8  64  36  60
 885248  8  96  72  80
 40033  9  34  16  33
 195841  9  58  22  57
 20683  10  27  19  24
 513000  10  80  45  75
 805688  11  93  30  92
 65728  12  40  31  33
 134379  12  51  38  43
 886464  12  96  54  90
 515375  15  80  54  71
 64232  17  39  26  36
 171288  17  55  24  54
 443889  17  76  38  73
 320264  18  68  32  66
 165464  20  54  38  48
 920673  20  97  33  96
 842751  23  94  63  84
 525824  24  80  62  66
 955016  24  98  63  89
 994688  29  99  60  92
 327763  30  67  51  58
 558441  30  81  57  72
 513856  34  78  52  72
 984067  35  98  59  92
 402597  42  69  56  61
 1016496  47  97  66  90
 1009736  50  96  59  93
 684019  51  82  64  75
 1024.84
                  
Re: Programming exercise
Message #17 Posted by Thomas Klemm on 28 Aug 2012, 3:33 a.m.,
in response to message #16 by Gerson W. Barbosa

Hi Gerson

Quote:
Is it similar to what you have used on Free42?

I don't think so. I'm using the following conditions:

  1. b < a
  2. c > a
  3. c3 <= a3 + b3

These give additional restrictions on c and d.

Here's a Python program that does (more or less) the same:

n = 20
for a in range(n):
    A = a**3
    for b in range(1, a):
        S = A + b**3
        m = trunc(pow(S, 1.0/3))
        for c in range(a+1, m+1):
            C = c**3
            d = 1
            while True:
                T = C + d**3
                if S <= T:
                    break
                d += 1
            if T == S:
                print S, a, b, c, d
    a += 1

And here the HP-42S program with some additional comments:

STO 00             ; input a
LBL 00             ; for a in ...
  RCL 00           
  3                
  Y^X              
  STO 02           ; A = a^3
  1                
  STO 01           ; b = 1
  LBL 01           ; for b in range(1, a)
    3              
    Y^X            
    RCL+ 02        
    STO 03         ; S = A + b^3
    3              
    1/X            
    Y^X            
    IP             
    STO 04         ; m = trunc(pow(S, 1.0/3))
    RCL 00         
    1              
    +              
    STO 05         ; c = a + 1
    LBL 02         ; for c in range(a+1, m+1)
      RCL 04       
      RCL 05       
      X>Y?         
      GTO 06       ; break if c > m
      3            
      Y^X          
      STO 07       ; C = c^3
      1            
      STO 06       ; d = 1
      LBL 03       ; while True
        RCL 06     
        3          
        Y^X        
        RCL+ 07    ; T = C + d^3
        RCL 03     
        X<=Y?      ; break if S <= T
        GTO 04     
        1          
        STO+ 06    ; d += 1
      GTO 03       
      LBL 04       
      X#Y?         ; next c if T # S
      GTO 05       
      CLA          ; if T == S
      AIP          
      |-"LF"       
      ARCL 00      
      |-","        
      ARCL 01      
      |-" "        
      ARCL 05      
      |-","        
      ARCL 06      
      AVIEW        ; print S, a, b, c, d
      LBL 05       
      1            
      STO+ 05      ; c += 1
    GTO 02         
    LBL 06         
    RCL 00         
    RCL 01         
    1              
    +              
    STO 01         ; b += 1
    X<Y?           
  GTO 01           
  1                
  STO+ 00          ; a += 1
GTO 00

HTH
Thomas

                        
Re: Programming exercise
Message #18 Posted by Gerson W. Barbosa on 28 Aug 2012, 8:47 a.m.,
in response to message #17 by Thomas Klemm

Hello Thomas,

Quote:
Quote:
Is it similar to what you have used on Free42?

I don't think so. I'm using the following conditions:

  1. b < a
  2. c > a
  3. c3 <= a3 + b3

Yes, the only similarity might be both algorithms don't need a large amount of memory and thus can be implemented on many programmable calculators. Also, larger numbers can be found given enough time. I estimate I would need about two or three minutes to find the first number on the HP-42S using my algorithm. Way better than my first solution but at least four times slower than yours. Well, the problem is more difficult than I initially though. No wonder it's associated with Hardy, Ramanujan and the like.

It is stated here that 

1729 is the smallest number which can be represented in two different ways as the sum of two cubes:

1729 = 13 + 123 
       93 + 103

The largest known similar number is:

885623890831 = 75113 + 77303 
               87593 + 59783
I think they mean "the largest known primitive solution" (http://oeis.org/A018850), as larger solutions can be found by multiplying the primitive ones by cubes (8, 27, 64, 125...).

Thanks for Python code and the comment to your HP-42S program.

Gerson.

      
Re: Programming exercise
Message #19 Posted by Bill (Smithville, NJ) on 23 Aug 2012, 6:38 a.m.,
in response to message #1 by Gerson W. Barbosa

There's a nice BBC program on this number (part of their Five Numbers series):

The First Taxicab Number

Bill

            
Re: Programming exercise
Message #20 Posted by Gerson W. Barbosa on 23 Aug 2012, 1:21 p.m.,
in response to message #19 by Bill (Smithville, NJ)

Bill, thanks for the link. Although I can listen to BBC live and some recorded programs, to that particular one I cannot. Any other link?

Gerson.

                  
Re: Programming exercise
Message #21 Posted by Bill (Smithville, NJ) on 24 Aug 2012, 7:01 a.m.,
in response to message #20 by Gerson W. Barbosa

Hi Gerson,

I couldn't get the program to play either.

I have the program in MP3 format somewhere.

Send me an email and I'll see if I can find it for you. The whole set of BBC series, Five Numbers, Another Fire Numbers, and A Further Five Numbers are funb to listen to. No heavy math, just fun programs.

Bill

      
Re: Programming exercise
Message #22 Posted by Gilles Carpentier on 23 Aug 2012, 8:39 a.m.,
in response to message #1 by Gerson W. Barbosa

Hi gerson

Here is a program for HP50 wich use GoferList Library http://www.hpcalc.org/details.php?id=6529


 3 INV ^ IP -> Nm
  
  'n^3' 'n' 1 Nm 1 SEQ -> Li
  
   Li 2  DROP Li  NSUB 1 + Nm SUB ADD  DOSUBS 
   Concat DUP Nub Diff SORT

Just type 9999 Plate? (or 99999 ...)

1999 or 9999 takes no time with emu48 100.000 takes 60 sec for the 10 solutions

All solutions are displayed sorted If ypu want only the smallest, add HEAD after the SORT or << MIN >> STREAM in place of SORT

Edited: 23 Aug 2012, 10:29 a.m.

            
Re: Programming exercise
Message #23 Posted by C.Ret on 23 Aug 2012, 1:06 p.m.,
in response to message #22 by Gilles Carpentier

Gilles you use a list to organize the seed of solutions. This gives me the basements of my HP-39gII implementation.

The HP39gII is easy to program and have a really large memory. So I use two list (L0 and L1) to store n=a3+b3 and a respectively. To discover (a b) (c d ) couples I use the convenient POS function which return the position in a list (or 0 if nothing found !).

EXPORT GERSON(m)
BEGIN
LOCAL a,b,c,d,n,r,p;
r:=3NTHROOT(m/2);
L0:={};
L1:={};

FOR a FROM 1 TO r DO
  b:=a;
  REPEAT
    b:=b+1;
    n:= a^3+b^3;
    p:= POS(L0,n);
    IF p>0 THEN
      c:=L1(p);
      d:=3 NTHROOT (n-c^3);
      PRINT(L0(p)+"="+{a,b}+"="+{c,d});
    ELSE
      L0:=CONCAT(L0,{n});
      L1:=CONCAT(L1,{a});
    END;
  UNTIL n>m ;
END;
END;

Please note the elegant and powerful NTHROOT function of this HP-39gII !

EDIT : Sorry for this confusing PRINTing : A better solution is : 

PRINT(n+"="+a+"+"+b+"="+c+"+"+d+"");
And you get a readable output full screen:

Edited: 23 Aug 2012, 1:38 p.m.

                  
Re: Programming exercise
Message #24 Posted by Tim Wessman on 23 Aug 2012, 2:05 p.m.,
in response to message #23 by C.Ret

Your program running on the real hardware 39gII takes ~7.3s for an input of 100000 returning 10 results.

Was that time for emu48 running real calc speed?

TW

Edited: 23 Aug 2012, 2:06 p.m.

                        
Re: Programming exercise
Message #25 Posted by Gilles Carpentier on 23 Aug 2012, 4:36 p.m.,
in response to message #24 by Tim Wessman

Hi Tim,

the algoritm is not the same but on real 50G hdw (the EMU48 is not accurate at all for 50G program) i get for all solutions :

1999 : 3.3 s
9999 : 16.3 s
99999 : 1504 s :O

So the 39gII is ~200 times faster here !. I will try to translate the CRET prog on my 50G for an accurate benchmark

I estimate in general the 39gII language 40 to 80 times faster than the UserRPL/50 depending what you do (and even far more quick with recursive program wich are incredibly slow in user RPL)

For sure the 39gII speed is impressive (and use almost no battery !)

Edited: 23 Aug 2012, 5:14 p.m.

                              
Re: Programming exercise
Message #26 Posted by C.Ret on 24 Aug 2012, 3:04 a.m.,
in response to message #25 by Gilles Carpentier

Here is a direct translation of the HP39gII version for HP28S. You still have to found any tricks to make it run faster on HP50g.

As for the HP39gII, the two large lists are store in main (global) memory. The local variables are exactly the same { m r a b n p c d } but are only set in separate nested sub-structures:

 DUP 2 / 3 INV ^ -> m r
   {} 'L0' STO
    {} 'L1' STO
    1 r FOR a
      a
      DO
        1 +
        DUP 3 ^ a 3 ^ +
        L0 OVER POS     -> b n p
         IF p THEN
            'L1' p GET
            n OVER 3 ^ - 3 INV ^ -> c d
             n a b c d 5 ->LIST DUP 1 DISP 
          ELSE
            L0 n + 'L0' STO
            L1 a + 'L1' STO
          END
          b n
        
      UNTIL m > END
      DROP
    NEXT
    CLMF
  

'GERSON' STO

2000  [GERSON]   ---->  1:{ 1729 9 10 1 12 }                in  ~1'15"

                        4:{ 1729 9 10 1 12 }
                        3:{ 4104 9 15 2 16 }
                        2:{ 13832 18 20 2 23.9999999999 }
20000 [GERSON]   ---->  1:{ 20683 19 24 10 26.9999999999 }  in  ~39'37" on real HP28S

Edited: 24 Aug 2012, 4:03 a.m.

            
Re: Programming exercise
Message #27 Posted by Thomas Klemm on 23 Aug 2012, 4:07 p.m.,
in response to message #22 by Gilles Carpentier

Just for the record: a solution for the HP-48GX without additional libraries:

\<< 3 XROOT IP \-> n
  \<< { } DUP 1 n
    FOR i 1 i
      FOR j
        i SQ i *
        j SQ j *
        + DUP2
        IF POS 0 >
        THEN
          ROT SWAP + SWAP
        ELSE
          +
        END
      NEXT
    NEXT
    DROP
  \>>
\>>

It's reasonable fast for 1999 or 9999 but gets rather slow for bigger values like 105. My emulator m48 just tells me: DON'T PANIC. But I wasn't patient enough to wait for the result.

Cheers
Thomas

                  
Re: Programming exercise
Message #28 Posted by Gilles Carpentier on 23 Aug 2012, 5:56 p.m.,
in response to message #27 by Thomas Klemm

Another way, use the Valentin Albillo algorithm but all with stack and modify to show all solutions :

 
 3 XROOT IP DUP DUPDUP * * 2 * O {}  -> r t n s
 
  0. t NDUPN DROP
  1. r FOR a
   a  r FOR b
    a DUPDUP * * b DUPDUP * * + DUP 'n' STO
    IF PICK THEN n 's' STO+ ELSE 1. n UNPICK END
  NEXT
 NEXT
 t DROPN s

9999 takes 7,3 s on emulator at real speed (probably less with a real 50G cal)

99999 results 'out of memory'

Use a string instead the stack could be interesting

Edited: 23 Aug 2012, 6:07 p.m.

                        
Re: Programming exercise
Message #29 Posted by Thomas Klemm on 24 Aug 2012, 5:04 p.m.,
in response to message #28 by Gilles Carpentier

Quote:
Use a string instead the stack could be interesting

Not exactly a string but the HP-42S allows for bit operations. So we can squeeze 36 bits into one register which is enough to find the smallest solutions. 

00 { 65 Byte Prgm }
   CLRG
   LBL 00
   RCL 00
   STO 01
05 ENTER
   X^2
   *
   STO 02
   LBL 01
10 RCL 01
   ENTER
   X^2
   *
   RCL+ 02
15 CLA
   AIP
   36
   RCL ST Y
   RCL ST Y
20 BASE/
   3
   +
   RCL IND ST X
   RDN
25 RDN
   MOD
   BIT?
   GTO 02
   +/-
30 1
   X<>Y
   ROTXY
   OR
   STO IND ST Y
35 GTO 03
   LBL 02
   AVIEW
   LBL 03
   DSE 01
40 GTO 01
   1
   STO+ 00
   GTO 00
44 END

Since Free42 allows to set the size to 9999 it's possible to find even more. But of course it's still rather limited compared to other calculators.

Cheers
Thomas

                  
Re: Programming exercise
Message #30 Posted by C.Ret on 24 Aug 2012, 2:58 a.m.,
in response to message #27 by Thomas Klemm

I like Thomas Klemm's version. Consise, good use of the stack and the lists. It also work like a charm on my old HP28S ! I get the { 1729 } answer alter about 1 min on it.

I have to admit that I am a bit lost with Gille's version despite they are both written in user RPL and look to follow relatively close algorithms.

Edited: 24 Aug 2012, 4:11 a.m.

                        
Re: Programming exercise
Message #31 Posted by Gilles Carpentier on 24 Aug 2012, 12:07 p.m.,
in response to message #30 by C.Ret

Cret, Some comments :

@ Get the cubic rooth of the max numbre

   3 INV ^ IP -> Nm 

@ Create Li : the list of all the cube < max number

   'n^3' 'n' 1 Nm 1 SEQ -> Li

@ For each element of the list
@   Get the tail of the list after the element position
@   Use ADD. 5 { 1 2 3 4 } ADD -> { 6 7 8 9 }
@   At the end of process, we get a list of list :
@   { { first add with all following}{ second add with...} ... { sum of the 2 last} }
@   The number 2 (Li 2) is to avoid an error for the last item.
@   With '2' DOSUBS works with 2 elements each time and then stop to the antepenultimate element

  Li 2  DROP Li  NSUB 1 + Nm SUB ADD  DOSUBS

@ 'Concat' the list so that { { a b c }{ d e }} -> { a b c d e }

   Concat  @ here we have a list with _all_ the a^3+b^3 combinaison
   DUP     @ but we want only the duplicates items ! so :   

@ Delete the duplicate of one list. The initial list is kept on the stack

  Nub    @ Ex { 1 4 5 4 1 } Nub -> {5}

@ Get all the duplicates by the difference of the 2 lists 
@ (there if no command "get the duplicate of a list" in Gofer
@ but the sequance <DUP Nub Diff> does the job 

  Diff     @ Ex { 1 4 5 4 1 }  { 5 } Diff -> { 1 4 }

  SORT     @ Here we get all the solutions sorted

I use a lot of this kind of programmation with Eulers projects.

in fact with global variables and few change, it s just 3 simple lines of code and few bytes 

3 XROOT IP 'Nm' STO 
1 Nm 1 Seq 3 ^ 'Li' STO
Li 2  DROP Li  NSUB 1 + Nm SUB ADD  DOSUBS
Concat DUP Nub Diff SORT

and i think this 3 lines could work :

1 3 XROOT IP 1 Seq 3 ^ 'Li' STO
Li 2  DROP Li  NSUB 1 + MAXR SUB ADD  DOSUBS
Concat DUP Nub Diff SORT

Edited: 24 Aug 2012, 12:46 p.m.

      
Re: Programming exercise
Message #32 Posted by Les Koller on 23 Aug 2012, 7:20 p.m.,
in response to message #1 by Gerson W. Barbosa

Shades of Ramanujan and Hardy?

            
Re: Programming exercise
Message #33 Posted by Gerson W. Barbosa on 23 Aug 2012, 7:42 p.m.,
in response to message #32 by Les Koller

Yes! :-)

http://en.wikipedia.org/wiki/Taxicab_number

http://mathworld.wolfram.com/TaxicabNumber.html

[ Return to Index | Top of Index ]

Go back to the main exhibit hall