The Q6600
Benchmarks Game

n-body Ada 2012 GNAT #5 program

source code

-- The Computer Language Benchmarks Game
-- https://salsa.debian.org/benchmarksgame-team/benchmarksgame/
--
-- Contributed by Pascal Obry on 2005/03/21
-- Modified by Brian Drummond on 2011/03/24
-- Updated by Jonathan Parker and Georg Bauhaus (May 2012)


with Ada.Command_Line; use Ada.Command_Line;
with Ada.Text_IO;      use Ada.Text_IO;
with Nbody_Pck;        use Nbody_Pck;
with Root;

procedure Nbody is

   subtype Real is Root.S_Real;

   package RIO is new Float_Io (Real);

   procedure Put
     (Item : Real; Fore : Field := 0; Aft : Field := 9;
      Exp  : Field := 0) renames RIO.Put;

   N : constant Integer := Integer'Value (Argument (1));

   Px, Py, Pz : Real := 0.0;

begin
   for I in Body_Name'Range loop
      Add_Momentum (I, Px, Py, Pz);
   end loop;

   Offset_Momentum (Sun, Px, Py, Pz);

   Put (Energy);
   New_Line;

   for K in 1 .. N loop
      Advance (0.01);
   end loop;

   Put (Energy);
   New_Line;

end Nbody;

with Ada.Numerics; use Ada.Numerics;
with Root; use Root;

package Nbody_Pck is

   subtype Real is Root.S_Real;

   Solar_Mass    : constant Real := 4.0 * Pi * Pi;
   Days_Per_Year : constant Real := 365.24;

   type Signed is range -2**15 .. 2**15-1; 
   subtype Body_Name is Signed range 0 .. 4;

   Jupiter : constant := 0;
   Saturn  : constant := 1;
   Neptune : constant := 2;
   Uranus  : constant := 3;
   Sun     : constant := 4;

   type Axes is (X, Y, Z);

   procedure Offset_Momentum
     (Planet     : in Body_Name;
      Px, Py, Pz : in Real);

   procedure Add_Momentum
     (Planet     : in     Body_Name;
      Px, Py, Pz : in out Real);

   function Energy return Real;

   procedure Advance (Dt : in Real);

private

   type Solar_System is array (Body_Name, Axes) of Real;
   pragma Convention (Ada, Solar_System);

   Position : Solar_System :=
     (Jupiter => (X  =>  4.84143144246472090e+00,
                  Y  => -1.16032004402742839e+00,
                  Z  => -1.03622044471123109e-01),
      Saturn  => (X  =>  8.34336671824457987e+00,
                  Y  =>  4.12479856412430479e+00,
                  Z  => -4.03523417114321381e-01),
      Uranus  => (X  =>  1.28943695621391310e+01,
                  y  => -1.51111514016986312e+01,
                  Z  => -2.23307578892655734e-01),
      Neptune => (X  =>  1.53796971148509165e+01,
                  Y  => -2.59193146099879641e+01,
                  Z  =>  1.79258772950371181e-01),
      Sun     => (X  =>  0.0,
                  Y  =>  0.0,
                  Z  =>  0.0));

   Velocity : Solar_System :=
     (Jupiter => (X  =>  1.66007664274403694e-03 * Days_Per_Year,
                  Y  =>  7.69901118419740425e-03 * Days_Per_Year,
                  Z  => -6.90460016972063023e-05 * Days_Per_Year),
      Saturn  => (X  => -2.76742510726862411e-03 * Days_Per_Year,
                  Y  =>  4.99852801234917238e-03 * Days_Per_Year,
                  Z  =>  2.30417297573763929e-05 * Days_Per_Year),
      Uranus  => (X  =>  2.96460137564761618e-03 * Days_Per_Year,
                  Y  =>  2.37847173959480950e-03 * Days_Per_Year,
                  Z  => -2.96589568540237556e-05 * Days_Per_Year),
      Neptune => (X  =>  2.68067772490389322e-03 * Days_Per_Year,
                  Y  =>  1.62824170038242295e-03 * Days_Per_Year,
                  Z  => -9.51592254519715870e-05 * Days_Per_Year),
      Sun     => (X  =>  0.0,
                  Y  =>  0.0,
                  Z  =>  0.0));

   type Body_Mass is array(Body_Name) of Real;

   Mass : constant Body_Mass :=
     (Jupiter => 9.54791938424326609e-04 * Solar_Mass,
      Saturn  => 2.85885980666130812e-04 * Solar_Mass,
      Uranus  => 4.36624404335156298e-05 * Solar_Mass,
      Neptune => 5.15138902046611451e-05 * Solar_Mass,
      Sun     => Solar_Mass);

end Nbody_Pck;

package body Nbody_Pck is

   procedure Offset_Momentum
     (Planet     : in Body_Name;
      Px, Py, Pz : in Real) is
   begin
      Velocity (Planet, X) := -Px / Solar_Mass;
      Velocity (Planet, Y) := -Py / Solar_Mass;
      Velocity (Planet, Z) := -Pz / Solar_Mass;
   end Offset_Momentum;

   procedure Add_Momentum
     (Planet     : in     Body_Name;
      Px, Py, Pz : in out Real) is
   begin
      Px := Px + Velocity (Planet, X) * Mass (Planet);
      Py := Py + Velocity (Planet, Y) * Mass (Planet);
      Pz := Pz + Velocity (Planet, Z) * Mass (Planet);
   end Add_Momentum;

   function Energy return Real is
      Dx, Dy, Dz, Distance : Real;
      E                    : Real := 0.0;
   begin
      for i in Body_Name loop
        E := E + 0.5 * Mass (i) *
          (Velocity (i, X) * Velocity (i, X)
         + Velocity (i, Y) * Velocity (i, Y)
         + Velocity (i, Z) * Velocity (i, Z));

        if i /= Body_Name'Last then
           for j in Body_Name'Succ (i) .. Body_Name'Last loop
              Dx := Position (i, X) - Position (j, X);
              Dy := Position (i, Y) - Position (j, Y);
              Dz := Position (i, Z) - Position (j, Z);
              Distance := Sqrt (Dx * Dx + Dy * Dy + Dz * Dz);
              E := E - (Mass (i) * Mass (j)) / Distance;
           end loop;
        end if;
      end loop;
      return E;
   end Energy;

   procedure Advance (Dt : in Real) is
      Dx, Dy, Dz, Mag, s : Real;
      Mass_i, Mass_j : Real;
   begin
      for i in Body_Name loop
         Mass_i := Mass(i);
         for j in Body_Name loop
            if j > i then
               Dx := Position (i, X) - Position (j, X);
               Dy := Position (i, Y) - Position (j, Y);
               Dz := Position (i, Z) - Position (j, Z);

               Mass_j := Mass(j);
               s      := SSE_Reciprocal_Sqrt (Dx*Dx + Dy*Dy + Dz*Dz);
               Mag    := (dt * s) * (s * s);

               Velocity (i, X) := Velocity (i, X) - Dx * Mass_j * Mag;
               Velocity (j, X) := Velocity (j, X) + Dx * Mass_i * Mag;
               Velocity (i, Y) := Velocity (i, Y) - Dy * Mass_j * Mag;
               Velocity (j, Y) := Velocity (j, Y) + Dy * Mass_i * Mag;
               Velocity (i, Z) := Velocity (i, Z) - Dz * Mass_j * Mag;
               Velocity (j, Z) := Velocity (j, Z) + Dz * Mass_i * Mag;
            end if;
         end loop;
      end loop;

      for i in Body_Name loop
         Position (i, X) := Position (i, X) + Dt * Velocity (i, X);
         Position (i, Y) := Position (i, Y) + Dt * Velocity (i, Y);
         Position (i, Z) := Position (i, Z) + Dt * Velocity (i, Z);
      end loop;
   end Advance;

end Nbody_Pck;

package Root is

   type S_Real is new Long_Float;

   pragma Assert (S_Real'Size = 64 and S_Real'digits > 13);
  
   type SSE_Vector is array (0 .. 1) of S_Real;

   function Sqrt (X : S_Real) return S_Real;

   function Sqrt_of_Reciprocal (X : S_Real) return S_Real;

   function SSE_Reciprocal_Sqrt (X : S_Real) return S_Real;
   -- Returns double precision 1.0 / Sqrt(X), for Long_Float X.

   pragma Inline (SSE_Reciprocal_Sqrt);

end Root;

package body root is

   -- "divpd" and "sqrtpd" are  double precision:

   type m128d is array (0 .. 1) of S_Real;
   for m128d'Alignment use 16;
   pragma Machine_Attribute (m128d, "vector_type");

   function ia32_Div (X, Y : m128d) return m128d;
   pragma Import (Intrinsic, ia32_Div, "__builtin_ia32_divpd");

   function ia32_Sqrt (X : m128d) return m128d;
   pragma Import (Intrinsic, ia32_Sqrt, "__builtin_ia32_sqrtpd");

   function Sqrt (X : S_Real) return S_Real is
   begin return ia32_Sqrt ((X, 1.0))(0);
   end Sqrt;

   function Sqrt_of_Reciprocal (X : S_Real) return S_Real is
      a : constant m128d := ia32_Div ((1.0, 1.0), (X, 1.0));
      b : constant m128d := ia32_Sqrt (a);
   begin
      return b(0);
   end Sqrt_of_Reciprocal;

   -- "rsqrtps" (Reciprocal Sqrt) operates on Float (single precision):

   type m128s is array (0 .. 3) of Float;
   for m128s'Alignment use 16;
   pragma Machine_Attribute (m128s, "vector_type");
   pragma Assert (Float'Digits < 7 and m128s'size = 128);

   function ia32_RSqrt (X : m128s) return m128s;
   pragma Import (Intrinsic, ia32_RSqrt, "__builtin_ia32_rsqrtps");

   function Recip_Sqrt (X : S_Real) return S_Real is
      Z :  constant m128s := ia32_RSqrt ((Float (X), others => 1.0));
      r : S_Real := S_Real (Z(0));
   begin
      for i in 1 .. 2 loop
         r := r * 1.5 - ((0.5 * X) * r) * (r * r);
      end loop;
      return r;
   end Recip_Sqrt;
   pragma Inline (Recip_Sqrt); 

   function SSE_Reciprocal_Sqrt (X : S_Real) return S_Real is
   begin
      if Abs X < 1.0e+30 and Abs X > 1.0e-30 then
         return Recip_Sqrt (X);
      else
         return Sqrt_of_Reciprocal (X);
      end if;
   end SSE_Reciprocal_Sqrt;

   x : constant m128d := (4.0, 6.0);
   y : constant m128d := (2.0, 2.0);
   pragma Assert (ia32_Div(x, y) = m128d'(2.0, 3.0));
   -- Minimal test, but a good idea when using pragma Import.

end root;


    

notes, command-line, and program output

NOTES:
64-bit Ubuntu quad core
GNATMAKE 9.3.0
gcc (Ubuntu 9.3.0-10ubuntu2) 9.3.0



Mon, 04 May 2020 19:18:11 GMT

MAKE:
gnatchop -r -w nbody.gnat-5.gnat
splitting nbody.gnat-5.gnat into:
   nbody.adb
   nbody_pck.ads
   nbody_pck.adb
   root.ads
   root.adb
gnatmake -O3 -fomit-frame-pointer -march=core2 -gnatNp -f nbody.adb -o nbody.gnat-5.gnat_run 
x86_64-linux-gnu-gcc-9 -c -O3 -fomit-frame-pointer -march=core2 -gnatNp nbody.adb
x86_64-linux-gnu-gcc-9 -c -O3 -fomit-frame-pointer -march=core2 -gnatNp nbody_pck.adb
x86_64-linux-gnu-gcc-9 -c -O3 -fomit-frame-pointer -march=core2 -gnatNp root.adb
x86_64-linux-gnu-gnatbind-9 -x nbody.ali
x86_64-linux-gnu-gnatlink-9 nbody.ali -O3 -fomit-frame-pointer -march=core2 -o nbody.gnat-5.gnat_run

4.90s to complete and log all make actions

COMMAND LINE:
./nbody.gnat-5.gnat_run 50000000

PROGRAM OUTPUT:
-0.169075164
-0.169059907