-- ----------------------------------------------------------------------- -- Title: test_relerr_2 -- Last Mod: Mon Oct 1 17:37:23 1990 -- Author: Vincent Broman -- Visibility: instantiated and used by a test driver -- Description: -- A generic procedure for testing relative accuracy of -- floating point functions of two arguments. -- Random pairs of arguments are chosen from two intervals and -- the results from the candidate function and a standard -- function are compared. Overflow is treated as in the -- generic_elementary_functions standard, except that -- if the standard function just barely overflows, and -- the candidate produces a finite result satisfying its -- laxer error criterion we assign an error of 0.0 to the trial. -- This may underestimate the overall error only if the -- long_real'safe_large <= real'safe_large. -- -- The results are returned in the form of the maximum and -- average magnitude of the relative error, divided by its error -- bound, as well as the argument which caused the worst error -- (useful for diagnosis). -- Exceptions: none -- ----------------------------------------------------------------------- with random_number, generic_primitive_functions; generic type real is digits <>; type long_real is digits <>; -- candidate for testing with function cand( x, y: in real) return real; -- standard for comparison with function std( x, y: in long_real) return long_real; -- required rel acc eps: long_real; -- function of the test std fn result returning a factor (defaulting to one) -- which increases the allowed error bound. with function acc_factor( res: long_real) return long_real; -- want exponential distribution on random args tested? exp_dist_arg1, exp_dist_arg2: boolean; procedure test_relerr_2( niter: in positive; -- nbr of tests in run minarg1, maxarg1: in real; -- bounds on domain minarg2, maxarg2: in real; -- bounds on domain avg_relerr: out long_real; -- L1 err norm / niter max_relerr: out long_real; -- Linf err norm worst_arg1, worst_arg2: out long_real); -- arg realizing max err procedure test_relerr_2( niter: in positive; -- nbr of tests in run minarg1, maxarg1: in real; -- bounds on domain minarg2, maxarg2: in real; -- bounds on domain avg_relerr: out long_real; -- L1 err norm / niter max_relerr: out long_real; -- Linf err norm worst_arg1, worst_arg2: out long_real) is -- arg realizing max err -- max_err_arg1: long_real := 0.0; max_err_arg2: long_real := 0.0; max_err: long_real := 0.0; sum_err: long_real := 0.0; relerr: long_real; arg1: real; arg2: real; -- following used only if exp_dist_arg1 true aminarg1: real; amaxarg1: real; negate1: boolean; lowexp1: integer; exprg1: real; -- following used only if exp_dist_arg2 true aminarg2: real; amaxarg2: real; negate2: boolean; lowexp2: integer; exprg2: real; package primf is new generic_primitive_functions( real, integer); use primf; procedure domain_scaling is -- begin if exp_dist_arg1 then aminarg1 := abs( minarg1); amaxarg1 := abs( maxarg1); negate1 := maxarg1 < 0.0; if aminarg1 = 0.0 then aminarg1 := real'small; end if; lowexp1 := exponent( aminarg1) - 1; exprg1 := real( exponent( amaxarg1) - lowexp1); end if; if exp_dist_arg2 then aminarg2 := abs( minarg2); amaxarg2 := abs( maxarg2); negate2 := maxarg2 < 0.0; if aminarg2 = 0.0 then aminarg2 := real'small; end if; lowexp2 := exponent( aminarg2) - 1; exprg2 := real( exponent( amaxarg2) - lowexp2); end if; end domain_scaling; function get_rand_arg1 return real is -- arg: real; mant: real; exd: integer; begin -- get_rand_arg1 loop begin if exp_dist_arg1 then mant := real( random_number.next) + 1.0; exd := integer( floor( exprg1 * real( random_number.next))); arg := scale( mant, lowexp1 + exd); if negate1 then arg := -arg; end if; else arg := minarg1 + (maxarg1 - minarg1) * real( random_number.next); end if; -- just perturb for variety arg := arg * (0.78 + 0.5 * real( random_number.next)); -- trim slop in the distribution, quit loop when satis if exp_dist_arg1 then exit when aminarg1 <= abs( arg) and abs( arg) <= amaxarg1; else exit when minarg1 <= arg and arg <= maxarg1; end if; exception when numeric_error| constraint_error=> null; end; end loop; return arg; end get_rand_arg1; function get_rand_arg2 return real is -- arg: real; mant: real; exd: integer; begin -- get_rand_arg2 loop begin if exp_dist_arg2 then mant := real( random_number.next) + 1.0; exd := integer( floor( exprg2 * real( random_number.next))); arg := scale( mant, lowexp2 + exd); if negate2 then arg := -arg; end if; else arg := minarg2 + (maxarg2 - minarg2) * real( random_number.next); end if; -- just perturb for variety arg := arg * (0.78 + 0.5 * real( random_number.next)); -- trim slop in the distribution, quit loop when satis if exp_dist_arg2 then exit when aminarg2 <= abs( arg) and abs( arg) <= amaxarg2; else exit when minarg2 <= arg and arg <= maxarg2; end if; exception when numeric_error| constraint_error=> null; end; end loop; return arg; end get_rand_arg2; function get_relerr( x, y: real) return long_real is -- -- Compute relative error of cand compared to std. -- Cannot handle the rare case where std barely overflows, -- but cand returns a number large enough to meet its laxer error criterion. -- ca, st: long_real; ca_ovfl: boolean := false; st_ovfl: boolean := false; big_real: long_real := long_real( real'safe_large * (1.0 - eps)); small_real: long_real := long_real( real'safe_small); total_miss: long_real := long_real'safe_large / long_real( niter); function epsi return long_real is begin return eps * acc_factor( st); end epsi; begin -- get_relerr begin ca := long_real( cand( x, y)); exception when numeric_error| constraint_error=> ca_ovfl := true; -- when argument_error we have a bug end; begin st := std( long_real( x), long_real( y)); exception when numeric_error| constraint_error=> st_ovfl := true; -- when argument_error must be caught by test_relerr_1. end; if st_ovfl then if ca_ovfl or else abs( ca) > big_real then return 0.0; -- benefit of the doubt else return total_miss; end if; elsif abs( st) > big_real then -- permit either ca_ovfl or small relerr if ca_ovfl then return 0.0; else return abs( ( ca - st) / st) / epsi; end if; elsif ca_ovfl then return total_miss; elsif st = 0.0 then -- require exactness for zero result if ca = 0.0 then return 0.0; else return total_miss; end if; elsif abs( st) < small_real / (1.0 - eps) then -- underflow permitted if st < 0.0 then st := - st; ca := - ca; end if; if 0.0 <= ca and ca <= small_real then return 0.0; else return abs( ( ca - st) / st) / epsi; end if; else -- normal case return abs( ( ca - st) / st) / epsi; end if; end get_relerr; begin -- test_relerr_2 domain_scaling; for iter in 1 .. niter loop arg1 := get_rand_arg1; arg2 := get_rand_arg2; begin relerr := get_relerr( arg1, arg2); if relerr > max_err then max_err := relerr; max_err_arg1 := long_real( arg1); max_err_arg2 := long_real( arg2); end if; sum_err := sum_err + relerr; exception when constraint_error| numeric_error| storage_error| program_error=> raise; when others=> -- assume this was argument_error from the math routines null; end; end loop; avg_relerr := sum_err / long_real( niter); max_relerr := max_err; worst_arg1 := max_err_arg1; worst_arg2 := max_err_arg2; end test_relerr_2;