with complex_types; use complex_types;
with elementary_functions; use elementary_functions;
with test_proto_sqrt_pkg; use test_proto_sqrt_pkg;

  function proto_sqrt(Z:complex) return complex is
  re : float;
  im : float;
  r  : float;
  x  : float := abs RE(z);
  y  : float := abs IM(z);
  a  : float;
  EPSILON : float;
  SQUARE_ROOT_EPSILON : float;

begin 
  EPSILON := FLOAT'EPSILON;
  while 1.0+EPSILON/2.0 /= 1.0 loop
    EPSILON := EPSILON/2.0;
  end loop;
  SQUARE_ROOT_EPSILON := 32.0*SQRT(EPSILON);
  if x > y then
    a := y/x;
    a := a*a;
    if a < SQUARE_ROOT_EPSILON then
      if RE(Z) > 0.0 then
        where := X_GT_Y_EPS_POS;
        re := SQRT(x+0.5*(0.5*y*y/x-(y*y/x)*a/8.0));
        im := SQRT(0.5*(0.5*y*y/x-(y*y/x)*a/8.0));
      else
        where := X_GT_Y_EPS_NEG;
        re := SQRT(0.5*(0.5*y*y/x-(y*y/x)*a/8.0));
        im := SQRT(x+0.5*(0.5*y*y/x-(y*y/x)*a/8.0));
      end if;
    else
      where := X_GT_Y;
      r  := x*sqrt(1.0+a);
      re := sqrt(0.5*(r+RE(Z)));
      im := sqrt(0.5*(r-RE(Z)));
    end if;
  else -- y > x
    a := x/y;
    a := a*a;
    if a < SQUARE_ROOT_EPSILON then -- y >> x
      where := Y_GT_X_EPS;
      r  := y + 0.5*x*x/y-0.125*(x*x/y)*a;
      re := SQRT(0.5*(y+(0.5*x*x/y+RE(Z))));
      im := SQRT(0.5*(y+(0.5*x*x/y-RE(Z))));
    else
      where := Y_GT_X;
      r  := y*sqrt(1.0+a);
      re := sqrt(0.5*(r+RE(Z)));
      im := sqrt(0.5*(r-RE(Z)));
    end if;
  end if;
  if IM(Z) < 0.0 then
    im := -im;
  end if;
  return compose_from_cartesian(re,im);
end proto_sqrt;