{: test that exp(ln(x)*y) causes a fatal error if x & y are large enough
   this tests for x & y -ve, small x and y just large enough to overflow }

program d08f( output );

var
   e : real;
   k : real;
   x : real;
   y : real;
   z : real;

function side : real;
begin
   side := k;
   k := k+1;
end; { side }

begin {d08f}

   k := 1;
   x := exp(1*ln(side));
   if abs(x-1) > 1.0e-37 then begin
      writeln( 'unexpected value of exp(ln(1)) is ', x);
   end;
   if k <> 2 then
      writeln( 'ln(x) in exp(ln(x)*y) has side effects' );

   k := 1;
   x := exp(side*ln(1));
   if abs(x-1) > 1.0e-37 then begin
      writeln( 'unexpected value of exp(ln(1)) is ', x);
   end;
   if k <> 2 then
      writeln( 'y in exp(ln(x)*y) for debug has side effects' );

   {find x, largest nr < 1}
   k := 1;
   repeat
      k := 0.5*k;
      x := 1 - k/2;
   until x >= 1;
   x := x-k;
   {writeln( 'k is ', k, ', x-k = x is ', x = x-k );}

   {find e, machine precision}
   e := 1;
   repeat
      e := 0.5*e;
      z := 1 + e/2;
   until z <= 1;
   writeln('real precision is ', e);

   {write( 'checking for underflow or denormalise ... ' );}
   k := 1;
   repeat
      y := k;
      k := 0.5*k;
      z := k*(1+e);
   until (k = 0) or (z = k);
   writeln('smallest real is ', y);

   z := 4*x/y; {assuming ieee floating point}
   writeln('max real is ', z);

   {find largest abs(y), st y<0 & y*ln(x) <= ln(z)}
   y := ln(z)/ln(x);
   k := abs(y)/2;
   if y*ln(x) < ln(z) then begin
      {writeln('y*ln(x) < ln(z)');}
      while (y-k)*ln(x) >= ln(z) do begin
         k := k/2;
      end;
      y := y-k;
   end
   else begin
      {writeln('ln(z)*y >= ln(x)');}
      while (y+k/2)*ln(x) < ln(z) do begin
         k := k/2;
      end;
      y := y+k;
   end;
   {writeln('y*ln(x) > ln(z) is ', y*ln(x) > ln(z),
           ', y*ln(x) - ln(z) is ', y*ln(x) - ln(z));}

   { this shouldn't overflow ... }
   writeln('exp(ln(1-', 1-x, ')*', y, ') is ', exp(ln(x)*y));

   {find smallest abs(y), st y<0 & y*ln(x) > ln(z)}
   y := ln(z)/ln(x);
   k := abs(y)/2;
   if y*ln(x) <= ln(z) then begin
      while (y-k/2)*ln(x) > ln(z) do begin
         k := k/2;
      end;
      y := y-k;
   end
   else begin
      while (y+k)*ln(x) <= ln(z) do begin
         k := k/2;
      end;
      y := y+k;
   end;

   writeln( 'checking for overflow, expecting to generate a fatal error' );
   writeln('exp(ln(1-', 1-x, ')*', y, ') is ',
           exp(ln(x)*y)); {!!! exp overflows}

   writeln( 'fails: exp(ln(z)*y) overflow not detected for large y' );
end. { d08f }
