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

   1 -- CXG2004.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 sin and cos 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 --      the following parts:
  35 --         Special value checks where the result is a known constant.
  36 --         Checks using an identity relationship.
  37 --
  38 -- SPECIAL REQUIREMENTS
  39 --      The Strict Mode for the numerical accuracy must be
  40 --      selected.  The method by which this mode is selected
  41 --      is implementation dependent.
  42 --
  43 -- APPLICABILITY CRITERIA:
  44 --      This test applies only to implementations supporting the
  45 --      Numerics Annex.
  46 --      This test only applies to the Strict Mode for numerical
  47 --      accuracy.
  48 --
  49 --
  50 -- CHANGE HISTORY:
  51 --      13 FEB 96   SAIC    Initial release for 2.1
  52 --      22 APR 96   SAIC    Changed to generic implementation.
  53 --      18 AUG 96   SAIC    Improvements to commentary.
  54 --      23 OCT 96   SAIC    Exact results are not required unless the
  55 --                          cycle is specified.
  56 --      28 FEB 97   PWB.CTA Removed checks where cycle 2.0*Pi is specified
  57 --      02 JUN 98   EDS     Revised calculations to ensure that X is exactly
  58 --                          three times Y per advice of numerics experts.
  59 --
  60 -- CHANGE NOTE:
  61 --      According to Ken Dritz, author of the Numerics Annex of the RM,
  62 --      one should never specify the cycle 2.0*Pi for the trigonometric
  63 --      functions.  In particular, if the machine number for the first
  64 --      argument is not an exact multiple of the machine number for the
  65 --      explicit cycle, then the specified exact results cannot be
  66 --      reasonably expected.  The affected checks in this test have been
  67 --      marked as comments, with the additional notation "pwb-math".
  68 --      Phil Brashear
  69 --!
  70
  71 --
  72 -- References:
  73 --
  74 -- Software Manual for the Elementary Functions
  75 -- William J. Cody, Jr. and William Waite
  76 -- Prentice-Hall, 1980
  77 --
  78 -- CRC Standard Mathematical Tables
  79 -- 23rd Edition
  80 --
  81 -- Implementation and Testing of Function Software
  82 -- W. J. Cody
  83 -- Problems and Methodologies in Mathematical Software Production
  84 -- editors P. C. Messina and A. Murli
  85 -- Lecture Notes in Computer Science   Volume 142
  86 -- Springer Verlag, 1982
  87 --
  88 -- The sin and cos checks are translated directly from
  89 -- the netlib FORTRAN code that was written by W. Cody.
  90 --
  91
  92 with System;
  93 with Report;
  94 with Ada.Numerics.Generic_Elementary_Functions;
  95 with Ada.Numerics.Elementary_Functions;
  96 procedure CXG2004 is
  97    Verbose : constant Boolean := False;
  98    Number_Samples : constant := 1000;
  99
 100    -- CRC Standard Mathematical Tables;  23rd Edition; pg 738
 101    Sqrt2 : constant :=
 102         1.41421_35623_73095_04880_16887_24209_69807_85696_71875_37695;
 103    Sqrt3 : constant :=
 104         1.73205_08075_68877_29352_74463_41505_87236_69428_05253_81039;
 105
 106    Pi : constant := Ada.Numerics.Pi;
 107
 108    generic
 109       type Real is digits <>;
 110    package Generic_Check is
 111       procedure Do_Test;
 112    end Generic_Check;
 113
 114    package body Generic_Check is
 115       package Elementary_Functions is new
 116            Ada.Numerics.Generic_Elementary_Functions (Real);
 117
 118       function Sin (X : Real) return Real renames
 119            Elementary_Functions.Sin;
 120       function Cos (X : Real) return Real renames
 121            Elementary_Functions.Cos;
 122       function Sin (X, Cycle : Real) return Real renames
 123            Elementary_Functions.Sin;
 124       function Cos (X, Cycle : Real) return Real renames
 125            Elementary_Functions.Cos;
 126
 127       Accuracy_Error_Reported : Boolean := False;
 128
 129       procedure Check (Actual, Expected : Real;
 130                        Test_Name : String;
 131                        MRE : Real) is
 132          Rel_Error,
 133          Abs_Error,
 134          Max_Error : Real;
 135       begin
 136
 137          -- In the case where the expected result is very small or 0
 138          -- we compute the maximum error as a multiple of Model_Epsilon instead
 139          -- of Model_Epsilon and Expected.
 140          Rel_Error := MRE * abs Expected * Real'Model_Epsilon;
 141          Abs_Error := MRE * Real'Model_Epsilon;
 142          if Rel_Error > Abs_Error then
 143             Max_Error := Rel_Error;
 144          else
 145             Max_Error := Abs_Error;
 146          end if;
 147
 148
 149          -- in addition to the relative error checks we apply the
 150          -- criteria of G.2.4(16)
 151          if abs (Actual) > 1.0 then
 152             Accuracy_Error_Reported := True;
 153             Report.Failed (Test_Name & " result > 1.0");
 154          elsif abs (Actual - Expected) > Max_Error then
 155             Accuracy_Error_Reported := True;
 156             Report.Failed (Test_Name &
 157                            " actual: " & Real'Image (Actual) &
 158                            " expected: " & Real'Image (Expected) &
 159                            " difference: " &
 160                            Real'Image (Actual - Expected) &
 161                            " mre:" &
 162                            Real'Image (Max_Error) );
 163          elsif Verbose then
 164             if Actual = Expected then
 165                Report.Comment (Test_Name & "  exact result");
 166             else
 167                Report.Comment (Test_Name & "  passed");
 168             end if;
 169          end if;
 170       end Check;
 171
 172
 173       procedure Sin_Check (A, B : Real;
 174                            Arg_Range : String) is
 175          -- test a selection of
 176          -- arguments selected from the range A to B.
 177          --
 178          -- This test uses the identity
 179          --   sin(x) = sin(x/3)*(3 - 4 * sin(x/3)**2)
 180          --
 181          -- Note that in this test we must take into account the
 182          -- error in the calculation of the expected result so
 183          -- the maximum relative error is larger than the
 184          -- accuracy required by the ARM.
 185
 186          X, Y, ZZ : Real;
 187          Actual, Expected : Real;
 188          MRE : Real;
 189          Ran : Real;
 190       begin
 191          Accuracy_Error_Reported := False;  -- reset
 192          for I in 1 .. Number_Samples loop
 193             -- Evenly distributed selection of arguments
 194             Ran := Real (I) / Real (Number_Samples);
 195
 196             -- make sure x and x/3 are both exactly representable
 197             -- on the machine.  See "Implementation and Testing of
 198             -- Function Software" page 44.
 199             X := (B - A) * Ran + A;
 200             Y := Real'Leading_Part
 201                       ( X/3.0,
 202                         Real'Machine_Mantissa - Real'Exponent (3.0) );
 203             X := Y * 3.0;
 204
 205             Actual := Sin (X);
 206
 207             ZZ := Sin(Y);
 208             Expected := ZZ * (3.0 - 4.0 * ZZ * ZZ);
 209
 210             -- note that since the expected value is computed, we
 211             -- must take the error in that computation into account.
 212             -- See Cody pp 139-141.
 213             MRE := 4.0;
 214
 215             Check (Actual, Expected,
 216                    "sin test of range" & Arg_Range &
 217                       Integer'Image (I),
 218                    MRE);
 219             exit when Accuracy_Error_Reported;
 220          end loop;
 221       exception
 222          when Constraint_Error =>
 223             Report.Failed
 224                ("Constraint_Error raised in sin check");
 225          when others =>
 226             Report.Failed ("exception in sin check");
 227       end Sin_Check;
 228
 229
 230
 231       procedure Cos_Check (A, B : Real;
 232                                   Arg_Range : String) is
 233          -- test a selection of
 234          -- arguments selected from the range A to B.
 235          --
 236          -- This test uses the identity
 237          --   cos(x) = cos(x/3)*(4 * cos(x/3)**2 - 3)
 238          --
 239          -- Note that in this test we must take into account the
 240          -- error in the calculation of the expected result so
 241          -- the maximum relative error is larger than the
 242          -- accuracy required by the ARM.
 243
 244          X, Y, ZZ : Real;
 245          Actual, Expected : Real;
 246          MRE : Real;
 247          Ran : Real;
 248       begin
 249          Accuracy_Error_Reported := False;  -- reset
 250          for I in 1 .. Number_Samples loop
 251             -- Evenly distributed selection of arguments
 252             Ran := Real (I) / Real (Number_Samples);
 253
 254             -- make sure x and x/3 are both exactly representable
 255             -- on the machine.  See "Implementation and Testing of
 256             -- Function Software" page 44.
 257             X := (B - A) * Ran + A;
 258             Y := Real'Leading_Part
 259                       ( X/3.0,
 260                         Real'Machine_Mantissa - Real'Exponent (3.0) );
 261             X := Y * 3.0;
 262
 263             Actual := Cos (X);
 264
 265             ZZ := Cos(Y);
 266             Expected := ZZ * (4.0 * ZZ * ZZ - 3.0);
 267
 268             -- note that since the expected value is computed, we
 269             -- must take the error in that computation into account.
 270             -- See Cody pp 141-143.
 271             MRE := 6.0;
 272
 273             Check (Actual, Expected,
 274                    "cos test of range" & Arg_Range &
 275                       Integer'Image (I),
 276                    MRE);
 277             exit when Accuracy_Error_Reported;
 278          end loop;
 279       exception
 280          when Constraint_Error =>
 281             Report.Failed
 282                ("Constraint_Error raised in cos check");
 283          when others =>
 284             Report.Failed ("exception in cos check");
 285       end Cos_Check;
 286
 287
 288       procedure Special_Angle_Checks is
 289          type Data_Point is
 290             record
 291                Degrees,
 292                Radians,
 293                Sine,
 294                Cosine         : Real;
 295                Sin_Result_Error,
 296                Cos_Result_Error : Boolean;
 297             end record;
 298
 299          type Test_Data_Type is array (Positive range <>) of Data_Point;
 300
 301          -- the values in the following table only involve static
 302          -- expressions to minimize any loss of precision.  However,
 303          -- there are two sources of error that must be accounted for
 304          -- in the following tests.
 305          -- First, when a cycle is not specified there can be a roundoff
 306          -- error in the value of Pi used.  This error does not apply
 307          -- when a cycle of 2.0 * Pi is explicitly provided.
 308          -- Second, the expected results that involve sqrt values also
 309          -- have a potential roundoff error.
 310          -- The amount of error due to error in the argument is computed
 311          -- as follows:
 312          --   sin(x+err) = sin(x)*cos(err) + cos(x)*sin(err)
 313          --              ~= sin(x) + err * cos(x)
 314          -- similarly for cos the error due to error in the argument is
 315          -- computed as follows:
 316          --   cos(x+err) = cos(x)*cos(err) - sin(x)*sin(err)
 317          --              ~= cos(x) - err * sin(x)
 318          -- In both cases the term "err" is bounded by 0.5 * argument.
 319
 320          Test_Data : constant Test_Data_Type := (
 321 --  degrees      radians          sine       cosine  sin_er cos_er    test #
 322   (  0.0,           0.0,          0.0,         1.0,  False, False ),    -- 1
 323   ( 30.0,        Pi/6.0,          0.5,   Sqrt3/2.0,  False, True  ),    -- 2
 324   ( 60.0,        Pi/3.0,    Sqrt3/2.0,         0.5,  True,  False ),    -- 3
 325   ( 90.0,        Pi/2.0,          1.0,         0.0,  False, False ),    -- 4
 326   (120.0,    2.0*Pi/3.0,    Sqrt3/2.0,        -0.5,  True,  False ),    -- 5
 327   (150.0,    5.0*Pi/6.0,          0.5,  -Sqrt3/2.0,  False, True  ),    -- 6
 328   (180.0,            Pi,          0.0,        -1.0,  False, False ),    -- 7
 329   (210.0,    7.0*Pi/6.0,         -0.5,  -Sqrt3/2.0,  False, True  ),    -- 8
 330   (240.0,    8.0*Pi/6.0,   -Sqrt3/2.0,        -0.5,  True,  False ),    -- 9
 331   (270.0,    9.0*Pi/6.0,         -1.0,         0.0,  False, False ),    -- 10
 332   (300.0,   10.0*Pi/6.0,   -Sqrt3/2.0,         0.5,  True,  False ),    -- 11
 333   (330.0,   11.0*Pi/6.0,         -0.5,   Sqrt3/2.0,  False, True  ),    -- 12
 334   (360.0,        2.0*Pi,          0.0,         1.0,  False, False ),    -- 13
 335   ( 45.0,        Pi/4.0,    Sqrt2/2.0,   Sqrt2/2.0,  True,  True  ),    -- 14
 336   (135.0,    3.0*Pi/4.0,    Sqrt2/2.0,  -Sqrt2/2.0,  True,  True  ),    -- 15
 337   (225.0,    5.0*Pi/4.0,   -Sqrt2/2.0,  -Sqrt2/2.0,  True,  True  ),    -- 16
 338   (315.0,    7.0*Pi/4.0,   -Sqrt2/2.0,   Sqrt2/2.0,  True,  True  ),    -- 17
 339   (405.0,    9.0*Pi/4.0,    Sqrt2/2.0,   Sqrt2/2.0,  True,  True  ) );  -- 18
 340
 341
 342          Y : Real;
 343          Sin_Arg_Err,
 344          Cos_Arg_Err,
 345          Sin_Result_Err,
 346          Cos_Result_Err  : Real;
 347       begin
 348          for I in Test_Data'Range loop
 349             -- compute error components
 350             Sin_Arg_Err := abs Test_Data (I).Cosine *
 351                            abs Test_Data (I).Radians / 2.0;
 352             Cos_Arg_Err := abs Test_Data (I).Sine *
 353                            abs Test_Data (I).Radians / 2.0;
 354
 355             if Test_Data (I).Sin_Result_Error then
 356                Sin_Result_Err := 0.5;
 357             else
 358                Sin_Result_Err := 0.0;
 359             end if;
 360
 361             if Test_Data (I).Cos_Result_Error then
 362                Cos_Result_Err := 1.0;
 363             else
 364                Cos_Result_Err := 0.0;
 365             end if;
 366
 367
 368
 369             Y := Sin (Test_Data (I).Radians);
 370             Check (Y, Test_Data (I).Sine,
 371                    "test" & Integer'Image (I) & " sin(r)",
 372                    2.0 + Sin_Arg_Err + Sin_Result_Err);
 373             Y := Cos (Test_Data (I).Radians);
 374             Check (Y, Test_Data (I).Cosine,
 375                    "test" & Integer'Image (I) & " cos(r)",
 376                    2.0 + Cos_Arg_Err + Cos_Result_Err);
 377             Y := Sin (Test_Data (I).Degrees, 360.0);
 378             Check (Y, Test_Data (I).Sine,
 379                    "test" & Integer'Image (I) & " sin(d,360)",
 380                    2.0 + Sin_Result_Err);
 381             Y := Cos (Test_Data (I).Degrees, 360.0);
 382             Check (Y, Test_Data (I).Cosine,
 383                    "test" & Integer'Image (I) & " cos(d,360)",
 384                    2.0 + Cos_Result_Err);
 385 --pwb-math            Y := Sin (Test_Data (I).Radians, 2.0*Pi);
 386 --pwb-math            Check (Y, Test_Data (I).Sine,
 387 --pwb-math                   "test" & Integer'Image (I) & " sin(r,2pi)",
 388 --pwb-math                   2.0 + Sin_Result_Err);
 389 --pwb-math            Y := Cos (Test_Data (I).Radians, 2.0*Pi);
 390 --pwb-math            Check (Y, Test_Data (I).Cosine,
 391 --pwb-math                   "test" & Integer'Image (I) & " cos(r,2pi)",
 392 --pwb-math                   2.0 + Cos_Result_Err);
 393          end loop;
 394       exception
 395          when Constraint_Error =>
 396             Report.Failed ("Constraint_Error raised in special angle test");
 397          when others =>
 398             Report.Failed ("exception in special angle test");
 399       end Special_Angle_Checks;
 400
 401
 402       -- check the rule of A.5.1(41);6.0 which requires that the
 403       -- result be exact if the mathematical result is 0.0, 1.0,
 404       -- or -1.0
 405       procedure Exact_Result_Checks is
 406          type Data_Point is
 407             record
 408                Degrees,
 409                Sine,
 410                Cosine         : Real;
 411             end record;
 412
 413          type Test_Data_Type is array (Positive range <>) of Data_Point;
 414          Test_Data : constant Test_Data_Type := (
 415             -- degrees       sine       cosine       test #
 416               (  0.0,         0.0,         1.0  ),    -- 1
 417               ( 90.0,         1.0,         0.0  ),    -- 2
 418               (180.0,         0.0,        -1.0  ),    -- 3
 419               (270.0,        -1.0,         0.0  ),    -- 4
 420               (360.0,         0.0,         1.0  ),    -- 5
 421               ( 90.0 + 360.0, 1.0,         0.0  ),    -- 6
 422               (180.0 + 360.0, 0.0,        -1.0  ),    -- 7
 423               (270.0 + 360.0,-1.0,         0.0  ),    -- 8
 424               (360.0 + 360.0, 0.0,         1.0  ) );  -- 9
 425
 426          Y : Real;
 427       begin
 428          for I in Test_Data'Range loop
 429             Y := Sin (Test_Data(I).Degrees, 360.0);
 430             if Y /= Test_Data(I).Sine then
 431                Report.Failed ("exact result for sin(" &
 432                   Real'Image (Test_Data(I).Degrees) &
 433                   ", 360.0) is not" &
 434                   Real'Image (Test_Data(I).Sine) &
 435                   "  Difference is " &
 436                   Real'Image (Y - Test_Data(I).Sine) );
 437             end if;
 438
 439             Y := Cos (Test_Data(I).Degrees, 360.0);
 440             if Y /= Test_Data(I).Cosine then
 441                Report.Failed ("exact result for cos(" &
 442                   Real'Image (Test_Data(I).Degrees) &
 443                   ", 360.0) is not" &
 444                   Real'Image (Test_Data(I).Cosine) &
 445                   "  Difference is " &
 446                   Real'Image (Y - Test_Data(I).Cosine) );
 447             end if;
 448          end loop;
 449       exception
 450          when Constraint_Error =>
 451             Report.Failed ("Constraint_Error raised in exact result check");
 452          when others =>
 453             Report.Failed ("exception in exact result check");
 454       end Exact_Result_Checks;
 455
 456
 457       procedure Do_Test is
 458       begin
 459          Special_Angle_Checks;
 460          Sin_Check (0.0, Pi/2.0, "0..pi/2");
 461          Sin_Check (6.0*Pi, 6.5*Pi, "6pi..6.5pi");
 462          Cos_Check (7.0*Pi, 7.5*Pi, "7pi..7.5pi");
 463          Exact_Result_Checks;
 464       end Do_Test;
 465    end Generic_Check;
 466
 467    -----------------------------------------------------------------------
 468    -----------------------------------------------------------------------
 469
 470    package Float_Check is new Generic_Check (Float);
 471
 472    -- check the floating point type with the most digits
 473    type A_Long_Float is digits System.Max_Digits;
 474    package A_Long_Float_Check is new Generic_Check (A_Long_Float);
 475
 476    -----------------------------------------------------------------------
 477    -----------------------------------------------------------------------
 478
 479
 480 begin
 481    Report.Test ("CXG2004",
 482                 "Check the accuracy of the sin and cos functions");
 483
 484    if Verbose then
 485       Report.Comment ("checking Standard.Float");
 486    end if;
 487
 488    Float_Check.Do_Test;
 489
 490    if Verbose then
 491       Report.Comment ("checking a digits" &
 492                       Integer'Image (System.Max_Digits) &
 493                       " floating point type");
 494    end if;
 495
 496    A_Long_Float_Check.Do_Test;
 497
 498    Report.Result;
 499 end CXG2004;