|
111
|
1 -- CXG2010.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 exp function returns
|
|
|
28 -- results that are within the error bound allowed.
|
|
|
29 --
|
|
|
30 -- TEST DESCRIPTION:
|
|
|
31 -- This test contains three test packages that are almost
|
|
|
32 -- identical. The first two packages differ only in the
|
|
|
33 -- floating point type that is being tested. The first
|
|
|
34 -- and third package differ only in whether the generic
|
|
|
35 -- elementary functions package or the pre-instantiated
|
|
|
36 -- package is used.
|
|
|
37 -- The test package is not generic so that the arguments
|
|
|
38 -- and expected results for some of the test values
|
|
|
39 -- can be expressed as universal real instead of being
|
|
|
40 -- computed at runtime.
|
|
|
41 --
|
|
|
42 -- SPECIAL REQUIREMENTS
|
|
|
43 -- The Strict Mode for the numerical accuracy must be
|
|
|
44 -- selected. The method by which this mode is selected
|
|
|
45 -- is implementation dependent.
|
|
|
46 --
|
|
|
47 -- APPLICABILITY CRITERIA:
|
|
|
48 -- This test applies only to implementations supporting the
|
|
|
49 -- Numerics Annex and where the Machine_Radix is 2, 4, 8, or 16.
|
|
|
50 -- This test only applies to the Strict Mode for numerical
|
|
|
51 -- accuracy.
|
|
|
52 --
|
|
|
53 --
|
|
|
54 -- CHANGE HISTORY:
|
|
|
55 -- 1 Mar 96 SAIC Initial release for 2.1
|
|
|
56 -- 2 Sep 96 SAIC Improved check routine
|
|
|
57 --
|
|
|
58 --!
|
|
|
59
|
|
|
60 --
|
|
|
61 -- References:
|
|
|
62 --
|
|
|
63 -- Software Manual for the Elementary Functions
|
|
|
64 -- William J. Cody, Jr. and William Waite
|
|
|
65 -- Prentice-Hall, 1980
|
|
|
66 --
|
|
|
67 -- CRC Standard Mathematical Tables
|
|
|
68 -- 23rd Edition
|
|
|
69 --
|
|
|
70 -- Implementation and Testing of Function Software
|
|
|
71 -- W. J. Cody
|
|
|
72 -- Problems and Methodologies in Mathematical Software Production
|
|
|
73 -- editors P. C. Messina and A. Murli
|
|
|
74 -- Lecture Notes in Computer Science Volume 142
|
|
|
75 -- Springer Verlag, 1982
|
|
|
76 --
|
|
|
77
|
|
|
78 --
|
|
|
79 -- Notes on derivation of error bound for exp(p)*exp(-p)
|
|
|
80 --
|
|
|
81 -- Let a = true value of exp(p) and ac be the computed value.
|
|
|
82 -- Then a = ac(1+e1), where |e1| <= 4*Model_Epsilon.
|
|
|
83 -- Similarly, let b = true value of exp(-p) and bc be the computed value.
|
|
|
84 -- Then b = bc(1+e2), where |e2| <= 4*ME.
|
|
|
85 --
|
|
|
86 -- The product of x and y is (x*y)(1+e3), where |e3| <= 1.0ME
|
|
|
87 --
|
|
|
88 -- Hence, the computed ab is [ac(1+e1)*bc(1+e2)](1+e3) =
|
|
|
89 -- (ac*bc)[1 + e1 + e2 + e3 + e1e2 + e1e3 + e2e3 + e1e2e3).
|
|
|
90 --
|
|
|
91 -- Throwing away the last four tiny terms, we have (ac*bc)(1 + eta),
|
|
|
92 --
|
|
|
93 -- where |eta| <= (4+4+1)ME = 9.0Model_Epsilon.
|
|
|
94
|
|
|
95 with System;
|
|
|
96 with Report;
|
|
|
97 with Ada.Numerics.Generic_Elementary_Functions;
|
|
|
98 with Ada.Numerics.Elementary_Functions;
|
|
|
99 procedure CXG2010 is
|
|
|
100 Verbose : constant Boolean := False;
|
|
|
101 Max_Samples : constant := 1000;
|
|
|
102 Accuracy_Error_Reported : Boolean := False;
|
|
|
103
|
|
|
104 package Float_Check is
|
|
|
105 subtype Real is Float;
|
|
|
106 procedure Do_Test;
|
|
|
107 end Float_Check;
|
|
|
108
|
|
|
109 package body Float_Check is
|
|
|
110 package Elementary_Functions is new
|
|
|
111 Ada.Numerics.Generic_Elementary_Functions (Real);
|
|
|
112 function Sqrt (X : Real) return Real renames
|
|
|
113 Elementary_Functions.Sqrt;
|
|
|
114 function Exp (X : Real) return Real renames
|
|
|
115 Elementary_Functions.Exp;
|
|
|
116
|
|
|
117
|
|
|
118 -- The following value is a lower bound on the accuracy
|
|
|
119 -- required. It is normally 0.0 so that the lower bound
|
|
|
120 -- is computed from Model_Epsilon. However, for tests
|
|
|
121 -- where the expected result is only known to a certain
|
|
|
122 -- amount of precision this bound takes on a non-zero
|
|
|
123 -- value to account for that level of precision.
|
|
|
124 Error_Low_Bound : Real := 0.0;
|
|
|
125
|
|
|
126 procedure Check (Actual, Expected : Real;
|
|
|
127 Test_Name : String;
|
|
|
128 MRE : Real) is
|
|
|
129 Max_Error : Real;
|
|
|
130 Rel_Error : Real;
|
|
|
131 Abs_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
|
|
|
135 -- instead 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 -- take into account the low bound on the error
|
|
|
145 if Max_Error < Error_Low_Bound then
|
|
|
146 Max_Error := Error_Low_Bound;
|
|
|
147 end if;
|
|
|
148
|
|
|
149 if abs (Actual - Expected) > Max_Error then
|
|
|
150 Accuracy_Error_Reported := True;
|
|
|
151 Report.Failed (Test_Name &
|
|
|
152 " actual: " & Real'Image (Actual) &
|
|
|
153 " expected: " & Real'Image (Expected) &
|
|
|
154 " difference: " & Real'Image (Actual - Expected) &
|
|
|
155 " max err:" & Real'Image (Max_Error) );
|
|
|
156 elsif Verbose then
|
|
|
157 if Actual = Expected then
|
|
|
158 Report.Comment (Test_Name & " exact result");
|
|
|
159 else
|
|
|
160 Report.Comment (Test_Name & " passed");
|
|
|
161 end if;
|
|
|
162 end if;
|
|
|
163 end Check;
|
|
|
164
|
|
|
165
|
|
|
166 procedure Argument_Range_Check_1 (A, B : Real;
|
|
|
167 Test : String) is
|
|
|
168 -- test a evenly distributed selection of
|
|
|
169 -- arguments selected from the range A to B.
|
|
|
170 -- Test using identity: EXP(X-V) = EXP(X) * EXP (-V)
|
|
|
171 -- The parameter One_Minus_Exp_Minus_V is the value
|
|
|
172 -- 1.0 - Exp (-V)
|
|
|
173 -- accurate to machine precision.
|
|
|
174 -- This procedure is a translation of part of Cody's test
|
|
|
175 X : Real;
|
|
|
176 Y : Real;
|
|
|
177 ZX, ZY : Real;
|
|
|
178 V : constant := 1.0 / 16.0;
|
|
|
179 One_Minus_Exp_Minus_V : constant := 6.058693718652421388E-2;
|
|
|
180
|
|
|
181 begin
|
|
|
182 Accuracy_Error_Reported := False;
|
|
|
183 for I in 1..Max_Samples loop
|
|
|
184 X := (B - A) * Real (I) / Real (Max_Samples) + A;
|
|
|
185 Y := X - V;
|
|
|
186 if Y < 0.0 then
|
|
|
187 X := Y + V;
|
|
|
188 end if;
|
|
|
189
|
|
|
190 ZX := Exp (X);
|
|
|
191 ZY := Exp (Y);
|
|
|
192
|
|
|
193 -- ZX := Exp(X) - Exp(X) * (1 - Exp(-V);
|
|
|
194 -- which simplifies to ZX := Exp (X-V);
|
|
|
195 ZX := ZX - ZX * One_Minus_Exp_Minus_V;
|
|
|
196
|
|
|
197 -- note that since the expected value is computed, we
|
|
|
198 -- must take the error in that computation into account.
|
|
|
199 Check (ZY, ZX,
|
|
|
200 "test " & Test & " -" &
|
|
|
201 Integer'Image (I) &
|
|
|
202 " exp (" & Real'Image (X) & ")",
|
|
|
203 9.0);
|
|
|
204 exit when Accuracy_Error_Reported;
|
|
|
205 end loop;
|
|
|
206 exception
|
|
|
207 when Constraint_Error =>
|
|
|
208 Report.Failed
|
|
|
209 ("Constraint_Error raised in argument range check 1");
|
|
|
210 when others =>
|
|
|
211 Report.Failed ("exception in argument range check 1");
|
|
|
212 end Argument_Range_Check_1;
|
|
|
213
|
|
|
214
|
|
|
215
|
|
|
216 procedure Argument_Range_Check_2 (A, B : Real;
|
|
|
217 Test : String) is
|
|
|
218 -- test a evenly distributed selection of
|
|
|
219 -- arguments selected from the range A to B.
|
|
|
220 -- Test using identity: EXP(X-V) = EXP(X) * EXP (-V)
|
|
|
221 -- The parameter One_Minus_Exp_Minus_V is the value
|
|
|
222 -- 1.0 - Exp (-V)
|
|
|
223 -- accurate to machine precision.
|
|
|
224 -- This procedure is a translation of part of Cody's test
|
|
|
225 X : Real;
|
|
|
226 Y : Real;
|
|
|
227 ZX, ZY : Real;
|
|
|
228 V : constant := 45.0 / 16.0;
|
|
|
229 -- 1/16 - Exp(45/16)
|
|
|
230 Coeff : constant := 2.4453321046920570389E-3;
|
|
|
231
|
|
|
232 begin
|
|
|
233 Accuracy_Error_Reported := False;
|
|
|
234 for I in 1..Max_Samples loop
|
|
|
235 X := (B - A) * Real (I) / Real (Max_Samples) + A;
|
|
|
236 Y := X - V;
|
|
|
237 if Y < 0.0 then
|
|
|
238 X := Y + V;
|
|
|
239 end if;
|
|
|
240
|
|
|
241 ZX := Exp (X);
|
|
|
242 ZY := Exp (Y);
|
|
|
243
|
|
|
244 -- ZX := Exp(X) * 1/16 - Exp(X) * Coeff;
|
|
|
245 -- where Coeff is 1/16 - Exp(45/16)
|
|
|
246 -- which simplifies to ZX := Exp (X-V);
|
|
|
247 ZX := ZX * 0.0625 - ZX * Coeff;
|
|
|
248
|
|
|
249 -- note that since the expected value is computed, we
|
|
|
250 -- must take the error in that computation into account.
|
|
|
251 Check (ZY, ZX,
|
|
|
252 "test " & Test & " -" &
|
|
|
253 Integer'Image (I) &
|
|
|
254 " exp (" & Real'Image (X) & ")",
|
|
|
255 9.0);
|
|
|
256 exit when Accuracy_Error_Reported;
|
|
|
257 end loop;
|
|
|
258 exception
|
|
|
259 when Constraint_Error =>
|
|
|
260 Report.Failed
|
|
|
261 ("Constraint_Error raised in argument range check 2");
|
|
|
262 when others =>
|
|
|
263 Report.Failed ("exception in argument range check 2");
|
|
|
264 end Argument_Range_Check_2;
|
|
|
265
|
|
|
266
|
|
|
267 procedure Do_Test is
|
|
|
268 begin
|
|
|
269
|
|
|
270 --- test 1 ---
|
|
|
271 declare
|
|
|
272 Y : Real;
|
|
|
273 begin
|
|
|
274 Y := Exp(1.0);
|
|
|
275 -- normal accuracy requirements
|
|
|
276 Check (Y, Ada.Numerics.e, "test 1 -- exp(1)", 4.0);
|
|
|
277 exception
|
|
|
278 when Constraint_Error =>
|
|
|
279 Report.Failed ("Constraint_Error raised in test 1");
|
|
|
280 when others =>
|
|
|
281 Report.Failed ("exception in test 1");
|
|
|
282 end;
|
|
|
283
|
|
|
284 --- test 2 ---
|
|
|
285 declare
|
|
|
286 Y : Real;
|
|
|
287 begin
|
|
|
288 Y := Exp(16.0) * Exp(-16.0);
|
|
|
289 Check (Y, 1.0, "test 2 -- exp(16)*exp(-16)", 9.0);
|
|
|
290 exception
|
|
|
291 when Constraint_Error =>
|
|
|
292 Report.Failed ("Constraint_Error raised in test 2");
|
|
|
293 when others =>
|
|
|
294 Report.Failed ("exception in test 2");
|
|
|
295 end;
|
|
|
296
|
|
|
297 --- test 3 ---
|
|
|
298 declare
|
|
|
299 Y : Real;
|
|
|
300 begin
|
|
|
301 Y := Exp (Ada.Numerics.Pi) * Exp (-Ada.Numerics.Pi);
|
|
|
302 Check (Y, 1.0, "test 3 -- exp(pi)*exp(-pi)", 9.0);
|
|
|
303 exception
|
|
|
304 when Constraint_Error =>
|
|
|
305 Report.Failed ("Constraint_Error raised in test 3");
|
|
|
306 when others =>
|
|
|
307 Report.Failed ("exception in test 3");
|
|
|
308 end;
|
|
|
309
|
|
|
310 --- test 4 ---
|
|
|
311 declare
|
|
|
312 Y : Real;
|
|
|
313 begin
|
|
|
314 Y := Exp(0.0);
|
|
|
315 Check (Y, 1.0, "test 4 -- exp(0.0)",
|
|
|
316 0.0); -- no error allowed
|
|
|
317 exception
|
|
|
318 when Constraint_Error =>
|
|
|
319 Report.Failed ("Constraint_Error raised in test 4");
|
|
|
320 when others =>
|
|
|
321 Report.Failed ("exception in test 4");
|
|
|
322 end;
|
|
|
323
|
|
|
324 --- test 5 ---
|
|
|
325 -- constants used here only have 19 digits of precision
|
|
|
326 if Real'Digits > 19 then
|
|
|
327 Error_Low_Bound := 0.00000_00000_00000_0001;
|
|
|
328 Report.Comment ("exp accuracy checked to 19 digits");
|
|
|
329 end if;
|
|
|
330
|
|
|
331 Argument_Range_Check_1 ( 1.0/Sqrt(Real(Real'Machine_Radix)),
|
|
|
332 1.0,
|
|
|
333 "5");
|
|
|
334 Error_Low_Bound := 0.0; -- reset
|
|
|
335
|
|
|
336 --- test 6 ---
|
|
|
337 -- constants used here only have 19 digits of precision
|
|
|
338 if Real'Digits > 19 then
|
|
|
339 Error_Low_Bound := 0.00000_00000_00000_0001;
|
|
|
340 Report.Comment ("exp accuracy checked to 19 digits");
|
|
|
341 end if;
|
|
|
342
|
|
|
343 Argument_Range_Check_2 (1.0,
|
|
|
344 Sqrt(Real(Real'Machine_Radix)),
|
|
|
345 "6");
|
|
|
346 Error_Low_Bound := 0.0; -- reset
|
|
|
347
|
|
|
348 end Do_Test;
|
|
|
349 end Float_Check;
|
|
|
350
|
|
|
351 -----------------------------------------------------------------------
|
|
|
352 -----------------------------------------------------------------------
|
|
|
353 -- check the floating point type with the most digits
|
|
|
354 type A_Long_Float is digits System.Max_Digits;
|
|
|
355
|
|
|
356
|
|
|
357 package A_Long_Float_Check is
|
|
|
358 subtype Real is A_Long_Float;
|
|
|
359 procedure Do_Test;
|
|
|
360 end A_Long_Float_Check;
|
|
|
361
|
|
|
362 package body A_Long_Float_Check is
|
|
|
363 package Elementary_Functions is new
|
|
|
364 Ada.Numerics.Generic_Elementary_Functions (Real);
|
|
|
365 function Sqrt (X : Real) return Real renames
|
|
|
366 Elementary_Functions.Sqrt;
|
|
|
367 function Exp (X : Real) return Real renames
|
|
|
368 Elementary_Functions.Exp;
|
|
|
369
|
|
|
370
|
|
|
371 -- The following value is a lower bound on the accuracy
|
|
|
372 -- required. It is normally 0.0 so that the lower bound
|
|
|
373 -- is computed from Model_Epsilon. However, for tests
|
|
|
374 -- where the expected result is only known to a certain
|
|
|
375 -- amount of precision this bound takes on a non-zero
|
|
|
376 -- value to account for that level of precision.
|
|
|
377 Error_Low_Bound : Real := 0.0;
|
|
|
378
|
|
|
379 procedure Check (Actual, Expected : Real;
|
|
|
380 Test_Name : String;
|
|
|
381 MRE : Real) is
|
|
|
382 Max_Error : Real;
|
|
|
383 Rel_Error : Real;
|
|
|
384 Abs_Error : Real;
|
|
|
385 begin
|
|
|
386 -- In the case where the expected result is very small or 0
|
|
|
387 -- we compute the maximum error as a multiple of Model_Epsilon
|
|
|
388 -- instead of Model_Epsilon and Expected.
|
|
|
389 Rel_Error := MRE * abs Expected * Real'Model_Epsilon;
|
|
|
390 Abs_Error := MRE * Real'Model_Epsilon;
|
|
|
391 if Rel_Error > Abs_Error then
|
|
|
392 Max_Error := Rel_Error;
|
|
|
393 else
|
|
|
394 Max_Error := Abs_Error;
|
|
|
395 end if;
|
|
|
396
|
|
|
397 -- take into account the low bound on the error
|
|
|
398 if Max_Error < Error_Low_Bound then
|
|
|
399 Max_Error := Error_Low_Bound;
|
|
|
400 end if;
|
|
|
401
|
|
|
402 if abs (Actual - Expected) > Max_Error then
|
|
|
403 Accuracy_Error_Reported := True;
|
|
|
404 Report.Failed (Test_Name &
|
|
|
405 " actual: " & Real'Image (Actual) &
|
|
|
406 " expected: " & Real'Image (Expected) &
|
|
|
407 " difference: " & Real'Image (Actual - Expected) &
|
|
|
408 " max err:" & Real'Image (Max_Error) );
|
|
|
409 elsif Verbose then
|
|
|
410 if Actual = Expected then
|
|
|
411 Report.Comment (Test_Name & " exact result");
|
|
|
412 else
|
|
|
413 Report.Comment (Test_Name & " passed");
|
|
|
414 end if;
|
|
|
415 end if;
|
|
|
416 end Check;
|
|
|
417
|
|
|
418
|
|
|
419 procedure Argument_Range_Check_1 (A, B : Real;
|
|
|
420 Test : String) is
|
|
|
421 -- test a evenly distributed selection of
|
|
|
422 -- arguments selected from the range A to B.
|
|
|
423 -- Test using identity: EXP(X-V) = EXP(X) * EXP (-V)
|
|
|
424 -- The parameter One_Minus_Exp_Minus_V is the value
|
|
|
425 -- 1.0 - Exp (-V)
|
|
|
426 -- accurate to machine precision.
|
|
|
427 -- This procedure is a translation of part of Cody's test
|
|
|
428 X : Real;
|
|
|
429 Y : Real;
|
|
|
430 ZX, ZY : Real;
|
|
|
431 V : constant := 1.0 / 16.0;
|
|
|
432 One_Minus_Exp_Minus_V : constant := 6.058693718652421388E-2;
|
|
|
433
|
|
|
434 begin
|
|
|
435 Accuracy_Error_Reported := False;
|
|
|
436 for I in 1..Max_Samples loop
|
|
|
437 X := (B - A) * Real (I) / Real (Max_Samples) + A;
|
|
|
438 Y := X - V;
|
|
|
439 if Y < 0.0 then
|
|
|
440 X := Y + V;
|
|
|
441 end if;
|
|
|
442
|
|
|
443 ZX := Exp (X);
|
|
|
444 ZY := Exp (Y);
|
|
|
445
|
|
|
446 -- ZX := Exp(X) - Exp(X) * (1 - Exp(-V);
|
|
|
447 -- which simplifies to ZX := Exp (X-V);
|
|
|
448 ZX := ZX - ZX * One_Minus_Exp_Minus_V;
|
|
|
449
|
|
|
450 -- note that since the expected value is computed, we
|
|
|
451 -- must take the error in that computation into account.
|
|
|
452 Check (ZY, ZX,
|
|
|
453 "test " & Test & " -" &
|
|
|
454 Integer'Image (I) &
|
|
|
455 " exp (" & Real'Image (X) & ")",
|
|
|
456 9.0);
|
|
|
457 exit when Accuracy_Error_Reported;
|
|
|
458 end loop;
|
|
|
459 exception
|
|
|
460 when Constraint_Error =>
|
|
|
461 Report.Failed
|
|
|
462 ("Constraint_Error raised in argument range check 1");
|
|
|
463 when others =>
|
|
|
464 Report.Failed ("exception in argument range check 1");
|
|
|
465 end Argument_Range_Check_1;
|
|
|
466
|
|
|
467
|
|
|
468
|
|
|
469 procedure Argument_Range_Check_2 (A, B : Real;
|
|
|
470 Test : String) is
|
|
|
471 -- test a evenly distributed selection of
|
|
|
472 -- arguments selected from the range A to B.
|
|
|
473 -- Test using identity: EXP(X-V) = EXP(X) * EXP (-V)
|
|
|
474 -- The parameter One_Minus_Exp_Minus_V is the value
|
|
|
475 -- 1.0 - Exp (-V)
|
|
|
476 -- accurate to machine precision.
|
|
|
477 -- This procedure is a translation of part of Cody's test
|
|
|
478 X : Real;
|
|
|
479 Y : Real;
|
|
|
480 ZX, ZY : Real;
|
|
|
481 V : constant := 45.0 / 16.0;
|
|
|
482 -- 1/16 - Exp(45/16)
|
|
|
483 Coeff : constant := 2.4453321046920570389E-3;
|
|
|
484
|
|
|
485 begin
|
|
|
486 Accuracy_Error_Reported := False;
|
|
|
487 for I in 1..Max_Samples loop
|
|
|
488 X := (B - A) * Real (I) / Real (Max_Samples) + A;
|
|
|
489 Y := X - V;
|
|
|
490 if Y < 0.0 then
|
|
|
491 X := Y + V;
|
|
|
492 end if;
|
|
|
493
|
|
|
494 ZX := Exp (X);
|
|
|
495 ZY := Exp (Y);
|
|
|
496
|
|
|
497 -- ZX := Exp(X) * 1/16 - Exp(X) * Coeff;
|
|
|
498 -- where Coeff is 1/16 - Exp(45/16)
|
|
|
499 -- which simplifies to ZX := Exp (X-V);
|
|
|
500 ZX := ZX * 0.0625 - ZX * Coeff;
|
|
|
501
|
|
|
502 -- note that since the expected value is computed, we
|
|
|
503 -- must take the error in that computation into account.
|
|
|
504 Check (ZY, ZX,
|
|
|
505 "test " & Test & " -" &
|
|
|
506 Integer'Image (I) &
|
|
|
507 " exp (" & Real'Image (X) & ")",
|
|
|
508 9.0);
|
|
|
509 exit when Accuracy_Error_Reported;
|
|
|
510 end loop;
|
|
|
511 exception
|
|
|
512 when Constraint_Error =>
|
|
|
513 Report.Failed
|
|
|
514 ("Constraint_Error raised in argument range check 2");
|
|
|
515 when others =>
|
|
|
516 Report.Failed ("exception in argument range check 2");
|
|
|
517 end Argument_Range_Check_2;
|
|
|
518
|
|
|
519
|
|
|
520 procedure Do_Test is
|
|
|
521 begin
|
|
|
522
|
|
|
523 --- test 1 ---
|
|
|
524 declare
|
|
|
525 Y : Real;
|
|
|
526 begin
|
|
|
527 Y := Exp(1.0);
|
|
|
528 -- normal accuracy requirements
|
|
|
529 Check (Y, Ada.Numerics.e, "test 1 -- exp(1)", 4.0);
|
|
|
530 exception
|
|
|
531 when Constraint_Error =>
|
|
|
532 Report.Failed ("Constraint_Error raised in test 1");
|
|
|
533 when others =>
|
|
|
534 Report.Failed ("exception in test 1");
|
|
|
535 end;
|
|
|
536
|
|
|
537 --- test 2 ---
|
|
|
538 declare
|
|
|
539 Y : Real;
|
|
|
540 begin
|
|
|
541 Y := Exp(16.0) * Exp(-16.0);
|
|
|
542 Check (Y, 1.0, "test 2 -- exp(16)*exp(-16)", 9.0);
|
|
|
543 exception
|
|
|
544 when Constraint_Error =>
|
|
|
545 Report.Failed ("Constraint_Error raised in test 2");
|
|
|
546 when others =>
|
|
|
547 Report.Failed ("exception in test 2");
|
|
|
548 end;
|
|
|
549
|
|
|
550 --- test 3 ---
|
|
|
551 declare
|
|
|
552 Y : Real;
|
|
|
553 begin
|
|
|
554 Y := Exp (Ada.Numerics.Pi) * Exp (-Ada.Numerics.Pi);
|
|
|
555 Check (Y, 1.0, "test 3 -- exp(pi)*exp(-pi)", 9.0);
|
|
|
556 exception
|
|
|
557 when Constraint_Error =>
|
|
|
558 Report.Failed ("Constraint_Error raised in test 3");
|
|
|
559 when others =>
|
|
|
560 Report.Failed ("exception in test 3");
|
|
|
561 end;
|
|
|
562
|
|
|
563 --- test 4 ---
|
|
|
564 declare
|
|
|
565 Y : Real;
|
|
|
566 begin
|
|
|
567 Y := Exp(0.0);
|
|
|
568 Check (Y, 1.0, "test 4 -- exp(0.0)",
|
|
|
569 0.0); -- no error allowed
|
|
|
570 exception
|
|
|
571 when Constraint_Error =>
|
|
|
572 Report.Failed ("Constraint_Error raised in test 4");
|
|
|
573 when others =>
|
|
|
574 Report.Failed ("exception in test 4");
|
|
|
575 end;
|
|
|
576
|
|
|
577 --- test 5 ---
|
|
|
578 -- constants used here only have 19 digits of precision
|
|
|
579 if Real'Digits > 19 then
|
|
|
580 Error_Low_Bound := 0.00000_00000_00000_0001;
|
|
|
581 Report.Comment ("exp accuracy checked to 19 digits");
|
|
|
582 end if;
|
|
|
583
|
|
|
584 Argument_Range_Check_1 ( 1.0/Sqrt(Real(Real'Machine_Radix)),
|
|
|
585 1.0,
|
|
|
586 "5");
|
|
|
587 Error_Low_Bound := 0.0; -- reset
|
|
|
588
|
|
|
589 --- test 6 ---
|
|
|
590 -- constants used here only have 19 digits of precision
|
|
|
591 if Real'Digits > 19 then
|
|
|
592 Error_Low_Bound := 0.00000_00000_00000_0001;
|
|
|
593 Report.Comment ("exp accuracy checked to 19 digits");
|
|
|
594 end if;
|
|
|
595
|
|
|
596 Argument_Range_Check_2 (1.0,
|
|
|
597 Sqrt(Real(Real'Machine_Radix)),
|
|
|
598 "6");
|
|
|
599 Error_Low_Bound := 0.0; -- reset
|
|
|
600
|
|
|
601 end Do_Test;
|
|
|
602 end A_Long_Float_Check;
|
|
|
603
|
|
|
604 -----------------------------------------------------------------------
|
|
|
605 -----------------------------------------------------------------------
|
|
|
606
|
|
|
607 package Non_Generic_Check is
|
|
|
608 procedure Do_Test;
|
|
|
609 subtype Real is Float;
|
|
|
610 end Non_Generic_Check;
|
|
|
611
|
|
|
612 package body Non_Generic_Check is
|
|
|
613
|
|
|
614 package Elementary_Functions renames
|
|
|
615 Ada.Numerics.Elementary_Functions;
|
|
|
616 function Sqrt (X : Real) return Real renames
|
|
|
617 Elementary_Functions.Sqrt;
|
|
|
618 function Exp (X : Real) return Real renames
|
|
|
619 Elementary_Functions.Exp;
|
|
|
620
|
|
|
621
|
|
|
622 -- The following value is a lower bound on the accuracy
|
|
|
623 -- required. It is normally 0.0 so that the lower bound
|
|
|
624 -- is computed from Model_Epsilon. However, for tests
|
|
|
625 -- where the expected result is only known to a certain
|
|
|
626 -- amount of precision this bound takes on a non-zero
|
|
|
627 -- value to account for that level of precision.
|
|
|
628 Error_Low_Bound : Real := 0.0;
|
|
|
629
|
|
|
630 procedure Check (Actual, Expected : Real;
|
|
|
631 Test_Name : String;
|
|
|
632 MRE : Real) is
|
|
|
633 Max_Error : Real;
|
|
|
634 Rel_Error : Real;
|
|
|
635 Abs_Error : Real;
|
|
|
636 begin
|
|
|
637 -- In the case where the expected result is very small or 0
|
|
|
638 -- we compute the maximum error as a multiple of Model_Epsilon
|
|
|
639 -- instead of Model_Epsilon and Expected.
|
|
|
640 Rel_Error := MRE * abs Expected * Real'Model_Epsilon;
|
|
|
641 Abs_Error := MRE * Real'Model_Epsilon;
|
|
|
642 if Rel_Error > Abs_Error then
|
|
|
643 Max_Error := Rel_Error;
|
|
|
644 else
|
|
|
645 Max_Error := Abs_Error;
|
|
|
646 end if;
|
|
|
647
|
|
|
648 -- take into account the low bound on the error
|
|
|
649 if Max_Error < Error_Low_Bound then
|
|
|
650 Max_Error := Error_Low_Bound;
|
|
|
651 end if;
|
|
|
652
|
|
|
653 if abs (Actual - Expected) > Max_Error then
|
|
|
654 Accuracy_Error_Reported := True;
|
|
|
655 Report.Failed (Test_Name &
|
|
|
656 " actual: " & Real'Image (Actual) &
|
|
|
657 " expected: " & Real'Image (Expected) &
|
|
|
658 " difference: " & Real'Image (Actual - Expected) &
|
|
|
659 " max err:" & Real'Image (Max_Error) );
|
|
|
660 elsif Verbose then
|
|
|
661 if Actual = Expected then
|
|
|
662 Report.Comment (Test_Name & " exact result");
|
|
|
663 else
|
|
|
664 Report.Comment (Test_Name & " passed");
|
|
|
665 end if;
|
|
|
666 end if;
|
|
|
667 end Check;
|
|
|
668
|
|
|
669
|
|
|
670 procedure Argument_Range_Check_1 (A, B : Real;
|
|
|
671 Test : String) is
|
|
|
672 -- test a evenly distributed selection of
|
|
|
673 -- arguments selected from the range A to B.
|
|
|
674 -- Test using identity: EXP(X-V) = EXP(X) * EXP (-V)
|
|
|
675 -- The parameter One_Minus_Exp_Minus_V is the value
|
|
|
676 -- 1.0 - Exp (-V)
|
|
|
677 -- accurate to machine precision.
|
|
|
678 -- This procedure is a translation of part of Cody's test
|
|
|
679 X : Real;
|
|
|
680 Y : Real;
|
|
|
681 ZX, ZY : Real;
|
|
|
682 V : constant := 1.0 / 16.0;
|
|
|
683 One_Minus_Exp_Minus_V : constant := 6.058693718652421388E-2;
|
|
|
684
|
|
|
685 begin
|
|
|
686 Accuracy_Error_Reported := False;
|
|
|
687 for I in 1..Max_Samples loop
|
|
|
688 X := (B - A) * Real (I) / Real (Max_Samples) + A;
|
|
|
689 Y := X - V;
|
|
|
690 if Y < 0.0 then
|
|
|
691 X := Y + V;
|
|
|
692 end if;
|
|
|
693
|
|
|
694 ZX := Exp (X);
|
|
|
695 ZY := Exp (Y);
|
|
|
696
|
|
|
697 -- ZX := Exp(X) - Exp(X) * (1 - Exp(-V);
|
|
|
698 -- which simplifies to ZX := Exp (X-V);
|
|
|
699 ZX := ZX - ZX * One_Minus_Exp_Minus_V;
|
|
|
700
|
|
|
701 -- note that since the expected value is computed, we
|
|
|
702 -- must take the error in that computation into account.
|
|
|
703 Check (ZY, ZX,
|
|
|
704 "test " & Test & " -" &
|
|
|
705 Integer'Image (I) &
|
|
|
706 " exp (" & Real'Image (X) & ")",
|
|
|
707 9.0);
|
|
|
708 exit when Accuracy_Error_Reported;
|
|
|
709 end loop;
|
|
|
710 exception
|
|
|
711 when Constraint_Error =>
|
|
|
712 Report.Failed
|
|
|
713 ("Constraint_Error raised in argument range check 1");
|
|
|
714 when others =>
|
|
|
715 Report.Failed ("exception in argument range check 1");
|
|
|
716 end Argument_Range_Check_1;
|
|
|
717
|
|
|
718
|
|
|
719
|
|
|
720 procedure Argument_Range_Check_2 (A, B : Real;
|
|
|
721 Test : String) is
|
|
|
722 -- test a evenly distributed selection of
|
|
|
723 -- arguments selected from the range A to B.
|
|
|
724 -- Test using identity: EXP(X-V) = EXP(X) * EXP (-V)
|
|
|
725 -- The parameter One_Minus_Exp_Minus_V is the value
|
|
|
726 -- 1.0 - Exp (-V)
|
|
|
727 -- accurate to machine precision.
|
|
|
728 -- This procedure is a translation of part of Cody's test
|
|
|
729 X : Real;
|
|
|
730 Y : Real;
|
|
|
731 ZX, ZY : Real;
|
|
|
732 V : constant := 45.0 / 16.0;
|
|
|
733 -- 1/16 - Exp(45/16)
|
|
|
734 Coeff : constant := 2.4453321046920570389E-3;
|
|
|
735
|
|
|
736 begin
|
|
|
737 Accuracy_Error_Reported := False;
|
|
|
738 for I in 1..Max_Samples loop
|
|
|
739 X := (B - A) * Real (I) / Real (Max_Samples) + A;
|
|
|
740 Y := X - V;
|
|
|
741 if Y < 0.0 then
|
|
|
742 X := Y + V;
|
|
|
743 end if;
|
|
|
744
|
|
|
745 ZX := Exp (X);
|
|
|
746 ZY := Exp (Y);
|
|
|
747
|
|
|
748 -- ZX := Exp(X) * 1/16 - Exp(X) * Coeff;
|
|
|
749 -- where Coeff is 1/16 - Exp(45/16)
|
|
|
750 -- which simplifies to ZX := Exp (X-V);
|
|
|
751 ZX := ZX * 0.0625 - ZX * Coeff;
|
|
|
752
|
|
|
753 -- note that since the expected value is computed, we
|
|
|
754 -- must take the error in that computation into account.
|
|
|
755 Check (ZY, ZX,
|
|
|
756 "test " & Test & " -" &
|
|
|
757 Integer'Image (I) &
|
|
|
758 " exp (" & Real'Image (X) & ")",
|
|
|
759 9.0);
|
|
|
760 exit when Accuracy_Error_Reported;
|
|
|
761 end loop;
|
|
|
762 exception
|
|
|
763 when Constraint_Error =>
|
|
|
764 Report.Failed
|
|
|
765 ("Constraint_Error raised in argument range check 2");
|
|
|
766 when others =>
|
|
|
767 Report.Failed ("exception in argument range check 2");
|
|
|
768 end Argument_Range_Check_2;
|
|
|
769
|
|
|
770
|
|
|
771 procedure Do_Test is
|
|
|
772 begin
|
|
|
773
|
|
|
774 --- test 1 ---
|
|
|
775 declare
|
|
|
776 Y : Real;
|
|
|
777 begin
|
|
|
778 Y := Exp(1.0);
|
|
|
779 -- normal accuracy requirements
|
|
|
780 Check (Y, Ada.Numerics.e, "test 1 -- exp(1)", 4.0);
|
|
|
781 exception
|
|
|
782 when Constraint_Error =>
|
|
|
783 Report.Failed ("Constraint_Error raised in test 1");
|
|
|
784 when others =>
|
|
|
785 Report.Failed ("exception in test 1");
|
|
|
786 end;
|
|
|
787
|
|
|
788 --- test 2 ---
|
|
|
789 declare
|
|
|
790 Y : Real;
|
|
|
791 begin
|
|
|
792 Y := Exp(16.0) * Exp(-16.0);
|
|
|
793 Check (Y, 1.0, "test 2 -- exp(16)*exp(-16)", 9.0);
|
|
|
794 exception
|
|
|
795 when Constraint_Error =>
|
|
|
796 Report.Failed ("Constraint_Error raised in test 2");
|
|
|
797 when others =>
|
|
|
798 Report.Failed ("exception in test 2");
|
|
|
799 end;
|
|
|
800
|
|
|
801 --- test 3 ---
|
|
|
802 declare
|
|
|
803 Y : Real;
|
|
|
804 begin
|
|
|
805 Y := Exp (Ada.Numerics.Pi) * Exp (-Ada.Numerics.Pi);
|
|
|
806 Check (Y, 1.0, "test 3 -- exp(pi)*exp(-pi)", 9.0);
|
|
|
807 exception
|
|
|
808 when Constraint_Error =>
|
|
|
809 Report.Failed ("Constraint_Error raised in test 3");
|
|
|
810 when others =>
|
|
|
811 Report.Failed ("exception in test 3");
|
|
|
812 end;
|
|
|
813
|
|
|
814 --- test 4 ---
|
|
|
815 declare
|
|
|
816 Y : Real;
|
|
|
817 begin
|
|
|
818 Y := Exp(0.0);
|
|
|
819 Check (Y, 1.0, "test 4 -- exp(0.0)",
|
|
|
820 0.0); -- no error allowed
|
|
|
821 exception
|
|
|
822 when Constraint_Error =>
|
|
|
823 Report.Failed ("Constraint_Error raised in test 4");
|
|
|
824 when others =>
|
|
|
825 Report.Failed ("exception in test 4");
|
|
|
826 end;
|
|
|
827
|
|
|
828 --- test 5 ---
|
|
|
829 -- constants used here only have 19 digits of precision
|
|
|
830 if Real'Digits > 19 then
|
|
|
831 Error_Low_Bound := 0.00000_00000_00000_0001;
|
|
|
832 Report.Comment ("exp accuracy checked to 19 digits");
|
|
|
833 end if;
|
|
|
834
|
|
|
835 Argument_Range_Check_1 ( 1.0/Sqrt(Real(Real'Machine_Radix)),
|
|
|
836 1.0,
|
|
|
837 "5");
|
|
|
838 Error_Low_Bound := 0.0; -- reset
|
|
|
839
|
|
|
840 --- test 6 ---
|
|
|
841 -- constants used here only have 19 digits of precision
|
|
|
842 if Real'Digits > 19 then
|
|
|
843 Error_Low_Bound := 0.00000_00000_00000_0001;
|
|
|
844 Report.Comment ("exp accuracy checked to 19 digits");
|
|
|
845 end if;
|
|
|
846
|
|
|
847 Argument_Range_Check_2 (1.0,
|
|
|
848 Sqrt(Real(Real'Machine_Radix)),
|
|
|
849 "6");
|
|
|
850 Error_Low_Bound := 0.0; -- reset
|
|
|
851
|
|
|
852 end Do_Test;
|
|
|
853 end Non_Generic_Check;
|
|
|
854
|
|
|
855 -----------------------------------------------------------------------
|
|
|
856 -----------------------------------------------------------------------
|
|
|
857
|
|
|
858 begin
|
|
|
859 Report.Test ("CXG2010",
|
|
|
860 "Check the accuracy of the exp function");
|
|
|
861
|
|
|
862 -- the test only applies to machines with a radix of 2,4,8, or 16
|
|
|
863 case Float'Machine_Radix is
|
|
|
864 when 2 | 4 | 8 | 16 => null;
|
|
|
865 when others =>
|
|
|
866 Report.Not_Applicable ("only applicable to binary radix");
|
|
|
867 Report.Result;
|
|
|
868 return;
|
|
|
869 end case;
|
|
|
870
|
|
|
871 if Verbose then
|
|
|
872 Report.Comment ("checking Standard.Float");
|
|
|
873 end if;
|
|
|
874
|
|
|
875 Float_Check.Do_Test;
|
|
|
876
|
|
|
877 if Verbose then
|
|
|
878 Report.Comment ("checking a digits" &
|
|
|
879 Integer'Image (System.Max_Digits) &
|
|
|
880 " floating point type");
|
|
|
881 end if;
|
|
|
882
|
|
|
883 A_Long_Float_Check.Do_Test;
|
|
|
884
|
|
|
885 if Verbose then
|
|
|
886 Report.Comment ("checking non-generic package");
|
|
|
887 end if;
|
|
|
888
|
|
|
889 Non_Generic_Check.Do_Test;
|
|
|
890
|
|
|
891 Report.Result;
|
|
|
892 end CXG2010;
|
