with Report;
with Ada.Real_Time;
with Text_IO;

procedure Whetstone is

   use Report;
   use Ada.Real_Time;
   use Text_IO;

   package IO is new Integer_IO (Integer);
   use IO;

   Start_Time, Stop_Time : Time;

   -- 
   -- Perform one million Whetstone instructions per cycle
   -- 

   no_of_cycles : natural := 1;   -- number of cycles of test to be timed
   rating : natural;              -- number of K Whetstones per sec


   i : integer := 10;     -- loop count weighting factor;

 
   type vector is array (integer range <>) of float;
   x1,x2,x3,x4,x,y,z,t,t1,t2 : float;
   e1 : vector (1..4);
   j,k,l,n1,n2,n3,n4,n5,n6,n7,n8,n9,n10,n11 : integer;
 

-- This is a standard Ada inplementation of the required math routines
-- that is Copyright Westinghouse Electric Corporation 1983,1984,1985.
-- These routines are provided for use by ACM SIGAda PIWG for making
-- measurements. These routines are copyrighted and may only be used for
-- performance measurements. These math routines must not be distributed
-- without this notice. No permission is granted to any party to modify,
-- to redistribute, to sell, to give away, or to otherwise use or
-- transmit these math routines without express written permission from
-- Westinghous Electric Corporation, c/o Jon Squire, P.O. Box 746 MS1615,
-- Baltimore, MD 21203.

  PI_2 : constant FLOAT := 1.5707963267949 ;
  PI : constant FLOAT := 2.0 * PI_2 ;

  function SIN ( X : FLOAT ) return FLOAT is -- Copyright Westinghouse 1985

    C1 : constant FLOAT := 1.57079631847 ;
    C3 : constant FLOAT := - 0.64596371106 ;
    C5 : constant FLOAT := 0.07968967928 ;
    C7 : constant FLOAT := - 0.00467376557 ;
    C9 : constant FLOAT := 0.00015148419 ;
    X_NORM : FLOAT ;
    X_INT : FLOAT ;
    X_2 : FLOAT ;
    Y : FLOAT ;
  begin
    X_NORM := X / PI_2 ;
    if abs ( X_NORM ) > 4.0 then  -- REDUCE TO  -2 PI .. 2 PI
      X_INT := FLOAT ( INTEGER( X_NORM / 4.0 )) ;
      X_NORM := X_NORM - 4.0 * X_INT ;
    end if ;
    if X_NORM > 2.0 then  -- REDUCE TO  -PI .. PI
      X_NORM := 2.0 - X_NORM ;
    elsif X_NORM < - 2.0 then
      X_NORM := - 2.0 - X_NORM ;
    end if ;
    if X_NORM > 1.0 then  -- REDUCE TO  -PI/2 .. PI/2
      X_NORM := 2.0 - X_NORM ;
    elsif X_NORM < - 1.0 then
      X_NORM := - 2.0 - X_NORM ;
    end if ;
    X_2 := X_NORM * X_NORM ;
    Y := ( C1 +( C3 +( C5 +( C7 + C9 * X_2 ) * X_2) * X_2) * X_2) * X_NORM ;
    return Y ;
  end SIN ;

  function COS ( X : FLOAT ) return FLOAT is

  begin
    return SIN ( X + PI_2 ) ;
  end COS ;


  function ATAN ( X : FLOAT ) return FLOAT is -- Copyright Westinghouse 1985

    C1 : constant FLOAT := 0.9999993329 ;
    C3 : constant FLOAT := - 0.3332985605 ;
    C5 : constant FLOAT := 0.1994653599 ;
    C7 : constant FLOAT := - 0.1390853351 ;
    C9 : constant FLOAT := 0.0964200441 ;
    C11 : constant FLOAT := - 0.0559098861 ;
    C13 : constant FLOAT := 0.0218612288 ;
    C15 : constant FLOAT := - 0.0040540580 ;
    A_2 : FLOAT ;
    Y : FLOAT ;
    A : FLOAT ;
  begin
    A := X ;
    if abs ( A ) > 1.0 then
      A := 1.0 / A ;
    end if ;
    A_2 := A * A ;
    Y := ( C1 +( C3 +( C5 +( C7 +( C9 +( C11 +( C13 + C15 * A_2 ) * A_2) * A_2
       ) * A_2) * A_2) * A_2) * A_2) * A ;
    if abs ( X ) >= 1.0 then
      if X < 0.0 then
        Y := - ( PI_2 + Y ) ;
      else
        Y := PI_2 - Y ;
      end if ;
    end if ;
    return Y ;
  end ATAN ;


  function SQRT ( X : FLOAT ) return FLOAT is -- Copyright Westinghouse 1985

    Y , ROOT_PWR , X_NORM : FLOAT ;
    A : constant FLOAT := 2.1902 ;
    B : constant FLOAT := - 3.0339 ;
    C : constant FLOAT := 1.5451 ;
  begin
    X_NORM := X ;
    ROOT_PWR := 1.0 ;
    if X <= 0.0 then
      return 0.0 ;
    end if ;
    if X > 1.0 then -- REDUCE TO 0.25 .. 1.0
      while X_NORM > 1.0 loop
        ROOT_PWR := ROOT_PWR * 2.0 ;
        X_NORM := X_NORM * 0.25 ;
      end loop ;
    else
      while X_NORM < 0.25 loop
        ROOT_PWR := ROOT_PWR * 0.5 ;
        X_NORM := X_NORM * 4.0 ;
      end loop ;
    end if ;
    Y := A + B / ( C + X_NORM ) ;
    Y := 0.5 * ( Y + X_NORM / Y ) ;
    Y := 0.5 * ( Y + X_NORM / Y ) ;
    Y := Y * ROOT_PWR ;
    return Y ;
  end SQRT ;

  function EXP ( X : FLOAT ) return FLOAT is -- Copyright Westinghouse 1985

    C1 : constant FLOAT := 9.99999900943303E-01 ;
    C2 : constant FLOAT := 5.00006347344554E-01 ;
    C3 : constant FLOAT := 1.66667985598315E-01 ;
    C4 : constant FLOAT := 4.16350120350139E-02 ;
    C5 : constant FLOAT := 8.32859610677671E-03 ;
    C6 : constant FLOAT := 1.43927433449119E-03 ;
    C7 : constant FLOAT := 2.04699933614437E-04 ;

--         4.01169746699903E-07 = MAX_ERROR APPROXIMATION-FUNCTION
    X1 : FLOAT ;
    Y : FLOAT ;
    E_PWR : FLOAT := 1.0 ;
    E : FLOAT := 2.71828182845905 ;
  begin
    X1 := abs ( X ) ;
    if X1 > 88.0 then
      return 0.0 ;
    end if ;
    while X1 >= 1.0 loop
      E_PWR := E_PWR * E * E ;
      X1 := X1 - 2.0 ;
    end loop ;
    Y := 1.0 + ( C1 +( C2 +( C3 +( C4 +( C5 +( C6 + C7 * X1 ) * X1) * X1) * X1
       ) * X1) * X1) * X1 ;
    Y := Y * E_PWR ;
    if X < 0.0 then
      Y := 1.0 / Y ;
    end if ;
    return Y ;
  end EXP ;

  function LOG10 ( X : FLOAT ) return FLOAT is -- Copyright Westinghouse 1985

    C1 : constant FLOAT := 0.868591718 ;
    C3 : constant FLOAT := 0.289335524 ;
    C5 : constant FLOAT := 0.177522071 ;
    C7 : constant FLOAT := 0.094376476 ;
    C9 : constant FLOAT := 0.191337714 ;
    C_R10 : constant FLOAT := 3.1622777 ;
    Y : FLOAT ;
    X_NORM : FLOAT ;
    X_LOG : FLOAT ;
    FRAC : FLOAT ;
    FRAC_2 : FLOAT ;
  begin
    X_LOG := 0.5 ;
    X_NORM := X ;
    if X <= 0.0 then
      return 0.0 ;
    end if ;
    if X >= 10.0 then
      while X_NORM >= 10.0  -- REDUCE TO 1.0 .. 10.0
          loop
        X_LOG := X_LOG + 1.0 ;
        X_NORM := X_NORM * 0.1 ;
      end loop ;
    else
      while X_NORM < 1.0  -- REDUCE TO 1.0 .. 10.0
          loop
        X_LOG := X_LOG - 1.0 ;
        X_NORM := X_NORM * 10.0 ;
      end loop ;
    end if ;
    FRAC := ( X_NORM - C_R10 ) / ( X_NORM + C_R10 ) ;
    FRAC_2 := FRAC * FRAC ;
    Y := ( C1 +( C3 +( C5 +( C7 + C9 * FRAC_2 ) * FRAC_2) * FRAC_2) * FRAC_2)
        * FRAC ;
    return Y + X_LOG ;
  end LOG10 ; -- end of copyrighted section

  function LOG ( X : FLOAT ) return FLOAT is 

  begin
    return 2.302585093 * LOG10 ( X ) ;
  end LOG ;

 

   procedure pa (e : in out vector) is
   -- tests computations with an array as a parameter
      j : integer;
      -- T,T2 : FLOAT are global variables
      begin
        J:=0;
        <<LAB>>
        e (1) := (e (1) + e (2) + e (3) - e (4)) * t;
        e (2) := (e (1) + e (2) - e (3) + e (4)) * t;
        e (3) := (e (1) - e (2) + e (3) + e (4)) * t;
        e (4) := (-e (1) + e (2) + e (3) + e (4)) / t2;
        j := j + 1;
        if J < 6 then
           goto LAB;
        end if;
   end pa;
 
 
   procedure p0 is
   -- tests computations with no parameters
   -- T1,T2 : FLOAT are global
   -- E1 : VECTOR (1..4) is global
   -- J,K,L : INTEGER are global
   begin
      e1 (j) := e1 (k);
      e1 (k) := e1 (l);
      e1 (l) := e1 (j);
   end P0;
 
 
   procedure P3 (X,Y: in out FLOAT;  Z : out FLOAT) is
   -- tests computations with simple identifiers as parameters
   -- T,T2 : FLOAT are global
   begin
      x := t * (x + y);
      y := t * (x + y);
      z := (x + y) / t2;
   end P3;
 
 
begin

   Test ("Whetstone", "Floating point benchmark");

   -- Set constants
   T := 0.499975;
   T1 := 0.50025;
   T2 := 2.0;

   -- Compute the execution frequency for the benchmark modules
   N1 := 0;      --Module 1 not executed
   N2 := 12 * I;
   N3 := 14 * I;
   N4 := 345*I;
   N5 := 0;      -- Module 5 not executed
   N6 := 210*I;
   N7 := 32*I;
   N8 := 899*I;
   N9 := 616*I;
   N10:= 0;      -- Module 10 not executed
   N11:= 93*I;
 
   Start_Time := Clock;  --Get Whetstone start time
 
cycle_loop:
   for cycle_no in 1 .. no_of_cycles loop

      -- Module 1 : computations with simple identifiers
      x1 := 1.0;
      x2 := -1.0;
      x3 := -1.0;
      x4 := -1.0;
      for i in 1 .. n1 loop
         x1 := (x1 + x2 + x3 - x4) * t;
         x2 := (x1 + x2 - x3 + x4) * t;
         x3 := (x1 + x2 + x3 + x4) * t;
         x4 := (-x1 + x2 + x3 + x4) * t;
      end loop;
      -- end Module 1
 
      -- Module 2: computations with array elements
      e1 (1) := 1.0;
      e1 (2) := -1.0;
      e1 (3) := -1.0;
      e1 (4) := -1.0;
      for i in 1 .. n2 loop
         e1 (1) := (e1 (1) + e1 (2) + e1 (3) - e1 (4)) * t;
         e1 (2) := (e1 (1) + e1 (2) - e1 (3) + e1 (4)) * t;
         e1 (3) := (e1 (1) - e1 (2) + e1 (3) + e1 (4)) * t;
         e1 (4) := (-e1 (1) + e1 (2) + e1 (3) + e1 (4)) * t;
      end loop;
      -- end Module 2
 
      -- Module 3 : passing an array as a parmeter
      for i in 1 .. n3 loop
         pa (e1);
      end loop;
      -- end Module 3
 
      -- Module 4 : performing conditional jumps
      j := 1;
      for I in 1..N4 loop
         if J = 1 then
            J := 2;
         else
            J := 3;
         end if;
         if J > 2 then
            J := 0;
         else
            J := 1;
         end if;
         if J < 1 then
            J := 1;
         else
            J := 0;
         end if;
      end loop;
      -- end Module 4
 
      -- Module 5 : omitted
 
      -- Module 6 : performing integer arithmetic
      j := 1;
      k := 2;
      l := 3;
      for i in 1 .. n6 loop
         j := j * (k - j) * (l - k);
         k := l*k - (l - j) * k;
         l := (l - k) * (k + j);
         e1 (l - 1) := float (j + k + l);
         e1 (k - 1) := float (j * k * l);
      end loop;
      -- end Module 6
 
      -- Module 7 : performing computations using trigonometric
      --            functions
      x := 0.5;
      y := 0.5;
      for i in 1 .. n7 loop
         x := t * atan (t2 * sin (x) * cos (x) / (cos (x+y) + cos (x-y)-1.0));
         y := t * atan (t2 * sin (y) * cos (y) / (cos (x+y) + cos (x-y)-1.0));
      end loop;
      -- end Module 7
 
      -- Module 8 : procedure calls with simple identifiers as
      --            parameters
      x := 1.0;
      y := 1.0;
      z := 1.0;
      for i in 1..n8 loop
         p3 (x,y,z);
      end loop;
      -- end Module 8
 
      -- Module 9 : array reference and procedure calls with no
      --            parameters
      j := 1;
      k := 2;
      l := 3;
      e1 (1) := 1.0;
      e1 (2) := 2.0;
      e1 (3) := 3.0;
      for i in 1 .. n9 loop
         p0;
      end loop;
      -- end Module 9
 
      -- Module 10 : integer arithmetic
      j := 2;
      k := 3;
      for i in 1 .. n10 loop
         j := j + k;
         k := k + j;
         j := k - j;
         k := k - j - j;
      end loop;
      -- end Module 10
 
      -- Module 11 : performing computations using standard
      --            mathematical functions
      x := 0.75;
      for i in 1 .. n11 loop
         x := sqrt (exp (log (x) / t1));
      end loop;
      -- end Moudle 11
 
   end loop cycle_loop;
 
   Stop_Time := Clock;

   Put ("Time taken = ");
   Put (Integer (1000.0 * To_Duration (Time_Span'(Stop_Time - Start_Time))));
   Put (" mSec");
   New_Line;

   Rating := Integer (1000.0 * Float (no_of_cycles) / 
      Float (To_Duration (Time_Span'(Stop_Time - Start_Time))));

   Put ("Whetstone rating = ");
   Put (Rating, 0);
   Put (" KWIPS");
   New_Line;

   Result;

end whetstone;