Sums of cubes (was Re: histrionics)

["Robert I. Eachus" <eachus@mitre.org>]

Ehud Lamm wrote:
> 
> If we are on the subject of VHLLs and SETL, let me ask a SETL question. I
> really don't knwo the language, only the basic idea - and since I can not
> find a PC version, it'll take some time for me to get to know it better
> (but learning something Dewar was involved in is likely to be a pleasure:
> SNOBOL, GNAT etc.)

    Good point.  It might be a nice idea to implement SETL in Ada so it
is
easily portable.  (This would be going full circle, since the first
validated
Ada implementation was in SETL!)
 
> How would you find all numbers like 407=4^3+0^3+7^3 that are equal to the
> sum of their digits' cubes?

    Interesting question, but I have found that programming in Ada works
best with a "sandwich" approach, or maybe you want to call it meet in
the middle.  To illustrate, note that one operation that we will be
doing a lot is cubing numbers.  This can be implemented in the form of a
table:

   procedure Sum_of_Cubes is 
     Cube_Table: array(0..9) of Integer :=
       (0, 1, 8, 27, 64, 125, 216, 343, 512, 729);
     -- There are other ways to initialize this table, but cut and paste
from
     -- APL and add commas is pretty quick.
     
     function Cubes(N: Integer) return Integer
       is begin return Cube_Table(N); end Cubes;
     pragma Inline(Cube);

   begin
     TBD;
   end Sum_of_Cubes;

   Okay, now to the algorithm.  There is no need to generalize to more
digits yet, and it is simple enough to write three loops.  I'll keep the
convention
that 000 and 001 are three digit numbers...

   procedure Sum_of_Cubes is
     Cube_Table: array(0..9) of Integer :=
       (0, 1, 8, 27, 64, 125, 216, 343, 512, 729);
     -- There are other ways to initialize this table, but cut and paste
from
     -- APL is very quick.
     
     function Cube(N: Integer) return Integer
       is begin return Cube_Table(N); end Cube;
     pragma Inline(Cube);

   begin
     for I in 0..9 loop
       for J in 0..9 loop
         for K in 0..9 loop
           if Cube(I) + Cube(J) + Cube(K) = 100*I + 10*J + K
           then Print_Result;
           end if;
         end loop;
       end loop;
     end loop;
   end Sum_of_Cubes;

   Now back to the low level, the computation on the right side of the
equality
can be strenth reduced.  A compiler might be smart enough to do it, but
there is no reason not to "do it ourselves" since the algorithm becomes
cleaner, IMHO.
 
   procedure Sum_of_Cubes is
     Cube_Table: array(0..9) of Integer :=
       (0, 1, 8, 27, 64, 125, 216, 343, 512, 729);
     Number: Integer := 0;     
     function Cube(N: Integer) is
       begin return Cube_Table(N); end Cube;
     pragma Inline(Cube);

   begin
     for I in 0..9 loop
       for J in 0..9 loop
         for K in 0..9 loop
           if Cube(I) + Cube(J) + Cube(K) = Number;
           then Print_Result;
           end if;
           Number := Number + 1;
         end loop;
       end loop;
     end loop;
   end Sum_of_Cubes;

   Should we do the same optimization for the sums of cubes?  It looks
like the
best solution/compromise is to do half of it:

   procedure Sum_of_Cubes is
     Cube_Table: array(0..9) of Integer :=
       (0, 1, 8, 27, 64, 125, 216, 343, 512, 729);
     Number: Integer := 0;
     Partial: Integer;     
     function Cube(N: Integer)return Integer
       is begin return Cube_Table(N); end Cube;
     pragma Inline(Cube);

   begin
     for I in 0..9 loop
       for J in 0..9 loop
         Partial := Cube(I) + Cube(J);
         for K in 0..9 loop
           if  Partial + Cube(K) = Number
           then Print_Result;
           end if;
           Number := Number + 1;
         end loop;
       end loop;
     end loop;
   end Sum_of_Cubes;

   Next step write Print_Result.  This could be done as an inline
procedure,
but it is now simple enough to do in place:

   with Ada.Integer_Test_IO; with Ada.Text_IO;
   procedure Sum_of_Cubes is
      subtype Digit is Integer range 0..9;
      Cube_Table: array(Digit) of Integer :=
           (0, 1, 8, 27, 64, 125, 216, 343, 512, 729);
     Number: Integer := 0;
     Partial: Integer;     
     function Cube(N: Integer) return Integer is
        begin return Cube_Table(N); end Cube;
     pragma Inline(Cube);

   begin
     for I in 0..9 loop
       for J in 0..9 loop
         Partial := Cube(I) + Cube(J);
         for K in 0..9 loop
           if  Partial + Cube(K) = Number;
           then Ada.Integer_Text_IO.Put(Number, 4);
           end if;
           Number := Number + 1;
         end loop;
       end loop;
     end loop;
   Ada.Text_IO.New_Line;
   end Sum_of_Cubes;

   A last cleanup that might help the compiler, but will certainly help
readers and maintainers.  There are too many instances of 0..9 in the
code.  Time for a subtype:

   with Ada.Integer_Text_IO; with Ada.Text_IO;
   procedure Sum_of_Cubes is
     subtype Digit is Integer range 0..9;
     Cube_Table: array(Digit) of Integer :=
           (0, 1, 8, 27, 64, 125, 216, 343, 512, 729);
     Number: Integer := 0;
     Partial: Integer;     
     function Cube(N: Digit) return Integer is
        begin return Cube_Table(N); end Cube;
     pragma Inline(Cube);

   begin
     for I in Digit loop
       for J in Digit loop
         Partial := Cube(I) + Cube(J);
         for K in Digit loop
           if  Partial + Cube(K) = Number
           then Ada.Integer_Text_IO.Put(Number, 4);
           end if;
           Number := Number + 1;
         end loop;
       end loop;
     end loop;
   Ada.Text_IO.New_Line;
   end Sum_of_Cubes;

   Modifying to check four digit numbers results in no new matches, but
is very easy (two new lines, two modified, if you don't count the name. 
Also, execution doesn't take long enough to measure.  Changing from
cubes to fourth power (but keeping four digits) results in an output of
    0    1 1634 8208 9474
and only requires changing one line (again ignoring name changes...) 
-- 

                                        Robert I. Eachus

with Standard_Disclaimer;
use  Standard_Disclaimer;
function Message (Text: in Clever_Ideas) return Better_Ideas is...