|
111
|
1 -- CXG2021.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 complex SIN and COS functions return
|
|
|
28 -- a result that is within the error bound allowed.
|
|
|
29 --
|
|
|
30 -- TEST DESCRIPTION:
|
|
|
31 -- This test consists of a generic package that is
|
|
|
32 -- instantiated to check complex numbers based upon
|
|
|
33 -- both Float and a long float type.
|
|
|
34 -- The test for each floating point type is divided into
|
|
|
35 -- several parts:
|
|
|
36 -- Special value checks where the result is a known constant.
|
|
|
37 -- Checks that use an identity for determining the result.
|
|
|
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 -- 27 Mar 96 SAIC Initial release for 2.1
|
|
|
53 -- 22 Aug 96 SAIC No longer skips test for systems with
|
|
|
54 -- more than 20 digits of precision.
|
|
|
55 --
|
|
|
56 --!
|
|
|
57
|
|
|
58 --
|
|
|
59 -- References:
|
|
|
60 --
|
|
|
61 -- W. J. Cody
|
|
|
62 -- CELEFUNT: A Portable Test Package for Complex Elementary Functions
|
|
|
63 -- Algorithm 714, Collected Algorithms from ACM.
|
|
|
64 -- Published in Transactions On Mathematical Software,
|
|
|
65 -- Vol. 19, No. 1, March, 1993, pp. 1-21.
|
|
|
66 --
|
|
|
67 -- CRC Standard Mathematical Tables
|
|
|
68 -- 23rd Edition
|
|
|
69 --
|
|
|
70
|
|
|
71 with System;
|
|
|
72 with Report;
|
|
|
73 with Ada.Numerics.Generic_Complex_Types;
|
|
|
74 with Ada.Numerics.Generic_Complex_Elementary_Functions;
|
|
|
75 procedure CXG2021 is
|
|
|
76 Verbose : constant Boolean := False;
|
|
|
77 -- Note that Max_Samples is the number of samples taken in
|
|
|
78 -- both the real and imaginary directions. Thus, for Max_Samples
|
|
|
79 -- of 100 the number of values checked is 10000.
|
|
|
80 Max_Samples : constant := 100;
|
|
|
81
|
|
|
82 E : constant := Ada.Numerics.E;
|
|
|
83 Pi : constant := Ada.Numerics.Pi;
|
|
|
84
|
|
|
85 generic
|
|
|
86 type Real is digits <>;
|
|
|
87 package Generic_Check is
|
|
|
88 procedure Do_Test;
|
|
|
89 end Generic_Check;
|
|
|
90
|
|
|
91 package body Generic_Check is
|
|
|
92 package Complex_Type is new
|
|
|
93 Ada.Numerics.Generic_Complex_Types (Real);
|
|
|
94 use Complex_Type;
|
|
|
95
|
|
|
96 package CEF is new
|
|
|
97 Ada.Numerics.Generic_Complex_Elementary_Functions (Complex_Type);
|
|
|
98
|
|
|
99 function Sin (X : Complex) return Complex renames CEF.Sin;
|
|
|
100 function Cos (X : Complex) return Complex renames CEF.Cos;
|
|
|
101
|
|
|
102 -- flag used to terminate some tests early
|
|
|
103 Accuracy_Error_Reported : Boolean := False;
|
|
|
104
|
|
|
105 -- The following value is a lower bound on the accuracy
|
|
|
106 -- required. It is normally 0.0 so that the lower bound
|
|
|
107 -- is computed from Model_Epsilon. However, for tests
|
|
|
108 -- where the expected result is only known to a certain
|
|
|
109 -- amount of precision this bound takes on a non-zero
|
|
|
110 -- value to account for that level of precision.
|
|
|
111 Error_Low_Bound : Real := 0.0;
|
|
|
112
|
|
|
113 -- the E_Factor is an additional amount added to the Expected
|
|
|
114 -- value prior to computing the maximum relative error.
|
|
|
115 -- This is needed because the error analysis (Cody pg 17-20)
|
|
|
116 -- requires this additional allowance.
|
|
|
117 procedure Check (Actual, Expected : Real;
|
|
|
118 Test_Name : String;
|
|
|
119 MRE : Real;
|
|
|
120 E_Factor : Real := 0.0) is
|
|
|
121 Max_Error : Real;
|
|
|
122 Rel_Error : Real;
|
|
|
123 Abs_Error : Real;
|
|
|
124 begin
|
|
|
125 -- In the case where the expected result is very small or 0
|
|
|
126 -- we compute the maximum error as a multiple of Model_Epsilon instead
|
|
|
127 -- of Model_Epsilon and Expected.
|
|
|
128 Rel_Error := MRE * Real'Model_Epsilon * (abs Expected + E_Factor);
|
|
|
129 Abs_Error := MRE * Real'Model_Epsilon;
|
|
|
130 if Rel_Error > Abs_Error then
|
|
|
131 Max_Error := Rel_Error;
|
|
|
132 else
|
|
|
133 Max_Error := Abs_Error;
|
|
|
134 end if;
|
|
|
135
|
|
|
136 -- take into account the low bound on the error
|
|
|
137 if Max_Error < Error_Low_Bound then
|
|
|
138 Max_Error := Error_Low_Bound;
|
|
|
139 end if;
|
|
|
140
|
|
|
141 if abs (Actual - Expected) > Max_Error then
|
|
|
142 Accuracy_Error_Reported := True;
|
|
|
143 Report.Failed (Test_Name &
|
|
|
144 " actual: " & Real'Image (Actual) &
|
|
|
145 " expected: " & Real'Image (Expected) &
|
|
|
146 " difference: " & Real'Image (Actual - Expected) &
|
|
|
147 " max err:" & Real'Image (Max_Error) &
|
|
|
148 " efactor:" & Real'Image (E_Factor) );
|
|
|
149 elsif Verbose then
|
|
|
150 if Actual = Expected then
|
|
|
151 Report.Comment (Test_Name & " exact result");
|
|
|
152 else
|
|
|
153 Report.Comment (Test_Name & " passed" &
|
|
|
154 " actual: " & Real'Image (Actual) &
|
|
|
155 " expected: " & Real'Image (Expected) &
|
|
|
156 " difference: " & Real'Image (Actual - Expected) &
|
|
|
157 " max err:" & Real'Image (Max_Error) &
|
|
|
158 " efactor:" & Real'Image (E_Factor) );
|
|
|
159 end if;
|
|
|
160 end if;
|
|
|
161 end Check;
|
|
|
162
|
|
|
163
|
|
|
164 procedure Check (Actual, Expected : Complex;
|
|
|
165 Test_Name : String;
|
|
|
166 MRE : Real;
|
|
|
167 R_Factor, I_Factor : Real := 0.0) is
|
|
|
168 begin
|
|
|
169 Check (Actual.Re, Expected.Re, Test_Name & " real part",
|
|
|
170 MRE, R_Factor);
|
|
|
171 Check (Actual.Im, Expected.Im, Test_Name & " imaginary part",
|
|
|
172 MRE, I_Factor);
|
|
|
173 end Check;
|
|
|
174
|
|
|
175
|
|
|
176 procedure Special_Value_Test is
|
|
|
177 -- In the following tests the expected result is accurate
|
|
|
178 -- to the machine precision so the minimum guaranteed error
|
|
|
179 -- bound can be used if the argument is exact.
|
|
|
180 -- Since the argument involves Pi, we must allow for this
|
|
|
181 -- inexact argument.
|
|
|
182 Minimum_Error : constant := 11.0;
|
|
|
183 begin
|
|
|
184 Check (Sin (Pi/2.0 + 0.0*i),
|
|
|
185 1.0 + 0.0*i,
|
|
|
186 "sin(pi/2+0i)",
|
|
|
187 Minimum_Error + 1.0);
|
|
|
188 Check (Cos (Pi/2.0 + 0.0*i),
|
|
|
189 0.0 + 0.0*i,
|
|
|
190 "cos(pi/2+0i)",
|
|
|
191 Minimum_Error + 1.0);
|
|
|
192 exception
|
|
|
193 when Constraint_Error =>
|
|
|
194 Report.Failed ("Constraint_Error raised in special value test");
|
|
|
195 when others =>
|
|
|
196 Report.Failed ("exception in special value test");
|
|
|
197 end Special_Value_Test;
|
|
|
198
|
|
|
199
|
|
|
200
|
|
|
201 procedure Exact_Result_Test is
|
|
|
202 No_Error : constant := 0.0;
|
|
|
203 begin
|
|
|
204 -- G.1.2(36);6.0
|
|
|
205 Check (Sin(0.0 + 0.0*i), 0.0 + 0.0 * i, "sin(0+0i)", No_Error);
|
|
|
206 Check (Cos(0.0 + 0.0*i), 1.0 + 0.0 * i, "cos(0+0i)", No_Error);
|
|
|
207 exception
|
|
|
208 when Constraint_Error =>
|
|
|
209 Report.Failed ("Constraint_Error raised in Exact_Result Test");
|
|
|
210 when others =>
|
|
|
211 Report.Failed ("exception in Exact_Result Test");
|
|
|
212 end Exact_Result_Test;
|
|
|
213
|
|
|
214
|
|
|
215 procedure Identity_Test (RA, RB, IA, IB : Real) is
|
|
|
216 -- Tests an identity over a range of values specified
|
|
|
217 -- by the 4 parameters. RA and RB denote the range for the
|
|
|
218 -- real part while IA and IB denote the range for the
|
|
|
219 -- imaginary part.
|
|
|
220 --
|
|
|
221 -- For this test we use the identity
|
|
|
222 -- Sin(Z) = Sin(Z-W) * Cos(W) + Cos(Z-W) * Sin(W)
|
|
|
223 -- and
|
|
|
224 -- Cos(Z) = Cos(Z-W) * Cos(W) - Sin(Z-W) * Sin(W)
|
|
|
225 --
|
|
|
226
|
|
|
227 X, Y : Real;
|
|
|
228 Z : Complex;
|
|
|
229 W : constant Complex := Compose_From_Cartesian(0.0625, 0.0625);
|
|
|
230 ZmW : Complex; -- Z - W
|
|
|
231 Sin_ZmW,
|
|
|
232 Cos_ZmW : Complex;
|
|
|
233 Actual1, Actual2 : Complex;
|
|
|
234 R_Factor : Real; -- additional real error factor
|
|
|
235 I_Factor : Real; -- additional imaginary error factor
|
|
|
236 Sin_W : constant Complex := (6.2581348413276935585E-2,
|
|
|
237 6.2418588008436587236E-2);
|
|
|
238 -- numeric stability is enhanced by using Cos(W) - 1.0 instead of
|
|
|
239 -- Cos(W) in the computation.
|
|
|
240 Cos_W_m_1 : constant Complex := (-2.5431314180235545803E-6,
|
|
|
241 -3.9062493377261771826E-3);
|
|
|
242
|
|
|
243
|
|
|
244 begin
|
|
|
245 if Real'Digits > 20 then
|
|
|
246 -- constants used here accurate to 20 digits. Allow 1
|
|
|
247 -- additional digit of error for computation.
|
|
|
248 Error_Low_Bound := 0.00000_00000_00000_0001;
|
|
|
249 Report.Comment ("accuracy checked to 19 digits");
|
|
|
250 end if;
|
|
|
251
|
|
|
252 Accuracy_Error_Reported := False; -- reset
|
|
|
253 for II in 0..Max_Samples loop
|
|
|
254 X := (RB - RA) * Real (II) / Real (Max_Samples) + RA;
|
|
|
255 for J in 0..Max_Samples loop
|
|
|
256 Y := (IB - IA) * Real (J) / Real (Max_Samples) + IA;
|
|
|
257
|
|
|
258 Z := Compose_From_Cartesian(X,Y);
|
|
|
259 ZmW := Z - W;
|
|
|
260 Sin_ZmW := Sin (ZmW);
|
|
|
261 Cos_ZmW := Cos (ZmW);
|
|
|
262
|
|
|
263 -- now for the first identity
|
|
|
264 -- Sin(Z) = Sin(Z-W) * Cos(W) + Cos(Z-W) * Sin(W)
|
|
|
265 -- = Sin(Z-W) * (1+(Cos(W)-1)) + Cos(Z-W) * Sin(W)
|
|
|
266 -- = Sin(Z-W) + Sin(Z-W)*(Cos(W)-1) + Cos(Z-W)*Sin(W)
|
|
|
267
|
|
|
268
|
|
|
269 Actual1 := Sin (Z);
|
|
|
270 Actual2 := Sin_ZmW + (Sin_ZmW * Cos_W_m_1 + Cos_ZmW * Sin_W);
|
|
|
271
|
|
|
272 -- The computation of the additional error factors are taken
|
|
|
273 -- from Cody pages 17-20.
|
|
|
274
|
|
|
275 R_Factor := abs (Re (Sin_ZmW) * Re (1.0 - Cos_W_m_1)) +
|
|
|
276 abs (Im (Sin_ZmW) * Im (1.0 - Cos_W_m_1)) +
|
|
|
277 abs (Re (Cos_ZmW) * Re (Sin_W)) +
|
|
|
278 abs (Re (Cos_ZmW) * Re (1.0 - Cos_W_m_1));
|
|
|
279
|
|
|
280 I_Factor := abs (Re (Sin_ZmW) * Im (1.0 - Cos_W_m_1)) +
|
|
|
281 abs (Im (Sin_ZmW) * Re (1.0 - Cos_W_m_1)) +
|
|
|
282 abs (Re (Cos_ZmW) * Im (Sin_W)) +
|
|
|
283 abs (Im (Cos_ZmW) * Re (1.0 - Cos_W_m_1));
|
|
|
284
|
|
|
285 Check (Actual1, Actual2,
|
|
|
286 "Identity_1_Test " & Integer'Image (II) &
|
|
|
287 Integer'Image (J) & ": Sin((" &
|
|
|
288 Real'Image (Z.Re) & ", " &
|
|
|
289 Real'Image (Z.Im) & ")) ",
|
|
|
290 11.0, R_Factor, I_Factor);
|
|
|
291
|
|
|
292 -- now for the second identity
|
|
|
293 -- Cos(Z) = Cos(Z-W) * Cos(W) - Sin(Z-W) * Sin(W)
|
|
|
294 -- = Cos(Z-W) * (1+(Cos(W)-1) - Sin(Z-W) * Sin(W)
|
|
|
295 Actual1 := Cos (Z);
|
|
|
296 Actual2 := Cos_ZmW + (Cos_ZmW * Cos_W_m_1 - Sin_ZmW * Sin_W);
|
|
|
297
|
|
|
298 -- The computation of the additional error factors are taken
|
|
|
299 -- from Cody pages 17-20.
|
|
|
300
|
|
|
301 R_Factor := abs (Re (Sin_ZmW) * Re (Sin_W)) +
|
|
|
302 abs (Im (Sin_ZmW) * Im (Sin_W)) +
|
|
|
303 abs (Re (Cos_ZmW) * Re (1.0 - Cos_W_m_1)) +
|
|
|
304 abs (Im (Cos_ZmW) * Im (1.0 - Cos_W_m_1));
|
|
|
305
|
|
|
306 I_Factor := abs (Re (Sin_ZmW) * Im (Sin_W)) +
|
|
|
307 abs (Im (Sin_ZmW) * Re (Sin_W)) +
|
|
|
308 abs (Re (Cos_ZmW) * Im (1.0 - Cos_W_m_1)) +
|
|
|
309 abs (Im (Cos_ZmW) * Re (1.0 - Cos_W_m_1));
|
|
|
310
|
|
|
311 Check (Actual1, Actual2,
|
|
|
312 "Identity_2_Test " & Integer'Image (II) &
|
|
|
313 Integer'Image (J) & ": Cos((" &
|
|
|
314 Real'Image (Z.Re) & ", " &
|
|
|
315 Real'Image (Z.Im) & ")) ",
|
|
|
316 11.0, R_Factor, I_Factor);
|
|
|
317
|
|
|
318 if Accuracy_Error_Reported then
|
|
|
319 -- only report the first error in this test in order to keep
|
|
|
320 -- lots of failures from producing a huge error log
|
|
|
321 Error_Low_Bound := 0.0; -- reset
|
|
|
322 return;
|
|
|
323 end if;
|
|
|
324 end loop;
|
|
|
325 end loop;
|
|
|
326
|
|
|
327 Error_Low_Bound := 0.0; -- reset
|
|
|
328 exception
|
|
|
329 when Constraint_Error =>
|
|
|
330 Report.Failed
|
|
|
331 ("Constraint_Error raised in Identity_Test" &
|
|
|
332 " for Z=(" & Real'Image (X) &
|
|
|
333 ", " & Real'Image (Y) & ")");
|
|
|
334 when others =>
|
|
|
335 Report.Failed ("exception in Identity_Test" &
|
|
|
336 " for Z=(" & Real'Image (X) &
|
|
|
337 ", " & Real'Image (Y) & ")");
|
|
|
338 end Identity_Test;
|
|
|
339
|
|
|
340
|
|
|
341 procedure Do_Test is
|
|
|
342 begin
|
|
|
343 Special_Value_Test;
|
|
|
344 Exact_Result_Test;
|
|
|
345 -- test regions where sin and cos have the same sign and
|
|
|
346 -- about the same magnitude. This will minimize subtraction
|
|
|
347 -- errors in the identities.
|
|
|
348 -- See Cody page 17.
|
|
|
349 Identity_Test (0.0625, 10.0, 0.0625, 10.0);
|
|
|
350 Identity_Test ( 16.0, 17.0, 16.0, 17.0);
|
|
|
351 end Do_Test;
|
|
|
352 end Generic_Check;
|
|
|
353
|
|
|
354 -----------------------------------------------------------------------
|
|
|
355 -----------------------------------------------------------------------
|
|
|
356 package Float_Check is new Generic_Check (Float);
|
|
|
357
|
|
|
358 -- check the floating point type with the most digits
|
|
|
359 type A_Long_Float is digits System.Max_Digits;
|
|
|
360 package A_Long_Float_Check is new Generic_Check (A_Long_Float);
|
|
|
361
|
|
|
362 -----------------------------------------------------------------------
|
|
|
363 -----------------------------------------------------------------------
|
|
|
364
|
|
|
365
|
|
|
366 begin
|
|
|
367 Report.Test ("CXG2021",
|
|
|
368 "Check the accuracy of the complex SIN and COS functions");
|
|
|
369
|
|
|
370 if Verbose then
|
|
|
371 Report.Comment ("checking Standard.Float");
|
|
|
372 end if;
|
|
|
373
|
|
|
374 Float_Check.Do_Test;
|
|
|
375
|
|
|
376 if Verbose then
|
|
|
377 Report.Comment ("checking a digits" &
|
|
|
378 Integer'Image (System.Max_Digits) &
|
|
|
379 " floating point type");
|
|
|
380 end if;
|
|
|
381
|
|
|
382 A_Long_Float_Check.Do_Test;
|
|
|
383
|
|
|
384
|
|
|
385 Report.Result;
|
|
|
386 end CXG2021;
|
