| -- CXG2021.A |
| -- |
| -- Grant of Unlimited Rights |
| -- |
| -- Under contracts F33600-87-D-0337, F33600-84-D-0280, MDA903-79-C-0687, |
| -- F08630-91-C-0015, and DCA100-97-D-0025, the U.S. Government obtained |
| -- unlimited rights in the software and documentation contained herein. |
| -- Unlimited rights are defined in DFAR 252.227-7013(a)(19). By making |
| -- this public release, the Government intends to confer upon all |
| -- recipients unlimited rights equal to those held by the Government. |
| -- These rights include rights to use, duplicate, release or disclose the |
| -- released technical data and computer software in whole or in part, in |
| -- any manner and for any purpose whatsoever, and to have or permit others |
| -- to do so. |
| -- |
| -- DISCLAIMER |
| -- |
| -- ALL MATERIALS OR INFORMATION HEREIN RELEASED, MADE AVAILABLE OR |
| -- DISCLOSED ARE AS IS. THE GOVERNMENT MAKES NO EXPRESS OR IMPLIED |
| -- WARRANTY AS TO ANY MATTER WHATSOEVER, INCLUDING THE CONDITIONS OF THE |
| -- SOFTWARE, DOCUMENTATION OR OTHER INFORMATION RELEASED, MADE AVAILABLE |
| -- OR DISCLOSED, OR THE OWNERSHIP, MERCHANTABILITY, OR FITNESS FOR A |
| -- PARTICULAR PURPOSE OF SAID MATERIAL. |
| --* |
| -- |
| -- OBJECTIVE: |
| -- Check that the complex SIN and COS functions return |
| -- a result that is within the error bound allowed. |
| -- |
| -- TEST DESCRIPTION: |
| -- This test consists of a generic package that is |
| -- instantiated to check complex numbers based upon |
| -- both Float and a long float type. |
| -- The test for each floating point type is divided into |
| -- several parts: |
| -- Special value checks where the result is a known constant. |
| -- Checks that use an identity for determining the result. |
| -- |
| -- SPECIAL REQUIREMENTS |
| -- The Strict Mode for the numerical accuracy must be |
| -- selected. The method by which this mode is selected |
| -- is implementation dependent. |
| -- |
| -- APPLICABILITY CRITERIA: |
| -- This test applies only to implementations supporting the |
| -- Numerics Annex. |
| -- This test only applies to the Strict Mode for numerical |
| -- accuracy. |
| -- |
| -- |
| -- CHANGE HISTORY: |
| -- 27 Mar 96 SAIC Initial release for 2.1 |
| -- 22 Aug 96 SAIC No longer skips test for systems with |
| -- more than 20 digits of precision. |
| -- |
| --! |
| |
| -- |
| -- References: |
| -- |
| -- W. J. Cody |
| -- CELEFUNT: A Portable Test Package for Complex Elementary Functions |
| -- Algorithm 714, Collected Algorithms from ACM. |
| -- Published in Transactions On Mathematical Software, |
| -- Vol. 19, No. 1, March, 1993, pp. 1-21. |
| -- |
| -- CRC Standard Mathematical Tables |
| -- 23rd Edition |
| -- |
| |
| with System; |
| with Report; |
| with Ada.Numerics.Generic_Complex_Types; |
| with Ada.Numerics.Generic_Complex_Elementary_Functions; |
| procedure CXG2021 is |
| Verbose : constant Boolean := False; |
| -- Note that Max_Samples is the number of samples taken in |
| -- both the real and imaginary directions. Thus, for Max_Samples |
| -- of 100 the number of values checked is 10000. |
| Max_Samples : constant := 100; |
| |
| E : constant := Ada.Numerics.E; |
| Pi : constant := Ada.Numerics.Pi; |
| |
| generic |
| type Real is digits <>; |
| package Generic_Check is |
| procedure Do_Test; |
| end Generic_Check; |
| |
| package body Generic_Check is |
| package Complex_Type is new |
| Ada.Numerics.Generic_Complex_Types (Real); |
| use Complex_Type; |
| |
| package CEF is new |
| Ada.Numerics.Generic_Complex_Elementary_Functions (Complex_Type); |
| |
| function Sin (X : Complex) return Complex renames CEF.Sin; |
| function Cos (X : Complex) return Complex renames CEF.Cos; |
| |
| -- flag used to terminate some tests early |
| Accuracy_Error_Reported : Boolean := False; |
| |
| -- The following value is a lower bound on the accuracy |
| -- required. It is normally 0.0 so that the lower bound |
| -- is computed from Model_Epsilon. However, for tests |
| -- where the expected result is only known to a certain |
| -- amount of precision this bound takes on a non-zero |
| -- value to account for that level of precision. |
| Error_Low_Bound : Real := 0.0; |
| |
| -- the E_Factor is an additional amount added to the Expected |
| -- value prior to computing the maximum relative error. |
| -- This is needed because the error analysis (Cody pg 17-20) |
| -- requires this additional allowance. |
| procedure Check (Actual, Expected : Real; |
| Test_Name : String; |
| MRE : Real; |
| E_Factor : Real := 0.0) is |
| Max_Error : Real; |
| Rel_Error : Real; |
| Abs_Error : Real; |
| begin |
| -- In the case where the expected result is very small or 0 |
| -- we compute the maximum error as a multiple of Model_Epsilon instead |
| -- of Model_Epsilon and Expected. |
| Rel_Error := MRE * Real'Model_Epsilon * (abs Expected + E_Factor); |
| Abs_Error := MRE * Real'Model_Epsilon; |
| if Rel_Error > Abs_Error then |
| Max_Error := Rel_Error; |
| else |
| Max_Error := Abs_Error; |
| end if; |
| |
| -- take into account the low bound on the error |
| if Max_Error < Error_Low_Bound then |
| Max_Error := Error_Low_Bound; |
| end if; |
| |
| if abs (Actual - Expected) > Max_Error then |
| Accuracy_Error_Reported := True; |
| Report.Failed (Test_Name & |
| " actual: " & Real'Image (Actual) & |
| " expected: " & Real'Image (Expected) & |
| " difference: " & Real'Image (Actual - Expected) & |
| " max err:" & Real'Image (Max_Error) & |
| " efactor:" & Real'Image (E_Factor) ); |
| elsif Verbose then |
| if Actual = Expected then |
| Report.Comment (Test_Name & " exact result"); |
| else |
| Report.Comment (Test_Name & " passed" & |
| " actual: " & Real'Image (Actual) & |
| " expected: " & Real'Image (Expected) & |
| " difference: " & Real'Image (Actual - Expected) & |
| " max err:" & Real'Image (Max_Error) & |
| " efactor:" & Real'Image (E_Factor) ); |
| end if; |
| end if; |
| end Check; |
| |
| |
| procedure Check (Actual, Expected : Complex; |
| Test_Name : String; |
| MRE : Real; |
| R_Factor, I_Factor : Real := 0.0) is |
| begin |
| Check (Actual.Re, Expected.Re, Test_Name & " real part", |
| MRE, R_Factor); |
| Check (Actual.Im, Expected.Im, Test_Name & " imaginary part", |
| MRE, I_Factor); |
| end Check; |
| |
| |
| procedure Special_Value_Test is |
| -- In the following tests the expected result is accurate |
| -- to the machine precision so the minimum guaranteed error |
| -- bound can be used if the argument is exact. |
| -- Since the argument involves Pi, we must allow for this |
| -- inexact argument. |
| Minimum_Error : constant := 11.0; |
| begin |
| Check (Sin (Pi/2.0 + 0.0*i), |
| 1.0 + 0.0*i, |
| "sin(pi/2+0i)", |
| Minimum_Error + 1.0); |
| Check (Cos (Pi/2.0 + 0.0*i), |
| 0.0 + 0.0*i, |
| "cos(pi/2+0i)", |
| Minimum_Error + 1.0); |
| exception |
| when Constraint_Error => |
| Report.Failed ("Constraint_Error raised in special value test"); |
| when others => |
| Report.Failed ("exception in special value test"); |
| end Special_Value_Test; |
| |
| |
| |
| procedure Exact_Result_Test is |
| No_Error : constant := 0.0; |
| begin |
| -- G.1.2(36);6.0 |
| Check (Sin(0.0 + 0.0*i), 0.0 + 0.0 * i, "sin(0+0i)", No_Error); |
| Check (Cos(0.0 + 0.0*i), 1.0 + 0.0 * i, "cos(0+0i)", No_Error); |
| exception |
| when Constraint_Error => |
| Report.Failed ("Constraint_Error raised in Exact_Result Test"); |
| when others => |
| Report.Failed ("exception in Exact_Result Test"); |
| end Exact_Result_Test; |
| |
| |
| procedure Identity_Test (RA, RB, IA, IB : Real) is |
| -- Tests an identity over a range of values specified |
| -- by the 4 parameters. RA and RB denote the range for the |
| -- real part while IA and IB denote the range for the |
| -- imaginary part. |
| -- |
| -- For this test we use the identity |
| -- Sin(Z) = Sin(Z-W) * Cos(W) + Cos(Z-W) * Sin(W) |
| -- and |
| -- Cos(Z) = Cos(Z-W) * Cos(W) - Sin(Z-W) * Sin(W) |
| -- |
| |
| X, Y : Real; |
| Z : Complex; |
| W : constant Complex := Compose_From_Cartesian(0.0625, 0.0625); |
| ZmW : Complex; -- Z - W |
| Sin_ZmW, |
| Cos_ZmW : Complex; |
| Actual1, Actual2 : Complex; |
| R_Factor : Real; -- additional real error factor |
| I_Factor : Real; -- additional imaginary error factor |
| Sin_W : constant Complex := (6.2581348413276935585E-2, |
| 6.2418588008436587236E-2); |
| -- numeric stability is enhanced by using Cos(W) - 1.0 instead of |
| -- Cos(W) in the computation. |
| Cos_W_m_1 : constant Complex := (-2.5431314180235545803E-6, |
| -3.9062493377261771826E-3); |
| |
| |
| begin |
| if Real'Digits > 20 then |
| -- constants used here accurate to 20 digits. Allow 1 |
| -- additional digit of error for computation. |
| Error_Low_Bound := 0.00000_00000_00000_0001; |
| Report.Comment ("accuracy checked to 19 digits"); |
| end if; |
| |
| Accuracy_Error_Reported := False; -- reset |
| for II in 0..Max_Samples loop |
| X := (RB - RA) * Real (II) / Real (Max_Samples) + RA; |
| for J in 0..Max_Samples loop |
| Y := (IB - IA) * Real (J) / Real (Max_Samples) + IA; |
| |
| Z := Compose_From_Cartesian(X,Y); |
| ZmW := Z - W; |
| Sin_ZmW := Sin (ZmW); |
| Cos_ZmW := Cos (ZmW); |
| |
| -- now for the first identity |
| -- Sin(Z) = Sin(Z-W) * Cos(W) + Cos(Z-W) * Sin(W) |
| -- = Sin(Z-W) * (1+(Cos(W)-1)) + Cos(Z-W) * Sin(W) |
| -- = Sin(Z-W) + Sin(Z-W)*(Cos(W)-1) + Cos(Z-W)*Sin(W) |
| |
| |
| Actual1 := Sin (Z); |
| Actual2 := Sin_ZmW + (Sin_ZmW * Cos_W_m_1 + Cos_ZmW * Sin_W); |
| |
| -- The computation of the additional error factors are taken |
| -- from Cody pages 17-20. |
| |
| R_Factor := abs (Re (Sin_ZmW) * Re (1.0 - Cos_W_m_1)) + |
| abs (Im (Sin_ZmW) * Im (1.0 - Cos_W_m_1)) + |
| abs (Re (Cos_ZmW) * Re (Sin_W)) + |
| abs (Re (Cos_ZmW) * Re (1.0 - Cos_W_m_1)); |
| |
| I_Factor := abs (Re (Sin_ZmW) * Im (1.0 - Cos_W_m_1)) + |
| abs (Im (Sin_ZmW) * Re (1.0 - Cos_W_m_1)) + |
| abs (Re (Cos_ZmW) * Im (Sin_W)) + |
| abs (Im (Cos_ZmW) * Re (1.0 - Cos_W_m_1)); |
| |
| Check (Actual1, Actual2, |
| "Identity_1_Test " & Integer'Image (II) & |
| Integer'Image (J) & ": Sin((" & |
| Real'Image (Z.Re) & ", " & |
| Real'Image (Z.Im) & ")) ", |
| 11.0, R_Factor, I_Factor); |
| |
| -- now for the second identity |
| -- Cos(Z) = Cos(Z-W) * Cos(W) - Sin(Z-W) * Sin(W) |
| -- = Cos(Z-W) * (1+(Cos(W)-1) - Sin(Z-W) * Sin(W) |
| Actual1 := Cos (Z); |
| Actual2 := Cos_ZmW + (Cos_ZmW * Cos_W_m_1 - Sin_ZmW * Sin_W); |
| |
| -- The computation of the additional error factors are taken |
| -- from Cody pages 17-20. |
| |
| R_Factor := abs (Re (Sin_ZmW) * Re (Sin_W)) + |
| abs (Im (Sin_ZmW) * Im (Sin_W)) + |
| abs (Re (Cos_ZmW) * Re (1.0 - Cos_W_m_1)) + |
| abs (Im (Cos_ZmW) * Im (1.0 - Cos_W_m_1)); |
| |
| I_Factor := abs (Re (Sin_ZmW) * Im (Sin_W)) + |
| abs (Im (Sin_ZmW) * Re (Sin_W)) + |
| abs (Re (Cos_ZmW) * Im (1.0 - Cos_W_m_1)) + |
| abs (Im (Cos_ZmW) * Re (1.0 - Cos_W_m_1)); |
| |
| Check (Actual1, Actual2, |
| "Identity_2_Test " & Integer'Image (II) & |
| Integer'Image (J) & ": Cos((" & |
| Real'Image (Z.Re) & ", " & |
| Real'Image (Z.Im) & ")) ", |
| 11.0, R_Factor, I_Factor); |
| |
| if Accuracy_Error_Reported then |
| -- only report the first error in this test in order to keep |
| -- lots of failures from producing a huge error log |
| Error_Low_Bound := 0.0; -- reset |
| return; |
| end if; |
| end loop; |
| end loop; |
| |
| Error_Low_Bound := 0.0; -- reset |
| exception |
| when Constraint_Error => |
| Report.Failed |
| ("Constraint_Error raised in Identity_Test" & |
| " for Z=(" & Real'Image (X) & |
| ", " & Real'Image (Y) & ")"); |
| when others => |
| Report.Failed ("exception in Identity_Test" & |
| " for Z=(" & Real'Image (X) & |
| ", " & Real'Image (Y) & ")"); |
| end Identity_Test; |
| |
| |
| procedure Do_Test is |
| begin |
| Special_Value_Test; |
| Exact_Result_Test; |
| -- test regions where sin and cos have the same sign and |
| -- about the same magnitude. This will minimize subtraction |
| -- errors in the identities. |
| -- See Cody page 17. |
| Identity_Test (0.0625, 10.0, 0.0625, 10.0); |
| Identity_Test ( 16.0, 17.0, 16.0, 17.0); |
| end Do_Test; |
| end Generic_Check; |
| |
| ----------------------------------------------------------------------- |
| ----------------------------------------------------------------------- |
| package Float_Check is new Generic_Check (Float); |
| |
| -- check the floating point type with the most digits |
| type A_Long_Float is digits System.Max_Digits; |
| package A_Long_Float_Check is new Generic_Check (A_Long_Float); |
| |
| ----------------------------------------------------------------------- |
| ----------------------------------------------------------------------- |
| |
| |
| begin |
| Report.Test ("CXG2021", |
| "Check the accuracy of the complex SIN and COS functions"); |
| |
| if Verbose then |
| Report.Comment ("checking Standard.Float"); |
| end if; |
| |
| Float_Check.Do_Test; |
| |
| if Verbose then |
| Report.Comment ("checking a digits" & |
| Integer'Image (System.Max_Digits) & |
| " floating point type"); |
| end if; |
| |
| A_Long_Float_Check.Do_Test; |
| |
| |
| Report.Result; |
| end CXG2021; |
