-- CXG2009.A
--
--                             Grant of Unlimited Rights
--
--     Under contracts F33600-87-D-0337, F33600-84-D-0280, MDA903-79-C-0687,
--     F08630-91-C-0015, and DCA100-97-D-0025, the U.S. Government obtained 
--     unlimited rights in the software and documentation contained herein.
--     Unlimited rights are defined in DFAR 252.227-7013(a)(19).  By making 
--     this public release, the Government intends to confer upon all 
--     recipients unlimited rights  equal to those held by the Government.  
--     These rights include rights to use, duplicate, release or disclose the 
--     released technical data and computer software in whole or in part, in 
--     any manner and for any purpose whatsoever, and to have or permit others 
--     to do so.
--
--                                    DISCLAIMER
--
--     ALL MATERIALS OR INFORMATION HEREIN RELEASED, MADE AVAILABLE OR
--     DISCLOSED ARE AS IS.  THE GOVERNMENT MAKES NO EXPRESS OR IMPLIED 
--     WARRANTY AS TO ANY MATTER WHATSOEVER, INCLUDING THE CONDITIONS OF THE
--     SOFTWARE, DOCUMENTATION OR OTHER INFORMATION RELEASED, MADE AVAILABLE 
--     OR DISCLOSED, OR THE OWNERSHIP, MERCHANTABILITY, OR FITNESS FOR A
--     PARTICULAR PURPOSE OF SAID MATERIAL.
--*
--
-- OBJECTIVE:
--      Check that the real sqrt and complex modulus functions
--      return results that are within the allowed
--      error bound.
--
-- TEST DESCRIPTION:
--      This test checks the accuracy of the sqrt and modulus functions 
--      by computing the norm of various vectors where the result
--      is known in advance.
--      This test uses real and complex math together as would an 
--      actual application.  Considerable use of generics is also
--      employed.
--
-- SPECIAL REQUIREMENTS
--      The Strict Mode for the numerical accuracy must be
--      selected.  The method by which this mode is selected
--      is implementation dependent.
--
-- APPLICABILITY CRITERIA:
--      This test applies only to implementations supporting the
--      Numerics Annex.
--      This test only applies to the Strict Mode for numerical
--      accuracy.
--
--
-- CHANGE HISTORY:
--      26 FEB 96   SAIC    Initial release for 2.1
--      22 AUG 96   SAIC    Revised Check procedure
--
--!

------------------------------------------------------------------------------

with System;
with Report;
with Ada.Numerics.Generic_Complex_Types;
with Ada.Numerics.Generic_Elementary_Functions;
procedure CXG2009 is
   Verbose : constant Boolean := False;

   --=====================================================================

   generic
      type Real is digits <>;
   package Generic_Real_Norm_Check is
      procedure Do_Test;
   end Generic_Real_Norm_Check;

   -----------------------------------------------------------------------

   package body Generic_Real_Norm_Check is
      type Vector is array (Integer range <>) of Real;

      package GEF is new Ada.Numerics.Generic_Elementary_Functions (Real);
      function Sqrt (X : Real) return Real renames GEF.Sqrt;

      function One_Norm (V : Vector) return Real is
      -- sum of absolute values of the elements of the vector
	 Result : Real := 0.0;
      begin
	 for I in V'Range loop
	    Result := Result + abs V(I);
	 end loop;
	 return Result;
      end One_Norm;

      function Inf_Norm (V : Vector) return Real is
      -- greatest absolute vector element
	 Result : Real := 0.0;
      begin
	 for I in V'Range loop
	    if abs V(I) > Result then
	       Result := abs V(I);
	    end if;
	 end loop;
	 return Result;
      end Inf_Norm;

      function Two_Norm (V : Vector) return Real is
      -- if greatest absolute vector element is 0 then return 0
      -- else return greatest * sqrt (sum((element / greatest) ** 2)))
      --   where greatest is Inf_Norm of the vector
	 Inf_N : Real;
	 Sum_Squares : Real;
	 Term : Real;
      begin
	 Inf_N := Inf_Norm (V);
	 if Inf_N = 0.0 then
	    return 0.0;
	 end if;
         Sum_Squares := 0.0;
	 for I in V'Range loop
	    Term := V (I) / Inf_N;
	    Sum_Squares := Sum_Squares + Term * Term;
	 end loop;
	 return Inf_N * Sqrt (Sum_Squares);
      end Two_Norm;


      procedure Check (Actual, Expected : Real;
		       Test_Name : String;
		       MRE : Real;
		       Vector_Length : Integer) is
         Rel_Error : Real;
         Abs_Error : Real;
         Max_Error : Real;
      begin
         -- In the case where the expected result is very small or 0
         -- we compute the maximum error as a multiple of Model_Epsilon instead
         -- of Model_Epsilon and Expected.
         Rel_Error := MRE * abs Expected * Real'Model_Epsilon;
         Abs_Error := MRE * Real'Model_Epsilon;
         if Rel_Error > Abs_Error then
            Max_Error := Rel_Error;
         else
            Max_Error := Abs_Error;
         end if; 

         if abs (Actual - Expected) > Max_Error then
            Report.Failed (Test_Name & 
	                     "  VectLength:" & 
                           Integer'Image (Vector_Length) &
                           " actual: " & Real'Image (Actual) &
                           " expected: " & Real'Image (Expected) &
                           " difference: " & 
                           Real'Image (Actual - Expected) &
                           " mre:" & Real'Image (Max_Error) );
         elsif Verbose then
            Report.Comment (Test_Name & " vector length" &
                            Integer'Image (Vector_Length));
	  end if;
      end Check;


      procedure Do_Test is
      begin
	 for Vector_Length in 1 .. 10 loop
	    declare
	       V  : Vector (1..Vector_Length) := (1..Vector_Length => 0.0);
	       V1 : Vector (1..Vector_Length) := (1..Vector_Length => 1.0);
	    begin
	       Check (One_Norm (V), 0.0, "one_norm (z)", 0.0, Vector_Length);
	       Check (Inf_Norm (V), 0.0, "inf_norm (z)", 0.0, Vector_Length);

	       for J in 1..Vector_Length loop
		 V := (1..Vector_Length => 0.0);
		 V (J) := 1.0;
	         Check (One_Norm (V), 1.0, "one_norm (010)", 
			0.0, Vector_Length);
	         Check (Inf_Norm (V), 1.0, "inf_norm (010)", 
			0.0, Vector_Length);
	         Check (Two_Norm (V), 1.0, "two_norm (010)", 
			0.0, Vector_Length);
	       end loop;

	       Check (One_Norm (V1), Real (Vector_Length), "one_norm (1)", 
		      0.0, Vector_Length);
	       Check (Inf_Norm (V1), 1.0, "inf_norm (1)", 
		      0.0, Vector_Length);

               -- error in computing Two_Norm and expected result
               -- are as follows  (ME is Model_Epsilon * Expected_Value):
               --   2ME from expected Sqrt
               --   2ME from Sqrt in Two_Norm times the error in the
               --   vector calculation.
               --   The vector calculation contains the following error
               --   based upon the length N of the vector:
               --      N*1ME from squaring terms in Two_Norm
               --      N*1ME from the division of each term in Two_Norm
               --      (N-1)*1ME from the sum of the terms
               -- This gives (2 + 2 * (N + N + (N-1)) ) * ME
               -- which simplifies to (2 + 2N + 2N + 2N - 2) * ME
               -- or 6*N*ME
	       Check (Two_Norm (V1), Sqrt (Real(Vector_Length)), 
                      "two_norm (1)", 
		      (Real (6 * Vector_Length)), 
		      Vector_Length);
	    exception
	       when others => Report.Failed ("exception for vector length" &
				Integer'Image (Vector_Length) );
	    end;
	 end loop;
      end Do_Test;
   end Generic_Real_Norm_Check;

   --=====================================================================

   generic
      type Real is digits <>;
   package Generic_Complex_Norm_Check is
      procedure Do_Test;
   end Generic_Complex_Norm_Check;

   -----------------------------------------------------------------------

   package body Generic_Complex_Norm_Check is
      package Complex_Types is new Ada.Numerics.Generic_Complex_Types (Real);
      use Complex_Types;
      type Vector is array (Integer range <>) of Complex;

      package GEF is new Ada.Numerics.Generic_Elementary_Functions (Real);
      function Sqrt (X : Real) return Real renames GEF.Sqrt;

      function One_Norm (V : Vector) return Real is
	 Result : Real := 0.0;
      begin
	 for I in V'Range loop
	    Result := Result + abs V(I);
	 end loop;
	 return Result;
      end One_Norm;

      function Inf_Norm (V : Vector) return Real is
	 Result : Real := 0.0;
      begin
	 for I in V'Range loop
	    if abs V(I) > Result then
	       Result := abs V(I);
	    end if;
	 end loop;
	 return Result;
      end Inf_Norm;

      function Two_Norm (V : Vector) return Real is
	 Inf_N : Real;
	 Sum_Squares : Real;
	 Term : Real;
      begin
	 Inf_N := Inf_Norm (V);
	 if Inf_N = 0.0 then
	    return 0.0;
	 end if;
         Sum_Squares := 0.0;
	 for I in V'Range loop
	    Term := abs (V (I) / Inf_N );
	    Sum_Squares := Sum_Squares + Term * Term;
	 end loop;
	 return Inf_N * Sqrt (Sum_Squares);
      end Two_Norm;


      procedure Check (Actual, Expected : Real;
		       Test_Name : String;
		       MRE : Real;
		       Vector_Length : Integer) is
         Rel_Error : Real;
         Abs_Error : Real;
         Max_Error : Real;
      begin
         -- In the case where the expected result is very small or 0
         -- we compute the maximum error as a multiple of Model_Epsilon instead
         -- of Model_Epsilon and Expected.
         Rel_Error := MRE * abs Expected * Real'Model_Epsilon;
         Abs_Error := MRE * Real'Model_Epsilon;
         if Rel_Error > Abs_Error then
            Max_Error := Rel_Error;
         else
            Max_Error := Abs_Error;
         end if; 

         if abs (Actual - Expected) > Max_Error then
            Report.Failed (Test_Name & 
	                     "  VectLength:" & 
                           Integer'Image (Vector_Length) &
                           " actual: " & Real'Image (Actual) &
                           " expected: " & Real'Image (Expected) &
                           " difference: " & 
                           Real'Image (Actual - Expected) &
                           " mre:" & Real'Image (Max_Error) );
         elsif Verbose then
            Report.Comment (Test_Name & " vector length" &
                            Integer'Image (Vector_Length));
	  end if;
      end Check;


      procedure Do_Test is
      begin
	 for Vector_Length in 1 .. 10 loop
	    declare
	       V  : Vector (1..Vector_Length) := 
                      (1..Vector_Length => (0.0, 0.0));
               X, Y : Vector (1..Vector_Length);
	    begin
	       Check (One_Norm (V), 0.0, "one_norm (z)", 0.0, Vector_Length);
	       Check (Inf_Norm (V), 0.0, "inf_norm (z)", 0.0, Vector_Length);

	       for J in 1..Vector_Length loop
		 X := (1..Vector_Length => (0.0, 0.0) );
                 Y := X;   -- X and Y are now both zeroed
		 X (J).Re := 1.0;
                 Y (J).Im := 1.0;
	         Check (One_Norm (X), 1.0, "one_norm (0x0)", 
			0.0, Vector_Length);
	         Check (Inf_Norm (X), 1.0, "inf_norm (0x0)", 
			0.0, Vector_Length);
	         Check (Two_Norm (X), 1.0, "two_norm (0x0)", 
			0.0, Vector_Length);
	         Check (One_Norm (Y), 1.0, "one_norm (0y0)", 
			0.0, Vector_Length);
	         Check (Inf_Norm (Y), 1.0, "inf_norm (0y0)", 
			0.0, Vector_Length);
	         Check (Two_Norm (Y), 1.0, "two_norm (0y0)", 
			0.0, Vector_Length);
	       end loop;

               V := (1..Vector_Length => (3.0, 4.0));

               -- error in One_Norm is 3*N*ME for abs computation +
               --  (N-1)*ME for the additions
               -- which gives (4N-1) * ME
	       Check (One_Norm (V), 5.0 * Real (Vector_Length), 
		      "one_norm ((3,4))", 
		      Real (4*Vector_Length - 1), 
		      Vector_Length);

               -- error in Inf_Norm is from abs of single element (3ME)
	       Check (Inf_Norm (V), 5.0, 
		      "inf_norm ((3,4))", 
		      3.0, 
		      Vector_Length);

               -- error in following comes from:
               --   2ME in sqrt of expected result
               --   3ME in Inf_Norm calculation
               --   2ME in sqrt of vector calculation
               --   vector calculation has following error
               --      3N*ME for abs
               --       N*ME for squaring
               --       N*ME for division
               --       (N-1)ME for sum
               -- this results in [2 + 3 + 2(6N-1) ] * ME
               -- or (12N + 3)ME
	       Check (Two_Norm (V), 5.0 * Sqrt (Real(Vector_Length)), 
                      "two_norm ((3,4))", 
		      (12.0 * Real (Vector_Length) + 3.0), 
		      Vector_Length);
	    exception
	       when others => Report.Failed ("exception for complex " &
                                             "vector length" &
	                                     Integer'Image (Vector_Length) );
	    end;
	 end loop;
      end Do_Test;
   end Generic_Complex_Norm_Check;

   --=====================================================================

   generic
      type Real is digits <>;
   package Generic_Norm_Check is
      procedure Do_Test;
   end Generic_Norm_Check;

   -----------------------------------------------------------------------

   package body Generic_Norm_Check is
      package RNC is new Generic_Real_Norm_Check (Real);
      package CNC is new Generic_Complex_Norm_Check (Real);
      procedure Do_Test is
      begin
         RNC.Do_Test;
         CNC.Do_Test;
      end Do_Test;
   end Generic_Norm_Check;

   --=====================================================================

   package Float_Check is new Generic_Norm_Check (Float);

   type A_Long_Float is digits System.Max_Digits;
   package A_Long_Float_Check is new Generic_Norm_Check (A_Long_Float);

   -----------------------------------------------------------------------

begin
   Report.Test ("CXG2009",
                "Check the accuracy of the real sqrt and complex " &
                " modulus functions"); 

   if Verbose then
      Report.Comment ("checking Standard.Float");
   end if;

   Float_Check.Do_Test;

   if Verbose then
      Report.Comment ("checking a digits" & 
                      Integer'Image (System.Max_Digits) &
                      " floating point type");
   end if;

   A_Long_Float_Check.Do_Test;

   Report.Result;
end CXG2009;