[Shootout-list] Statistics Ada bench

[Pascal Obry]

Here is a fixed version of the statistic bench for Ada.

-- $Id: moments-gnat.code,v 1.7 2005/03/15 06:18:00 bfulgham Exp $
-- http://dada.perl.it/shootout/
-- Ada 95 code by C.C.
-- Modified by Pascal Obry on 2005/03/18

-- Annotated Ada Reference Manual ISO/IEC 8652:1995: http://www.ada-auth.org/

with System, Ada.Numerics.Generic_Elementary_Functions;
with Ada.Unchecked_Deallocation;
with Ada.Text_IO; use Ada.Text_IO;

procedure Moments is

   type Real is digits Positive'Max (15, System.Max_Digits);

   package AF is new Ada.Numerics.Generic_Elementary_Functions (Real);
   package RIO is new Ada.Text_IO.Float_IO (Num => Real);

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

   generic
       type Item_Type is private;
       with function "<=" (X, Y : Item_Type) return Boolean;
       type Sequence is array (Integer range <>) of Item_Type;
   package Sort_Pkg is
      procedure Quick_Sort (S : in out Sequence);
   end Sort_Pkg;  --  Copied from Southampton Ada MPICH supercomputer bindings

   package body Sort_Pkg is

      procedure Quick_Sort (S : in out Sequence) is

         procedure Sort_Recursive (Lwb, Upb : Integer) is
            Pivot : Item_Type := S (Upb);
            Front : Integer := Lwb;
            Back  : Integer := Upb;
            Temp  : Item_Type;
         begin
            if Lwb < Upb then

               while Front <= Back loop

                  while not (Pivot <= S (Front)) loop
                     Front := Front + 1;
                  end loop;

                  while not  (S (Back) <= Pivot) loop
                     Back := Back - 1;
                  end loop;

                  if Front <= Back then
                     Temp := S (Front);
                     S (Front) := S (Back);
                     S (Back) := Temp;
                     Front := Front + 1;
                     Back  := Back - 1;
                  end if;
               end loop;

               Sort_Recursive (Lwb, Back);
               Sort_Recursive (Front, Upb);
            end if;
         end Sort_Recursive;

      begin
         Sort_Recursive (S'First, S'Last);
      end Quick_Sort;

   end Sort_Pkg;

   type Real_Array is array (Integer range <>) of Real;
   type Real_Array_Ptr is access Real_Array;

   procedure Deallocate_Real_Array is new Ada.Unchecked_Deallocation
     (Object => Real_Array, Name => Real_Array_Ptr);

   package Sort is new Sort_Pkg (Real, "<=" => "<=", Sequence => Real_Array);

   Data                     : Real_Array_Ptr := new Real_Array (1 .. 4_096);
   Temp_Array               : Real_Array_Ptr;
   Dev, D_2, Mean, Median   : Real;
   Standard_Deviation       : Real;
   Sum, Avg_Abs_Deviation   : Real    := 0.0;
   Variance, Skew, Kurtosis : Real    := 0.0;
   M                        : Natural := 0;

begin
   while not End_Of_File loop
      M := M + 1;

      if M > Data'Last then
         Temp_Array := new Real_Array (1 .. Data'Last * 4);
         Temp_Array (Data'Range) := Data.all;
         Deallocate_Real_Array (Data);
         Data := Temp_Array;
      end if;

      RIO.Get (Item => Data (M));
      Sum := Sum + Data (M);
   end loop;

   Mean := Sum / Real (M);

   for K in 1 .. M loop
      Dev               := Data (K) - Mean;
      Avg_Abs_Deviation := Avg_Abs_Deviation + abs Dev;
      D_2               := Dev * Dev;
      Variance          := Variance + D_2;
      Skew              := Skew + (D_2 * Dev);
      Kurtosis          := Kurtosis + (D_2 * D_2);
   end loop;

   Avg_Abs_Deviation  := Avg_Abs_Deviation / Real (M);
   Variance           := Variance / Real (M - 1);
   Standard_Deviation := AF.Sqrt (Variance);

   if Variance < 10.0 * Real'Model_Epsilon then
      Put_Line ("> Reduced accuracy results: 0 = ((Variance/10 + 1) - 1)");
   else
      Skew := Skew / (Real (M) * Variance * Standard_Deviation);
      Kurtosis := -3.0 + Kurtosis / (Real (M) * Variance * Variance);
   end if;

   Sort.Quick_Sort (S => Data (1 .. M));

   if 1 = (M mod 2) then
      Median := Data ((M + 1) / 2);
   else
      Median := (Data (M / 2) + Data (1 + M / 2)) / 2.0;
   end if;

   Put_Line ("n:                 " & Integer'Image (M));
   Put ("median:             "); Put (Median);              New_Line;
   Put ("mean:               "); Put (Mean);                New_Line;
   Put ("average_deviation:  "); Put (Avg_Abs_Deviation);   New_Line;
   Put ("standard_deviation: "); Put (Standard_Deviation);  New_Line;
   Put ("variance:           "); Put (Variance);            New_Line;
   Put ("skew:               "); Put (Skew);                New_Line;
   Put ("kurtosis:           "); Put (Kurtosis);            New_Line;
end Moments;

Pascal.

-- 

--|------------------------------------------------------
--| Pascal Obry                           Team-Ada Member
--| 45, rue Gabriel Peri - 78114 Magny Les Hameaux FRANCE
--|------------------------------------------------------
--|              http://www.obry.org
--| "The best way to travel is by means of imagination"
--|
--| gpg --keyserver wwwkeys.pgp.net --recv-key C1082595