old-autovect-branch/gcc/testsuite/ada/acats/tests/cxg/cxg2009.a

   1 -- CXG2009.A
   2 --
   3 --                             Grant of Unlimited Rights
   4 --
   5 --     Under contracts F33600-87-D-0337, F33600-84-D-0280, MDA903-79-C-0687,
   6 --     F08630-91-C-0015, and DCA100-97-D-0025, the U.S. Government obtained
   7 --     unlimited rights in the software and documentation contained herein.
   8 --     Unlimited rights are defined in DFAR 252.227-7013(a)(19).  By making
   9 --     this public release, the Government intends to confer upon all
  10 --     recipients unlimited rights  equal to those held by the Government.
  11 --     These rights include rights to use, duplicate, release or disclose the
  12 --     released technical data and computer software in whole or in part, in
  13 --     any manner and for any purpose whatsoever, and to have or permit others
  14 --     to do so.
  15 --
  16 --                                    DISCLAIMER
  17 --
  18 --     ALL MATERIALS OR INFORMATION HEREIN RELEASED, MADE AVAILABLE OR
  19 --     DISCLOSED ARE AS IS.  THE GOVERNMENT MAKES NO EXPRESS OR IMPLIED
  20 --     WARRANTY AS TO ANY MATTER WHATSOEVER, INCLUDING THE CONDITIONS OF THE
  21 --     SOFTWARE, DOCUMENTATION OR OTHER INFORMATION RELEASED, MADE AVAILABLE
  22 --     OR DISCLOSED, OR THE OWNERSHIP, MERCHANTABILITY, OR FITNESS FOR A
  23 --     PARTICULAR PURPOSE OF SAID MATERIAL.
  24 --*
  25 --
  26 -- OBJECTIVE:
  27 --      Check that the real sqrt and complex modulus functions
  28 --      return results that are within the allowed
  29 --      error bound.
  30 --
  31 -- TEST DESCRIPTION:
  32 --      This test checks the accuracy of the sqrt and modulus functions
  33 --      by computing the norm of various vectors where the result
  34 --      is known in advance.
  35 --      This test uses real and complex math together as would an
  36 --      actual application.  Considerable use of generics is also
  37 --      employed.
  38 --
  39 -- SPECIAL REQUIREMENTS
  40 --      The Strict Mode for the numerical accuracy must be
  41 --      selected.  The method by which this mode is selected
  42 --      is implementation dependent.
  43 --
  44 -- APPLICABILITY CRITERIA:
  45 --      This test applies only to implementations supporting the
  46 --      Numerics Annex.
  47 --      This test only applies to the Strict Mode for numerical
  48 --      accuracy.
  49 --
  50 --
  51 -- CHANGE HISTORY:
  52 --      26 FEB 96   SAIC    Initial release for 2.1
  53 --      22 AUG 96   SAIC    Revised Check procedure
  54 --
  55 --!
  56
  57 ------------------------------------------------------------------------------
  58
  59 with System;
  60 with Report;
  61 with Ada.Numerics.Generic_Complex_Types;
  62 with Ada.Numerics.Generic_Elementary_Functions;
  63 procedure CXG2009 is
  64    Verbose : constant Boolean := False;
  65
  66    --=====================================================================
  67
  68    generic
  69       type Real is digits <>;
  70    package Generic_Real_Norm_Check is
  71       procedure Do_Test;
  72    end Generic_Real_Norm_Check;
  73
  74    -----------------------------------------------------------------------
  75
  76    package body Generic_Real_Norm_Check is
  77       type Vector is array (Integer range <>) of Real;
  78
  79       package GEF is new Ada.Numerics.Generic_Elementary_Functions (Real);
  80       function Sqrt (X : Real) return Real renames GEF.Sqrt;
  81
  82       function One_Norm (V : Vector) return Real is
  83       -- sum of absolute values of the elements of the vector
  84          Result : Real := 0.0;
  85       begin
  86          for I in V'Range loop
  87             Result := Result + abs V(I);
  88          end loop;
  89          return Result;
  90       end One_Norm;
  91
  92       function Inf_Norm (V : Vector) return Real is
  93       -- greatest absolute vector element
  94          Result : Real := 0.0;
  95       begin
  96          for I in V'Range loop
  97             if abs V(I) > Result then
  98                Result := abs V(I);
  99             end if;
 100          end loop;
 101          return Result;
 102       end Inf_Norm;
 103
 104       function Two_Norm (V : Vector) return Real is
 105       -- if greatest absolute vector element is 0 then return 0
 106       -- else return greatest * sqrt (sum((element / greatest) ** 2)))
 107       --   where greatest is Inf_Norm of the vector
 108          Inf_N : Real;
 109          Sum_Squares : Real;
 110          Term : Real;
 111       begin
 112          Inf_N := Inf_Norm (V);
 113          if Inf_N = 0.0 then
 114             return 0.0;
 115          end if;
 116          Sum_Squares := 0.0;
 117          for I in V'Range loop
 118             Term := V (I) / Inf_N;
 119             Sum_Squares := Sum_Squares + Term * Term;
 120          end loop;
 121          return Inf_N * Sqrt (Sum_Squares);
 122       end Two_Norm;
 123
 124
 125       procedure Check (Actual, Expected : Real;
 126                        Test_Name : String;
 127                        MRE : Real;
 128                        Vector_Length : Integer) is
 129          Rel_Error : Real;
 130          Abs_Error : Real;
 131          Max_Error : Real;
 132       begin
 133          -- In the case where the expected result is very small or 0
 134          -- we compute the maximum error as a multiple of Model_Epsilon instead
 135          -- of Model_Epsilon and Expected.
 136          Rel_Error := MRE * abs Expected * Real'Model_Epsilon;
 137          Abs_Error := MRE * Real'Model_Epsilon;
 138          if Rel_Error > Abs_Error then
 139             Max_Error := Rel_Error;
 140          else
 141             Max_Error := Abs_Error;
 142          end if;
 143
 144          if abs (Actual - Expected) > Max_Error then
 145             Report.Failed (Test_Name &
 146                              "  VectLength:" &
 147                            Integer'Image (Vector_Length) &
 148                            " actual: " & Real'Image (Actual) &
 149                            " expected: " & Real'Image (Expected) &
 150                            " difference: " &
 151                            Real'Image (Actual - Expected) &
 152                            " mre:" & Real'Image (Max_Error) );
 153          elsif Verbose then
 154             Report.Comment (Test_Name & " vector length" &
 155                             Integer'Image (Vector_Length));
 156           end if;
 157       end Check;
 158
 159
 160       procedure Do_Test is
 161       begin
 162          for Vector_Length in 1 .. 10 loop
 163             declare
 164                V  : Vector (1..Vector_Length) := (1..Vector_Length => 0.0);
 165                V1 : Vector (1..Vector_Length) := (1..Vector_Length => 1.0);
 166             begin
 167                Check (One_Norm (V), 0.0, "one_norm (z)", 0.0, Vector_Length);
 168                Check (Inf_Norm (V), 0.0, "inf_norm (z)", 0.0, Vector_Length);
 169
 170                for J in 1..Vector_Length loop
 171                  V := (1..Vector_Length => 0.0);
 172                  V (J) := 1.0;
 173                  Check (One_Norm (V), 1.0, "one_norm (010)",
 174                         0.0, Vector_Length);
 175                  Check (Inf_Norm (V), 1.0, "inf_norm (010)",
 176                         0.0, Vector_Length);
 177                  Check (Two_Norm (V), 1.0, "two_norm (010)",
 178                         0.0, Vector_Length);
 179                end loop;
 180
 181                Check (One_Norm (V1), Real (Vector_Length), "one_norm (1)",
 182                       0.0, Vector_Length);
 183                Check (Inf_Norm (V1), 1.0, "inf_norm (1)",
 184                       0.0, Vector_Length);
 185
 186                -- error in computing Two_Norm and expected result
 187                -- are as follows  (ME is Model_Epsilon * Expected_Value):
 188                --   2ME from expected Sqrt
 189                --   2ME from Sqrt in Two_Norm times the error in the
 190                --   vector calculation.
 191                --   The vector calculation contains the following error
 192                --   based upon the length N of the vector:
 193                --      N*1ME from squaring terms in Two_Norm
 194                --      N*1ME from the division of each term in Two_Norm
 195                --      (N-1)*1ME from the sum of the terms
 196                -- This gives (2 + 2 * (N + N + (N-1)) ) * ME
 197                -- which simplifies to (2 + 2N + 2N + 2N - 2) * ME
 198                -- or 6*N*ME
 199                Check (Two_Norm (V1), Sqrt (Real(Vector_Length)),
 200                       "two_norm (1)",
 201                       (Real (6 * Vector_Length)),
 202                       Vector_Length);
 203             exception
 204                when others => Report.Failed ("exception for vector length" &
 205                                 Integer'Image (Vector_Length) );
 206             end;
 207          end loop;
 208       end Do_Test;
 209    end Generic_Real_Norm_Check;
 210
 211    --=====================================================================
 212
 213    generic
 214       type Real is digits <>;
 215    package Generic_Complex_Norm_Check is
 216       procedure Do_Test;
 217    end Generic_Complex_Norm_Check;
 218
 219    -----------------------------------------------------------------------
 220
 221    package body Generic_Complex_Norm_Check is
 222       package Complex_Types is new Ada.Numerics.Generic_Complex_Types (Real);
 223       use Complex_Types;
 224       type Vector is array (Integer range <>) of Complex;
 225
 226       package GEF is new Ada.Numerics.Generic_Elementary_Functions (Real);
 227       function Sqrt (X : Real) return Real renames GEF.Sqrt;
 228
 229       function One_Norm (V : Vector) return Real is
 230          Result : Real := 0.0;
 231       begin
 232          for I in V'Range loop
 233             Result := Result + abs V(I);
 234          end loop;
 235          return Result;
 236       end One_Norm;
 237
 238       function Inf_Norm (V : Vector) return Real is
 239          Result : Real := 0.0;
 240       begin
 241          for I in V'Range loop
 242             if abs V(I) > Result then
 243                Result := abs V(I);
 244             end if;
 245          end loop;
 246          return Result;
 247       end Inf_Norm;
 248
 249       function Two_Norm (V : Vector) return Real is
 250          Inf_N : Real;
 251          Sum_Squares : Real;
 252          Term : Real;
 253       begin
 254          Inf_N := Inf_Norm (V);
 255          if Inf_N = 0.0 then
 256             return 0.0;
 257          end if;
 258          Sum_Squares := 0.0;
 259          for I in V'Range loop
 260             Term := abs (V (I) / Inf_N );
 261             Sum_Squares := Sum_Squares + Term * Term;
 262          end loop;
 263          return Inf_N * Sqrt (Sum_Squares);
 264       end Two_Norm;
 265
 266
 267       procedure Check (Actual, Expected : Real;
 268                        Test_Name : String;
 269                        MRE : Real;
 270                        Vector_Length : Integer) is
 271          Rel_Error : Real;
 272          Abs_Error : Real;
 273          Max_Error : Real;
 274       begin
 275          -- In the case where the expected result is very small or 0
 276          -- we compute the maximum error as a multiple of Model_Epsilon instead
 277          -- of Model_Epsilon and Expected.
 278          Rel_Error := MRE * abs Expected * Real'Model_Epsilon;
 279          Abs_Error := MRE * Real'Model_Epsilon;
 280          if Rel_Error > Abs_Error then
 281             Max_Error := Rel_Error;
 282          else
 283             Max_Error := Abs_Error;
 284          end if;
 285
 286          if abs (Actual - Expected) > Max_Error then
 287             Report.Failed (Test_Name &
 288                              "  VectLength:" &
 289                            Integer'Image (Vector_Length) &
 290                            " actual: " & Real'Image (Actual) &
 291                            " expected: " & Real'Image (Expected) &
 292                            " difference: " &
 293                            Real'Image (Actual - Expected) &
 294                            " mre:" & Real'Image (Max_Error) );
 295          elsif Verbose then
 296             Report.Comment (Test_Name & " vector length" &
 297                             Integer'Image (Vector_Length));
 298           end if;
 299       end Check;
 300
 301
 302       procedure Do_Test is
 303       begin
 304          for Vector_Length in 1 .. 10 loop
 305             declare
 306                V  : Vector (1..Vector_Length) :=
 307                       (1..Vector_Length => (0.0, 0.0));
 308                X, Y : Vector (1..Vector_Length);
 309             begin
 310                Check (One_Norm (V), 0.0, "one_norm (z)", 0.0, Vector_Length);
 311                Check (Inf_Norm (V), 0.0, "inf_norm (z)", 0.0, Vector_Length);
 312
 313                for J in 1..Vector_Length loop
 314                  X := (1..Vector_Length => (0.0, 0.0) );
 315                  Y := X;   -- X and Y are now both zeroed
 316                  X (J).Re := 1.0;
 317                  Y (J).Im := 1.0;
 318                  Check (One_Norm (X), 1.0, "one_norm (0x0)",
 319                         0.0, Vector_Length);
 320                  Check (Inf_Norm (X), 1.0, "inf_norm (0x0)",
 321                         0.0, Vector_Length);
 322                  Check (Two_Norm (X), 1.0, "two_norm (0x0)",
 323                         0.0, Vector_Length);
 324                  Check (One_Norm (Y), 1.0, "one_norm (0y0)",
 325                         0.0, Vector_Length);
 326                  Check (Inf_Norm (Y), 1.0, "inf_norm (0y0)",
 327                         0.0, Vector_Length);
 328                  Check (Two_Norm (Y), 1.0, "two_norm (0y0)",
 329                         0.0, Vector_Length);
 330                end loop;
 331
 332                V := (1..Vector_Length => (3.0, 4.0));
 333
 334                -- error in One_Norm is 3*N*ME for abs computation +
 335                --  (N-1)*ME for the additions
 336                -- which gives (4N-1) * ME
 337                Check (One_Norm (V), 5.0 * Real (Vector_Length),
 338                       "one_norm ((3,4))",
 339                       Real (4*Vector_Length - 1),
 340                       Vector_Length);
 341
 342                -- error in Inf_Norm is from abs of single element (3ME)
 343                Check (Inf_Norm (V), 5.0,
 344                       "inf_norm ((3,4))",
 345                       3.0,
 346                       Vector_Length);
 347
 348                -- error in following comes from:
 349                --   2ME in sqrt of expected result
 350                --   3ME in Inf_Norm calculation
 351                --   2ME in sqrt of vector calculation
 352                --   vector calculation has following error
 353                --      3N*ME for abs
 354                --       N*ME for squaring
 355                --       N*ME for division
 356                --       (N-1)ME for sum
 357                -- this results in [2 + 3 + 2(6N-1) ] * ME
 358                -- or (12N + 3)ME
 359                Check (Two_Norm (V), 5.0 * Sqrt (Real(Vector_Length)),
 360                       "two_norm ((3,4))",
 361                       (12.0 * Real (Vector_Length) + 3.0),
 362                       Vector_Length);
 363             exception
 364                when others => Report.Failed ("exception for complex " &
 365                                              "vector length" &
 366                                              Integer'Image (Vector_Length) );
 367             end;
 368          end loop;
 369       end Do_Test;
 370    end Generic_Complex_Norm_Check;
 371
 372    --=====================================================================
 373
 374    generic
 375       type Real is digits <>;
 376    package Generic_Norm_Check is
 377       procedure Do_Test;
 378    end Generic_Norm_Check;
 379
 380    -----------------------------------------------------------------------
 381
 382    package body Generic_Norm_Check is
 383       package RNC is new Generic_Real_Norm_Check (Real);
 384       package CNC is new Generic_Complex_Norm_Check (Real);
 385       procedure Do_Test is
 386       begin
 387          RNC.Do_Test;
 388          CNC.Do_Test;
 389       end Do_Test;
 390    end Generic_Norm_Check;
 391
 392    --=====================================================================
 393
 394    package Float_Check is new Generic_Norm_Check (Float);
 395
 396    type A_Long_Float is digits System.Max_Digits;
 397    package A_Long_Float_Check is new Generic_Norm_Check (A_Long_Float);
 398
 399    -----------------------------------------------------------------------
 400
 401 begin
 402    Report.Test ("CXG2009",
 403                 "Check the accuracy of the real sqrt and complex " &
 404                 " modulus functions");
 405
 406    if Verbose then
 407       Report.Comment ("checking Standard.Float");
 408    end if;
 409
 410    Float_Check.Do_Test;
 411
 412    if Verbose then
 413       Report.Comment ("checking a digits" &
 414                       Integer'Image (System.Max_Digits) &
 415                       " floating point type");
 416    end if;
 417
 418    A_Long_Float_Check.Do_Test;
 419
 420    Report.Result;
 421 end CXG2009;