--  -----------------------------------------------------------------------
--  Title:       test_relerr_1
--  Last Mod:    Tue Oct  2 14:38:16 1990
--  Author:      Vincent Broman <broman@schroeder>
--  Visibility:  instantiated and used by a test driver
--  Description: 
--		 A generic procedure for testing relative accuracy of
--		 floating point functions of one argument.
--		 Random arguments are chosen from an interval 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 <>;
   with function cand( x: in real) return real;	-- candidate for testing
   with function std( x: in long_real) return long_real; -- std. for comparison
   eps: long_real;			-- required rel acc
   exp_dist_arg: boolean;		-- want expon dist on random args tested
procedure test_relerr_1( niter: in positive; -- nbr of tests in run
			 minarg, maxarg: in real; -- bounds on domain
			 avg_relerr: out long_real; -- L1 err norm / niter
			 max_relerr: out long_real; -- Linf err norm
			 worst_arg: out long_real); -- arg realizing max err

procedure test_relerr_1( niter: in positive; -- nbr of tests in run
			 minarg, maxarg: in real; -- bounds on domain
			 avg_relerr: out long_real; -- L1 err norm / niter
			 max_relerr: out long_real; -- Linf err norm
			 worst_arg: out long_real) is -- arg realizing max err
-- 
   max_err_arg: long_real := 0.0;
   max_err: long_real := 0.0;
   sum_err: long_real := 0.0;
   relerr: long_real;
   arg: real;

   -- following used only if exp_dist_arg true
   aminarg: real;
   amaxarg: real;
   negate: boolean;
   lowexp: integer;
   exprg: real;
   

   package primf is new generic_primitive_functions( real, integer);  use primf;


   procedure domain_scaling is
-- 

   begin
      if exp_dist_arg then
	 aminarg := abs( minarg);
	 amaxarg := abs( maxarg);
	 negate := maxarg < 0.0;
	 if aminarg = 0.0 then
	    aminarg := real'small;
	 end if;
	 lowexp := exponent( aminarg) - 1;
	 exprg := real( exponent( amaxarg) - lowexp);
      end if;
   end domain_scaling;

   function get_rand_arg return real is
-- 
      arg: real;
      mant: real;
      exd: integer;

   begin				-- get_rand_arg
      loop
	 begin
	    if exp_dist_arg then
	       mant := real( random_number.next) + 1.0;
	       exd := integer( floor( exprg * real( random_number.next)));
	       arg := scale( mant, lowexp + exd);
	       if negate then
		  arg := -arg;
	       end if;
	    else
	       arg := minarg + (maxarg - minarg) * 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_arg then
	       exit when aminarg <= abs( arg) and abs( arg) <= amaxarg;
	    else
	       exit when minarg <= arg and arg <= maxarg;
	    end if;
	 exception
	    when numeric_error| constraint_error=>
	       null;
	 end;
      end loop;
      return arg;
   end get_rand_arg;
   
   function get_relerr( x: 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);

   begin				-- get_relerr
      begin
      	 ca := long_real( cand( x));
      exception
	 when numeric_error| constraint_error=>
	    ca_ovfl := true;
	 -- when argument_error we have a bug
      end;
      begin
	 st := std( long_real( x));
      exception
	 when numeric_error| constraint_error=>
	    st_ovfl := true;
	 -- when argument_error we have a bug
      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) / eps;
	 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) / eps;
	 end if;
      else
	 -- normal case
      	 return abs( ( ca - st) / st) / eps;
      end if;
   end get_relerr;

begin					-- test_relerr_1
   domain_scaling;
   for iter in 1 .. niter loop
      arg := get_rand_arg;
      begin
	 relerr := get_relerr( arg);
	 if relerr > max_err then
	    max_err := relerr;
	    max_err_arg := long_real( arg);
	 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;			-- only for two-arg [co]tangents
      end;
   end loop;
   avg_relerr := sum_err / long_real( niter);
   max_relerr := max_err;
   worst_arg := max_err_arg;
end test_relerr_1;