Subversion Repositories openrisc

-- 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;