with float_text_io; use float_text_io;
with text_io; use text_io;
with elementary_functions; use elementary_functions;
procedure test_fit is
  x : float;
  y : float;
  z : float;

  function distr(x:float) return float is
  begin
    return x*exp(-x);
  end distr;

  function intrg(x:float) return float is
  begin
    return 1.0-(x+1.0)*exp(-x);
  end intrg;

  function pois(r:float) return float is -- 0 <= r <= 1
    x1, x2, y1, y2 : float;
  begin
    if r<=0.0 then return 0.0; end if;
    if r>0.9999 then return 12.0; end if;
    x1 := 1.0;
    y2 := r;
    for i in 1..20 loop -- just to prevent run away
      y1 := 1.0-(x1+1.0)*exp(-x1);
      x2 := x1+(y2-y1)/(x1*exp(-x1));
      exit when abs(x2-x1)/(abs(x2)+abs(x1))< float'epsilon;
      exit when abs(y2-y1)/(abs(y2)+abs(y1))< float'epsilon;
      x1 := x2;
    end loop;
    return x2;
  end pois;

begin
    x := 0.0;
    y := pois(x);
    z := x-intrg(y);
    put(x); put(" pois="); put(y); put(" err="); put(z); new_line;

    x := 0.00001;
    y := pois(x);
    z := x-intrg(y);
    put(x); put(" pois="); put(y); put(" err="); put(z); new_line;

    x := 0.0001;
    y := pois(x);
    z := x-intrg(y);
    put(x); put(" pois="); put(y); put(" err="); put(z); new_line;

    x := 0.001;
    y := pois(x);
    z := x-intrg(y);
    put(x); put(" pois="); put(y); put(" err="); put(z); new_line;

  for I in 1..100 loop
    x := float(I)/100.0;
    y := pois(x);
    z := x-intrg(y);
    put(x); put(" pois="); put(y); put(" err="); put(z); new_line;
  end loop;

    x := 0.995;
    y := pois(x);
    z := x-intrg(y);
    put(x); put(" pois="); put(y); put(" err="); put(z); new_line;

    x := 0.998;
    y := pois(x);
    z := x-intrg(y);
    put(x); put(" pois="); put(y); put(" err="); put(z); new_line;

    x := 0.999;
    y := pois(x);
    z := x-intrg(y);
    put(x); put(" pois="); put(y); put(" err="); put(z); new_line;

    x := 0.9995;
    y := pois(x);
    z := x-intrg(y);
    put(x); put(" pois="); put(y); put(" err="); put(z); new_line;

    x := 0.9999;
    y := pois(x);
    z := x-intrg(y);
    put(x); put(" pois="); put(y); put(" err="); put(z); new_line;

    x := 0.99991;
    y := pois(x);
    z := x-intrg(y);
    put(x); put(" pois="); put(y); put(" err="); put(z); new_line;

end test_fit;