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 59with System; 60with Report; 61with Ada.Numerics.Generic_Complex_Types; 62with Ada.Numerics.Generic_Elementary_Functions; 63procedure 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 401begin 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; 421end CXG2009; 422