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

   1 -- CXG2014.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 SINH and COSH functions return
  28 --      results that are within the error bound allowed.
  29 --
  30 -- TEST DESCRIPTION:
  31 --      This test consists of a generic package that is
  32 --      instantiated to check both Float and a long float type.
  33 --      The test for each floating point type is divided into
  34 --      several parts:
  35 --         Special value checks where the result is a known constant.
  36 --         Checks that use an identity for determining the result.
  37 --         Exception checks.
  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 --      15 Mar 96   SAIC    Initial release for 2.1
  53 --      03 Jun 98   EDS     In line 80, change 1000 to 1024, making it a model
  54 --                          number.  Add Taylor Series terms in line 281.
  55 --      15 Feb 99   RLB     Repaired Subtraction_Error_Test to avoid precision
  56 --                          problems.
  57 --!
  58
  59 --
  60 -- References:
  61 --
  62 -- Software Manual for the Elementary Functions
  63 -- William J. Cody, Jr. and William Waite
  64 -- Prentice-Hall, 1980
  65 --
  66 -- CRC Standard Mathematical Tables
  67 -- 23rd Edition
  68 --
  69 -- Implementation and Testing of Function Software
  70 -- W. J. Cody
  71 -- Problems and Methodologies in Mathematical Software Production
  72 -- editors P. C. Messina and A. Murli
  73 -- Lecture Notes in Computer Science   Volume 142
  74 -- Springer Verlag, 1982
  75 --
  76
  77 with System;
  78 with Report;
  79 with Ada.Numerics.Generic_Elementary_Functions;
  80 procedure CXG2014 is
  81    Verbose : constant Boolean := False;
  82    Max_Samples : constant := 1024;
  83
  84    E  : constant := Ada.Numerics.E;
  85    Cosh1 : constant := (E + 1.0 / E) / 2.0;    -- cosh(1.0)
  86
  87    generic
  88       type Real is digits <>;
  89    package Generic_Check is
  90       procedure Do_Test;
  91    end Generic_Check;
  92
  93    package body Generic_Check is
  94       package Elementary_Functions is new
  95            Ada.Numerics.Generic_Elementary_Functions (Real);
  96       function Sinh (X : Real) return Real renames
  97            Elementary_Functions.Sinh;
  98       function Cosh (X : Real) return Real renames
  99            Elementary_Functions.Cosh;
 100       function Log (X : Real) return Real renames
 101            Elementary_Functions.Log;
 102
 103       -- flag used to terminate some tests early
 104       Accuracy_Error_Reported : Boolean := False;
 105
 106
 107       procedure Check (Actual, Expected : Real;
 108                        Test_Name : String;
 109                        MRE : Real) is
 110          Max_Error : Real;
 111          Rel_Error : Real;
 112          Abs_Error : Real;
 113       begin
 114          -- In the case where the expected result is very small or 0
 115          -- we compute the maximum error as a multiple of Model_Small instead
 116          -- of Model_Epsilon and Expected.
 117          Rel_Error := MRE * abs Expected * Real'Model_Epsilon;
 118          Abs_Error := MRE * Real'Model_Small;
 119          if Rel_Error > Abs_Error then
 120             Max_Error := Rel_Error;
 121          else
 122             Max_Error := Abs_Error;
 123          end if;
 124
 125          if abs (Actual - Expected) > Max_Error then
 126             Accuracy_Error_Reported := True;
 127             Report.Failed (Test_Name &
 128                            " actual: " & Real'Image (Actual) &
 129                            " expected: " & Real'Image (Expected) &
 130                            " difference: " & Real'Image (Actual - Expected) &
 131                            " max err:" & Real'Image (Max_Error) );
 132          elsif Verbose then
 133             if Actual = Expected then
 134                Report.Comment (Test_Name & "  exact result");
 135             else
 136                Report.Comment (Test_Name & "  passed");
 137             end if;
 138          end if;
 139       end Check;
 140
 141
 142       procedure Special_Value_Test is
 143          -- In the following tests the expected result is accurate
 144          -- to the machine precision so the minimum guaranteed error
 145          -- bound can be used.
 146          Minimum_Error : constant := 8.0;
 147       begin
 148          Check (Sinh (1.0),
 149                 (E - 1.0 / E) / 2.0,
 150                 "sinh(1)",
 151                 Minimum_Error);
 152          Check (Cosh (1.0),
 153                 Cosh1,
 154                 "cosh(1)",
 155                 Minimum_Error);
 156          Check (Sinh (2.0),
 157                 (E * E - (1.0 / (E * E))) / 2.0,
 158                 "sinh(2)",
 159                 Minimum_Error);
 160          Check (Cosh (2.0),
 161                 (E * E + (1.0 / (E * E))) / 2.0,
 162                 "cosh(2)",
 163                 Minimum_Error);
 164          Check (Sinh (-1.0),
 165                 (1.0 / E - E) / 2.0,
 166                 "sinh(-1)",
 167                 Minimum_Error);
 168       exception
 169          when Constraint_Error =>
 170             Report.Failed ("Constraint_Error raised in special value test");
 171          when others =>
 172             Report.Failed ("exception in special value test");
 173       end Special_Value_Test;
 174
 175
 176
 177       procedure Exact_Result_Test is
 178          No_Error : constant := 0.0;
 179       begin
 180          -- A.5.1(38);6.0
 181          Check (Sinh (0.0),  0.0, "sinh(0)", No_Error);
 182          Check (Cosh (0.0),  1.0, "cosh(0)", No_Error);
 183       exception
 184          when Constraint_Error =>
 185             Report.Failed ("Constraint_Error raised in Exact_Result Test");
 186          when others =>
 187             Report.Failed ("exception in Exact_Result Test");
 188       end Exact_Result_Test;
 189
 190
 191       procedure Identity_1_Test is
 192       -- For the Sinh test use the identity
 193       --    2 * Sinh(x) * Cosh(1) = Sinh(x+1) + Sinh (x-1)
 194       -- which is transformed to
 195       --    Sinh(x) = ((Sinh(x+1) + Sinh(x-1)) * C
 196       -- where C = 1/(2*Cosh(1))
 197       --
 198       -- For the Cosh test use the identity
 199       --    2 * Cosh(x) * Cosh(1) = Cosh(x+1) + Cosh(x-1)
 200       -- which is transformed to
 201       --    Cosh(x) = C * (Cosh(x+1) + Cosh(x-1))
 202       -- where C is the same as above
 203       --
 204       -- see Cody pg 230-231 for details on the error analysis.
 205       -- The net result is a relative error bound of 16 * Model_Epsilon.
 206
 207          A : constant := 3.0;
 208             -- large upper bound but not so large as to cause Cosh(B)
 209             -- to overflow
 210          B : constant Real := Log(Real'Safe_Last) - 2.0;
 211          X_Minus_1, X, X_Plus_1 : Real;
 212          Actual1, Actual2 : Real;
 213          C : constant := 1.0 / (2.0 * Cosh1);
 214       begin
 215          Accuracy_Error_Reported := False;  -- reset
 216          for I in 1..Max_Samples loop
 217             -- make sure there is no error in x-1, x, and x+1
 218             X_Plus_1 :=  (B - A) * Real (I) / Real (Max_Samples) + A;
 219             X_Plus_1  := Real'Machine (X_Plus_1);
 220             X         := Real'Machine (X_Plus_1 - 1.0);
 221             X_Minus_1 := Real'Machine (X - 1.0);
 222
 223             -- Sinh(x) = ((Sinh(x+1) + Sinh(x-1)) * C
 224             Actual1 := Sinh(X);
 225             Actual2 := C * (Sinh(X_Plus_1) + Sinh(X_Minus_1));
 226
 227             Check (Actual1, Actual2,
 228                    "Identity_1_Test " & Integer'Image (I) & ": sinh(" &
 229                    Real'Image (X) & ") ",
 230                    16.0);
 231
 232             -- Cosh(x) = C * (Cosh(x+1) + Cosh(x-1))
 233             Actual1 := Cosh (X);
 234             Actual2 := C * (Cosh(X_Plus_1) + Cosh (X_Minus_1));
 235             Check (Actual1, Actual2,
 236                    "Identity_1_Test " & Integer'Image (I) & ": cosh(" &
 237                    Real'Image (X) & ") ",
 238                    16.0);
 239
 240             if Accuracy_Error_Reported then
 241               -- only report the first error in this test in order to keep
 242               -- lots of failures from producing a huge error log
 243               return;
 244             end if;
 245
 246          end loop;
 247
 248       exception
 249          when Constraint_Error =>
 250             Report.Failed
 251                ("Constraint_Error raised in Identity_1_Test" &
 252                 " for X=" & Real'Image (X));
 253          when others =>
 254             Report.Failed ("exception in Identity_1_Test" &
 255                 " for X=" & Real'Image (X));
 256       end Identity_1_Test;
 257
 258
 259
 260       procedure Subtraction_Error_Test is
 261       -- This test detects the error resulting from subtraction if
 262       -- the obvious algorithm was used for computing sinh.  That is,
 263       -- it it is computed as (e**x - e**-x)/2.
 264       -- We check the result by using a Taylor series expansion that
 265       -- will produce a result accurate to the machine precision for
 266       -- the range under test.
 267       --
 268       -- The maximum relative error bound for this test is
 269       --  8 for the sinh operation and 7 for the Taylor series
 270       -- for a total of 15 * Model_Epsilon
 271          A : constant := 0.0;
 272          B : constant := 0.5;
 273          X : Real;
 274          X_Squared : Real;
 275          Actual, Expected : Real;
 276       begin
 277          if Real'digits > 15 then
 278              return; -- The approximation below is not accurate beyond
 279                      -- 15 digits. Adding more terms makes the error
 280                      -- larger, so it makes the test worse for more normal
 281                      -- values. Thus, we skip this subtest for larger than
 282                      -- 15 digits.
 283          end if;
 284          Accuracy_Error_Reported := False;  -- reset
 285          for I in 1..Max_Samples loop
 286             X :=  (B - A) * Real (I) / Real (Max_Samples) + A;
 287             X_Squared := X * X;
 288
 289             Actual := Sinh(X);
 290
 291             -- The Taylor series regrouped a bit
 292             Expected :=
 293                X * (1.0 + (X_Squared / 6.0) *
 294                           (1.0 + (X_Squared/20.0) *
 295                                  (1.0 + (X_Squared/42.0) *
 296                                         (1.0 + (X_Squared/72.0) *
 297                                                (1.0 + (X_Squared/110.0) *
 298                                                       (1.0 + (X_Squared/156.0)
 299                    ))))));
 300
 301             Check (Actual, Expected,
 302                    "Subtraction_Error_Test " & Integer'Image (I) & ": sinh(" &
 303                    Real'Image (X) & ") ",
 304                    15.0);
 305
 306             if Accuracy_Error_Reported then
 307               -- only report the first error in this test in order to keep
 308               -- lots of failures from producing a huge error log
 309               return;
 310             end if;
 311
 312          end loop;
 313
 314       exception
 315          when Constraint_Error =>
 316             Report.Failed
 317                ("Constraint_Error raised in Subtraction_Error_Test");
 318          when others =>
 319             Report.Failed ("exception in Subtraction_Error_Test");
 320       end Subtraction_Error_Test;
 321
 322
 323       procedure Exception_Test is
 324          X1, X2 : Real := 0.0;
 325       begin
 326          -- this part of the test is only applicable if 'Machine_Overflows
 327          -- is true.
 328          if Real'Machine_Overflows then
 329
 330             begin
 331               X1 := Sinh (Real'Safe_Last / 2.0);
 332               Report.Failed ("no exception for sinh overflow");
 333             exception
 334                when Constraint_Error => null;
 335                when others =>
 336                   Report.Failed ("wrong exception sinh overflow");
 337             end;
 338
 339             begin
 340               X2 := Cosh (Real'Safe_Last / 2.0);
 341               Report.Failed ("no exception for cosh overflow");
 342             exception
 343                when Constraint_Error => null;
 344                when others =>
 345                   Report.Failed ("wrong exception cosh overflow");
 346             end;
 347
 348          end if;
 349
 350          -- optimizer thwarting
 351          if Report.Ident_Bool (False) then
 352             Report.Comment (Real'Image (X1 + X2));
 353          end if;
 354       end Exception_Test;
 355
 356
 357       procedure Do_Test is
 358       begin
 359          Special_Value_Test;
 360          Exact_Result_Test;
 361          Identity_1_Test;
 362          Subtraction_Error_Test;
 363          Exception_Test;
 364       end Do_Test;
 365    end Generic_Check;
 366
 367    -----------------------------------------------------------------------
 368    -----------------------------------------------------------------------
 369    package Float_Check is new Generic_Check (Float);
 370
 371    -- check the floating point type with the most digits
 372    type A_Long_Float is digits System.Max_Digits;
 373    package A_Long_Float_Check is new Generic_Check (A_Long_Float);
 374
 375    -----------------------------------------------------------------------
 376    -----------------------------------------------------------------------
 377
 378
 379 begin
 380    Report.Test ("CXG2014",
 381                 "Check the accuracy of the SINH and COSH functions");
 382
 383    if Verbose then
 384       Report.Comment ("checking Standard.Float");
 385    end if;
 386
 387    Float_Check.Do_Test;
 388
 389    if Verbose then
 390       Report.Comment ("checking a digits" &
 391                       Integer'Image (System.Max_Digits) &
 392                       " floating point type");
 393    end if;
 394
 395    A_Long_Float_Check.Do_Test;
 396
 397
 398    Report.Result;
 399 end CXG2014;