--     ALGORITHM 404 COLLECTED ALGORITHMS FROM ACM.
--     ALGORITHM APPEARED IN COMM. ACM, VOL. 14, NO. 01,
--     P. 048.
--    converted to Ada using existing Complex_Types and
--    Complex_Elementary_Functions packages by Jon Squire  8/25/94
--    some goto's remain.

with Text_IO; use Text_IO; -- can be replaced by raising exception
with Complex_Types; Use Complex_Types;
with Elementary_Functions; use Elementary_Functions;
with Complex_Elementary_Functions; use Complex_Elementary_Functions;

function Complex_Gamma(ZZ : Complex) return Complex is
  ZM, T, TT, SUM, TERM, DEN, A ,CGAMMA : Complex;
  -- COEFFICIENTS IN STIRLING*S APPROXIMATION FOR LN(GAMMA(T))
  C : array(1..12) of float := (1.0/12.0, 
          -1.0/360.0, 1.0/1260.0, -1.0/1680.0, 1.0/1188.0, -691.0/360360.0, 
          1.0/156.0, -3617.0/122400.0, 43867.0/244188.0, -174611.0/125400.0, 
          77683.0/5796.0, 0.0);
  REFLEK : Boolean := true;
  PI : Complex := (3.141593,0.0);
  Z : Complex := ZZ;
  X  : Float := RE(Z);
  Y  : Float := IM(Z);
  XDIST : Float;
  M  : Integer;
  -- TOL := LIMIT OF PRECISION OF COMPUTER SYSTEM IN SINGLE PRECISI
  TOL : Float := 1.0E-7;
begin
  -- DETERMINE WHETHER Z IS TOO CLOSE TO A POLE
  -- CHECK WHETHER TOO CLOSE TO ORIGIN
  if X < TOL then -- 20
    -- FIND THE NEAREST POLE AND COMPUTE DISTANCE TO IT
    XDIST := X-Float(Integer(X)); --Float(Integer(X-0.5));
    ZM := (XDIST,Y);
    if ABS(ZM) < TOL then --10
      -- IF Z IS TOO CLOSE TO A POLE, PRINT ERROR MESSAGE AND RETURN
      -- WITH CGAMMA := (1.E7,0.0E0)
      -- WRITE(IOUT,900) Z
      Put_line("ARGUMENT OF GAMMA FUNCTION IS TOO CLOSE TO A POLE");
      CGAMMA := (1.0E7,0.0);
      return CGAMMA;
    end if; --10
    if X < 0.0 then
      -- FOR REAL(Z) NEGATIVE EMPLOY THE REFLECTION FORMULA
      -- GAMMA(Z) := PI/(SIN(PI*Z)*GAMMA(1-Z))
      -- AND COMPUTE GAMMA(1-Z).  NOTE REFLEK IS A TAG TO INDICATE THA
      -- THIS RELATION MUST BE USED LATER.
      REFLEK := false;
      Z := (1.0,0.0)-Z;
      X := 1.0-X;
      Y := -Y;
    end if;
  end if; --20
  -- IF Z IS NOT TOO CLOSE TO A POLE, MAKE REAL(Z)>10 AND ARG(Z)<P
  M := 0;
  while X < 10.0 loop
      X := X + 1.0;
      M := M + 1;
  end loop;
  while abs Y >= X loop
      X := X + 1.0;
      M := M + 1;
  end loop;
  T := (X,Y);
  TT := T*T;
  DEN := T;
  SUM := (T-(0.5,0.0))*log(T)-T+(0.5* log(2.0*3.14159),0.0);
  for J in 1..11 loop --70
    TERM := C(J)/DEN;
    -- TEST REAL AND IMAGINARY PARTS OF LN(GAMMA(Z)) SEPARATELY FOR
    -- CONVERGENCE.  IF Z IS REAL SKIP IMAGINARY PART OF CHECK.
      if RE(SUM) /= 0.0 and then abs(RE(TERM)/RE(SUM)) < TOL then --80
      if Y = 0.0 then GOTO L100; end if;
      if IM(SUM) /= 0.0 and then ABS(IM(TERM)/IM(SUM)) < TOL then GOTO L100; end if;
    end if; --80
    SUM := SUM + TERM;
    DEN := DEN*TT;
  end loop; --70
  -- STIRLING*S SERIES DID NOT CONVERGE.  PRINT ERROR MESSAGE AND
  -- PROCEDE.
--90    WRITE(IOUT,910) Z
  Put_Line("No Convergence");
  -- RECURSION RELATION USED TO OBTAIN LN(GAMMA(Z))
  -- LN(GAMMA(Z)) := LN(GAMMA(Z+M)/(Z*(Z+1)*...*(Z+M-1)))
  -- := LN(GAMMA(Z+M)-LN(Z)-LN(Z+1)-...-LN(Z+M
<<L100>> 
  if M /= 0 then --120
    for I in 1..M loop
      A := (Float(I)-1.0,0.0);
      SUM := SUM-log(Z+A);
    end loop;
  end if; -- 120
  -- CHECK TO SEE IF REFLECTION FORMULA SHOULD BE USED
  if not REFLEK then --130
      SUM := log(PI/sin(PI*Z))-SUM;
      Z := (1.0,0.0) -Z;
  end if; --130
  CGAMMA := exp(SUM);
  return CGAMMA;
--910   FORMAT(44H ERROR - STIRLING*S SERIES HAS NOT CONVERGED/14X,4HZ := ,
--     12E14.7)
end Complex_Gamma;