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 92with System; 93with Report; 94with Ada.Numerics.Generic_Elementary_Functions; 95with Ada.Numerics.Elementary_Functions; 96procedure 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 480begin 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; 499end CXG2004; 500