gcc/testsuite/ada/acats/tests/cxg/cxg2018.a

   1 -- CXG2018.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 complex EXP function returns
  28 --      a result that is within the error bound allowed.
  29 --
  30 -- TEST DESCRIPTION:
  31 --      This test consists of a generic package that is
  32 --      instantiated to check complex numbers based upon
  33 --      both Float and a long float type.
  34 --      The test for each floating point type is divided into
  35 --      several parts:
  36 --         Special value checks where the result is a known constant.
  37 --         Checks that use an identity for determining the result.
  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 --      21 Mar 96   SAIC    Initial release for 2.1
  53 --      17 Aug 96   SAIC    Incorporated reviewer comments.
  54 --      27 Aug 99   RLB     Repair on the error result of checks.
  55 --      02 Apr 03   RLB     Added code to discard excess precision in the
  56 --                          construction of the test value for the
  57 --                          Identity_Test.
  58 --
  59 --!
  60
  61 --
  62 -- References:
  63 --
  64 -- W. J. Cody
  65 -- CELEFUNT: A Portable Test Package for Complex Elementary Functions
  66 -- Algorithm 714, Collected Algorithms from ACM.
  67 -- Published in Transactions On Mathematical Software,
  68 -- Vol. 19, No. 1, March, 1993, pp. 1-21.
  69 --
  70 -- CRC Standard Mathematical Tables
  71 -- 23rd Edition
  72 --
  73
  74 with System;
  75 with Report;
  76 with Ada.Numerics.Generic_Complex_Types;
  77 with Ada.Numerics.Generic_Complex_Elementary_Functions;
  78 procedure CXG2018 is
  79    Verbose : constant Boolean := False;
  80    -- Note that Max_Samples is the number of samples taken in
  81    -- both the real and imaginary directions.  Thus, for Max_Samples
  82    -- of 100 the number of values checked is 10000.
  83    Max_Samples : constant := 100;
  84
  85    E  : constant := Ada.Numerics.E;
  86    Pi : constant := Ada.Numerics.Pi;
  87
  88    generic
  89       type Real is digits <>;
  90    package Generic_Check is
  91       procedure Do_Test;
  92    end Generic_Check;
  93
  94    package body Generic_Check is
  95       package Complex_Type is new
  96            Ada.Numerics.Generic_Complex_Types (Real);
  97       use Complex_Type;
  98
  99       package CEF is new
 100            Ada.Numerics.Generic_Complex_Elementary_Functions (Complex_Type);
 101
 102       function Exp (X : Complex) return Complex renames CEF.Exp;
 103       function Exp (X : Imaginary) return Complex renames CEF.Exp;
 104
 105       -- flag used to terminate some tests early
 106       Accuracy_Error_Reported : Boolean := False;
 107
 108
 109       -- The following value is a lower bound on the accuracy
 110       -- required.  It is normally 0.0 so that the lower bound
 111       -- is computed from Model_Epsilon.  However, for tests
 112       -- where the expected result is only known to a certain
 113       -- amount of precision this bound takes on a non-zero
 114       -- value to account for that level of precision.
 115       Error_Low_Bound : Real := 0.0;
 116
 117       procedure Check (Actual, Expected : Real;
 118                        Test_Name : String;
 119                        MRE : Real) is
 120          Max_Error : Real;
 121          Rel_Error : Real;
 122          Abs_Error : Real;
 123       begin
 124          -- In the case where the expected result is very small or 0
 125          -- we compute the maximum error as a multiple of Model_Small instead
 126          -- of Model_Epsilon and Expected.
 127          Rel_Error := MRE * abs Expected * Real'Model_Epsilon;
 128          Abs_Error := MRE * Real'Model_Small;
 129          if Rel_Error > Abs_Error then
 130             Max_Error := Rel_Error;
 131          else
 132             Max_Error := Abs_Error;
 133          end if;
 134
 135          -- take into account the low bound on the error
 136          if Max_Error < Error_Low_Bound then
 137             Max_Error := Error_Low_Bound;
 138          end if;
 139
 140          if abs (Actual - Expected) > Max_Error then
 141             Accuracy_Error_Reported := True;
 142             Report.Failed (Test_Name &
 143                            " actual: " & Real'Image (Actual) &
 144                            " expected: " & Real'Image (Expected) &
 145                            " difference: " & Real'Image (Actual - Expected) &
 146                            " max err:" & Real'Image (Max_Error) );
 147          elsif Verbose then
 148             if Actual = Expected then
 149                Report.Comment (Test_Name & "  exact result");
 150             else
 151                Report.Comment (Test_Name & "  passed");
 152             end if;
 153          end if;
 154       end Check;
 155
 156
 157       procedure Check (Actual, Expected : Complex;
 158                        Test_Name : String;
 159                        MRE : Real) is
 160       begin
 161          Check (Actual.Re, Expected.Re, Test_Name & " real part", MRE);
 162          Check (Actual.Im, Expected.Im, Test_Name & " imaginary part", MRE);
 163       end Check;
 164
 165
 166       procedure Special_Value_Test is
 167          -- In the following tests the expected result is accurate
 168          -- to the machine precision so the minimum guaranteed error
 169          -- bound can be used.
 170          --
 171          -- The error bounds given assumed z is  exact.  When using
 172          -- pi there is an extra error of 1.0ME.
 173          -- The pi inside the exp call requires that the complex
 174          -- component have an extra error allowance of 1.0*angle*ME.
 175          -- Thus for pi/2,the Minimum_Error_I is
 176          -- (2.0 + 1.0(pi/2))ME <= 3.6ME.
 177          -- For pi, it is (2.0 + 1.0*pi)ME <= 5.2ME,
 178          -- and for 2pi, it is (2.0 + 1.0(2pi))ME <= 8.3ME.
 179
 180          -- The addition of 1 or i to a result is so that neither of
 181          -- the components of an expected result is 0.  This is so
 182          -- that a reasonable relative error is allowed.
 183          Minimum_Error_C : constant := 7.0;   -- for exp(Complex)
 184          Minimum_Error_I : constant := 2.0;   -- for exp(Imaginary)
 185       begin
 186          Check (Exp (1.0 + 0.0*i) + i,
 187                 E + i,
 188                 "exp(1+0i)",
 189                 Minimum_Error_C);
 190          Check (Exp ((Pi / 2.0) * i) + 1.0,
 191                 1.0 + 1.0*i,
 192                 "exp(pi/2*i)",
 193                 3.6);
 194          Check (Exp (Pi * i) + i,
 195                 -1.0 + 1.0*i,
 196                 "exp(pi*i)",
 197                 5.2);
 198          Check (Exp (Pi * 2.0 * i) + i,
 199                 1.0 + i,
 200                 "exp(2pi*i)",
 201                 8.3);
 202       exception
 203          when Constraint_Error =>
 204             Report.Failed ("Constraint_Error raised in special value test");
 205          when others =>
 206             Report.Failed ("exception in special value test");
 207       end Special_Value_Test;
 208
 209
 210
 211       procedure Exact_Result_Test is
 212          No_Error : constant := 0.0;
 213       begin
 214          -- G.1.2(36);6.0
 215          Check (Exp(0.0 + 0.0*i),  1.0 + 0.0 * i, "exp(0+0i)", No_Error);
 216          Check (Exp(      0.0*i),  1.0 + 0.0 * i, "exp(0i)", No_Error);
 217       exception
 218          when Constraint_Error =>
 219             Report.Failed ("Constraint_Error raised in Exact_Result Test");
 220          when others =>
 221             Report.Failed ("exception in Exact_Result Test");
 222       end Exact_Result_Test;
 223
 224
 225       procedure Identity_Test (A, B : Real) is
 226       -- For this test we use the identity
 227       --    Exp(Z) = Exp(Z-W) * Exp (W)
 228       -- where W = (1+i)/16
 229       --
 230       -- The second part of this test checks the identity
 231       --    Exp(Z) * Exp(-Z)  = 1
 232       --
 233
 234          X, Y : Complex;
 235          Actual1, Actual2 : Complex;
 236          W : constant Complex := (0.0625, 0.0625);
 237          -- the following constant was taken from the CELEFUNC EXP test.
 238          -- This is the value EXP(W) - 1
 239          C : constant Complex := (6.2416044877018563681e-2,
 240                                   6.6487597751003112768e-2);
 241       begin
 242          if Real'Digits > 20 then
 243             -- constant ExpW is accurate to 20 digits.
 244             -- The low bound is 19 * 10**-20
 245             Error_Low_Bound := 0.00000_00000_00019;
 246             Report.Comment ("complex exp accuracy checked to 20 digits");
 247          end if;
 248
 249          Accuracy_Error_Reported := False;  -- reset
 250          for II in 1..Max_Samples loop
 251             X.Re :=  Real'Machine ((B - A) * Real (II) / Real (Max_Samples)
 252                         + A);
 253             for J in 1..Max_Samples loop
 254                X.Im :=  Real'Machine ((B - A) * Real (J) / Real (Max_Samples)
 255                            + A);
 256
 257                Actual1 := Exp(X);
 258
 259                -- Exp(X) = Exp(X-W) * Exp (W)
 260                --        = Exp(X-W) * (1 - (1-Exp(W))
 261                --        = Exp(X-W) * (1 + (Exp(W) - 1))
 262                --        = Exp(X-W) * (1 + C)
 263                Y := X - W;
 264                Actual2 := Exp(Y);
 265                Actual2 := Actual2 + Actual2 * C;
 266
 267                Check (Actual1, Actual2,
 268                       "Identity_1_Test " & Integer'Image (II) &
 269                          Integer'Image (J) & ": Exp((" &
 270                          Real'Image (X.Re) & ", " &
 271                          Real'Image (X.Im) & ")) ",
 272                       20.0);   -- 2 exp and 1 multiply and 1 add = 2*7+1*5+1
 273                     -- Note: The above is not strictly correct, as multiply
 274                     -- has a box error, rather than a relative error.
 275                     -- Supposedly, the interval is chosen to avoid the need
 276                     -- to worry about this.
 277
 278                -- Exp(X) * Exp(-X) + i  = 1 + i
 279                -- The addition of i is to allow a reasonable relative
 280                -- error in the imaginary part
 281                Actual2 := (Actual1 * Exp(-X)) + i;
 282                Check (Actual2, (1.0, 1.0),
 283                       "Identity_2_Test " & Integer'Image (II) &
 284                          Integer'Image (J) & ": Exp((" &
 285                          Real'Image (X.Re) & ", " &
 286                          Real'Image (X.Im) & ")) ",
 287                       20.0);   -- 2 exp and 1 multiply and one add = 2*7+1*5+1
 288
 289                if Accuracy_Error_Reported then
 290                  -- only report the first error in this test in order to keep
 291                  -- lots of failures from producing a huge error log
 292                  return;
 293                end if;
 294             end loop;
 295          end loop;
 296          Error_Low_Bound := 0.0;
 297       exception
 298          when Constraint_Error =>
 299             Report.Failed
 300                ("Constraint_Error raised in Identity_Test" &
 301                 " for X=(" & Real'Image (X.Re) &
 302                 ", " & Real'Image (X.Im) & ")");
 303          when others =>
 304             Report.Failed ("exception in Identity_Test" &
 305                 " for X=(" & Real'Image (X.Re) &
 306                 ", " & Real'Image (X.Im) & ")");
 307       end Identity_Test;
 308
 309
 310
 311       procedure Do_Test is
 312       begin
 313          Special_Value_Test;
 314          Exact_Result_Test;
 315             -- test regions where we can avoid cancellation error problems
 316             -- See Cody page 10.
 317          Identity_Test (0.0625, 1.0);
 318          Identity_Test (15.0, 17.0);
 319          Identity_Test (1.625, 3.0);
 320       end Do_Test;
 321    end Generic_Check;
 322
 323    -----------------------------------------------------------------------
 324    -----------------------------------------------------------------------
 325    package Float_Check is new Generic_Check (Float);
 326
 327    -- check the floating point type with the most digits
 328    type A_Long_Float is digits System.Max_Digits;
 329    package A_Long_Float_Check is new Generic_Check (A_Long_Float);
 330
 331    -----------------------------------------------------------------------
 332    -----------------------------------------------------------------------
 333
 334
 335 begin
 336    Report.Test ("CXG2018",
 337                 "Check the accuracy of the complex EXP function");
 338
 339    if Verbose then
 340       Report.Comment ("checking Standard.Float");
 341    end if;
 342
 343    Float_Check.Do_Test;
 344
 345    if Verbose then
 346       Report.Comment ("checking a digits" &
 347                       Integer'Image (System.Max_Digits) &
 348                       " floating point type");
 349    end if;
 350
 351    A_Long_Float_Check.Do_Test;
 352
 353
 354    Report.Result;
 355 end CXG2018;