------------------------------------------------------------------------------
--                                                                          --
--                         GNAT RUN-TIME COMPONENTS                         --
--                                                                          --
--            ADA.NUMERICS.GENERIC_COMPLEX_ELEMENTARY_FUNCTIONS             --
--                                                                          --
--                                 B o d y                                  --
--                                                                          --
--          Copyright (C) 1992-2013, Free Software Foundation, Inc.         --
--                                                                          --
-- GNAT is free software;  you can  redistribute it  and/or modify it under --
-- terms of the  GNU General Public License as published  by the Free Soft- --
-- ware  Foundation;  either version 3,  or (at your option) any later ver- --
-- sion.  GNAT is distributed in the hope that it will be useful, but WITH- --
-- OUT ANY WARRANTY;  without even the  implied warranty of MERCHANTABILITY --
-- or FITNESS FOR A PARTICULAR PURPOSE.                                     --
--                                                                          --
-- As a special exception under Section 7 of GPL version 3, you are granted --
-- additional permissions described in the GCC Runtime Library Exception,   --
-- version 3.1, as published by the Free Software Foundation.               --
--                                                                          --
-- You should have received a copy of the GNU General Public License and    --
-- a copy of the GCC Runtime Library Exception along with this program;     --
-- see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see    --
-- <http://www.gnu.org/licenses/>.                                          --
--                                                                          --
-- GNAT was originally developed  by the GNAT team at  New York University. --
-- Extensive contributions were provided by Ada Core Technologies Inc.      --
--                                                                          --
------------------------------------------------------------------------------

with Ada.Numerics.Generic_Elementary_Functions;

package body Ada.Numerics.Generic_Complex_Elementary_Functions is

   package Elementary_Functions is new
      Ada.Numerics.Generic_Elementary_Functions (Real'Base);
   use Elementary_Functions;

   PI      : constant := 3.14159_26535_89793_23846_26433_83279_50288_41971;
   PI_2    : constant := PI / 2.0;
   Sqrt_Two : constant := 1.41421_35623_73095_04880_16887_24209_69807_85696;
   Log_Two : constant := 0.69314_71805_59945_30941_72321_21458_17656_80755;

   subtype T is Real'Base;

   Epsilon                 : constant T := 2.0      ** (1 - T'Model_Mantissa);
   Square_Root_Epsilon     : constant T := Sqrt_Two ** (1 - T'Model_Mantissa);
   Inv_Square_Root_Epsilon : constant T := Sqrt_Two ** (T'Model_Mantissa - 1);
   Root_Root_Epsilon       : constant T := Sqrt_Two **
                                                 ((1 - T'Model_Mantissa) / 2);
   Log_Inverse_Epsilon_2   : constant T := T (T'Model_Mantissa - 1) / 2.0;

   Complex_Zero : constant Complex := (0.0,  0.0);
   Complex_One  : constant Complex := (1.0,  0.0);
   Complex_I    : constant Complex := (0.0,  1.0);
   Half_Pi      : constant Complex := (PI_2, 0.0);

   --------
   -- ** --
   --------

   function "**" (Left : Complex; Right : Complex) return Complex is
   begin
      if Re (Right) = 0.0
        and then Im (Right) = 0.0
        and then Re (Left)  = 0.0
        and then Im (Left)  = 0.0
      then
         raise Argument_Error;

      elsif Re (Left) = 0.0
        and then Im (Left) = 0.0
        and then Re (Right) < 0.0
      then
         raise Constraint_Error;

      elsif Re (Left) = 0.0 and then Im (Left) = 0.0 then
         return Left;

      elsif Right = (0.0, 0.0)  then
         return Complex_One;

      elsif Re (Right) = 0.0 and then Im (Right) = 0.0 then
         return 1.0 + Right;

      elsif Re (Right) = 1.0 and then Im (Right) = 0.0 then
         return Left;

      else
         return Exp (Right * Log (Left));
      end if;
   end "**";

   function "**" (Left : Real'Base; Right : Complex) return Complex is
   begin
      if Re (Right) = 0.0 and then Im (Right) = 0.0 and then Left = 0.0 then
         raise Argument_Error;

      elsif Left = 0.0 and then Re (Right) < 0.0 then
         raise Constraint_Error;

      elsif Left = 0.0 then
         return Compose_From_Cartesian (Left, 0.0);

      elsif Re (Right) = 0.0 and then Im (Right) = 0.0 then
         return Complex_One;

      elsif Re (Right) = 1.0 and then Im (Right) = 0.0 then
         return Compose_From_Cartesian (Left, 0.0);

      else
         return Exp (Log (Left) * Right);
      end if;
   end "**";

   function "**" (Left : Complex; Right : Real'Base) return Complex is
   begin
      if Right = 0.0
        and then Re (Left) = 0.0
        and then Im (Left) = 0.0
      then
         raise Argument_Error;

      elsif Re (Left) = 0.0
        and then Im (Left) = 0.0
        and then Right < 0.0
      then
         raise Constraint_Error;

      elsif Re (Left) = 0.0 and then Im (Left) = 0.0 then
         return Left;

      elsif Right = 0.0 then
         return Complex_One;

      elsif Right = 1.0 then
         return Left;

      else
         return Exp (Right * Log (Left));
      end if;
   end "**";

   ------------
   -- Arccos --
   ------------

   function Arccos (X : Complex) return Complex is
      Result : Complex;

   begin
      if X = Complex_One then
         return Complex_Zero;

      elsif abs Re (X) < Square_Root_Epsilon and then
         abs Im (X) < Square_Root_Epsilon
      then
         return Half_Pi - X;

      elsif abs Re (X) > Inv_Square_Root_Epsilon or else
            abs Im (X) > Inv_Square_Root_Epsilon
      then
         return -2.0 * Complex_I * Log (Sqrt ((1.0 + X) / 2.0) +
                            Complex_I * Sqrt ((1.0 - X) / 2.0));
      end if;

      Result := -Complex_I * Log (X + Complex_I * Sqrt (1.0 - X * X));

      if Im (X) = 0.0
        and then abs Re (X) <= 1.00
      then
         Set_Im (Result, Im (X));
      end if;

      return Result;
   end Arccos;

   -------------
   -- Arccosh --
   -------------

   function Arccosh (X : Complex) return Complex is
      Result : Complex;

   begin
      if X = Complex_One then
         return Complex_Zero;

      elsif abs Re (X) < Square_Root_Epsilon and then
         abs Im (X) < Square_Root_Epsilon
      then
         Result := Compose_From_Cartesian (-Im (X), -PI_2 + Re (X));

      elsif abs Re (X) > Inv_Square_Root_Epsilon or else
            abs Im (X) > Inv_Square_Root_Epsilon
      then
         Result := Log_Two + Log (X);

      else
         Result := 2.0 * Log (Sqrt ((1.0 + X) / 2.0) +
                              Sqrt ((X - 1.0) / 2.0));
      end if;

      if Re (Result) <= 0.0 then
         Result := -Result;
      end if;

      return Result;
   end Arccosh;

   ------------
   -- Arccot --
   ------------

   function Arccot (X : Complex) return Complex is
      Xt : Complex;

   begin
      if abs Re (X) < Square_Root_Epsilon and then
         abs Im (X) < Square_Root_Epsilon
      then
         return Half_Pi - X;

      elsif abs Re (X) > 1.0 / Epsilon or else
            abs Im (X) > 1.0 / Epsilon
      then
         Xt := Complex_One  /  X;

         if Re (X) < 0.0 then
            Set_Re (Xt, PI - Re (Xt));
            return Xt;
         else
            return Xt;
         end if;
      end if;

      Xt := Complex_I * Log ((X - Complex_I) / (X + Complex_I)) / 2.0;

      if Re (Xt) < 0.0 then
         Xt := PI + Xt;
      end if;

      return Xt;
   end Arccot;

   --------------
   -- Arccoth --
   --------------

   function Arccoth (X : Complex) return Complex is
      R : Complex;

   begin
      if X = (0.0, 0.0) then
         return Compose_From_Cartesian (0.0, PI_2);

      elsif abs Re (X) < Square_Root_Epsilon
         and then abs Im (X) < Square_Root_Epsilon
      then
         return PI_2 * Complex_I + X;

      elsif abs Re (X) > 1.0 / Epsilon or else
            abs Im (X) > 1.0 / Epsilon
      then
         if Im (X) > 0.0 then
            return (0.0, 0.0);
         else
            return PI * Complex_I;
         end if;

      elsif Im (X) = 0.0 and then Re (X) = 1.0 then
         raise Constraint_Error;

      elsif Im (X) = 0.0 and then Re (X) = -1.0 then
         raise Constraint_Error;
      end if;

      begin
         R := Log ((1.0 + X) / (X - 1.0)) / 2.0;

      exception
         when Constraint_Error =>
            R := (Log (1.0 + X) - Log (X - 1.0)) / 2.0;
      end;

      if Im (R) < 0.0 then
         Set_Im (R, PI + Im (R));
      end if;

      if Re (X) = 0.0 then
         Set_Re (R, Re (X));
      end if;

      return R;
   end Arccoth;

   ------------
   -- Arcsin --
   ------------

   function Arcsin (X : Complex) return Complex is
      Result : Complex;

   begin
      --  For very small argument, sin (x) = x

      if abs Re (X) < Square_Root_Epsilon and then
         abs Im (X) < Square_Root_Epsilon
      then
         return X;

      elsif abs Re (X) > Inv_Square_Root_Epsilon or else
            abs Im (X) > Inv_Square_Root_Epsilon
      then
         Result := -Complex_I * (Log (Complex_I * X) + Log (2.0 * Complex_I));

         if Im (Result) > PI_2 then
            Set_Im (Result, PI - Im (X));

         elsif Im (Result) < -PI_2 then
            Set_Im (Result, -(PI + Im (X)));
         end if;

         return Result;
      end if;

      Result := -Complex_I * Log (Complex_I * X + Sqrt (1.0 - X * X));

      if Re (X) = 0.0 then
         Set_Re (Result, Re (X));

      elsif Im (X) = 0.0
        and then abs Re (X) <= 1.00
      then
         Set_Im (Result, Im (X));
      end if;

      return Result;
   end Arcsin;

   -------------
   -- Arcsinh --
   -------------

   function Arcsinh (X : Complex) return Complex is
      Result : Complex;

   begin
      if abs Re (X) < Square_Root_Epsilon and then
         abs Im (X) < Square_Root_Epsilon
      then
         return X;

      elsif abs Re (X) > Inv_Square_Root_Epsilon or else
            abs Im (X) > Inv_Square_Root_Epsilon
      then
         Result := Log_Two + Log (X); -- may have wrong sign

         if (Re (X) < 0.0 and then Re (Result) > 0.0)
           or else (Re (X) > 0.0 and then Re (Result) < 0.0)
         then
            Set_Re (Result, -Re (Result));
         end if;

         return Result;
      end if;

      Result := Log (X + Sqrt (1.0 + X * X));

      if Re (X) = 0.0 then
         Set_Re (Result, Re (X));
      elsif Im  (X) = 0.0 then
         Set_Im (Result, Im  (X));
      end if;

      return Result;
   end Arcsinh;

   ------------
   -- Arctan --
   ------------

   function Arctan (X : Complex) return Complex is
   begin
      if abs Re (X) < Square_Root_Epsilon and then
         abs Im (X) < Square_Root_Epsilon
      then
         return X;

      else
         return -Complex_I * (Log (1.0 + Complex_I * X)
                            - Log (1.0 - Complex_I * X)) / 2.0;
      end if;
   end Arctan;

   -------------
   -- Arctanh --
   -------------

   function Arctanh (X : Complex) return Complex is
   begin
      if abs Re (X) < Square_Root_Epsilon and then
         abs Im (X) < Square_Root_Epsilon
      then
         return X;
      else
         return (Log (1.0 + X) - Log (1.0 - X)) / 2.0;
      end if;
   end Arctanh;

   ---------
   -- Cos --
   ---------

   function Cos (X : Complex) return Complex is
   begin
      return
        Compose_From_Cartesian
          (Cos (Re (X))  * Cosh (Im (X)),
           -(Sin (Re (X)) * Sinh (Im (X))));
   end Cos;

   ----------
   -- Cosh --
   ----------

   function Cosh (X : Complex) return Complex is
   begin
      return
        Compose_From_Cartesian
          (Cosh (Re (X)) * Cos (Im (X)),
           Sinh (Re (X)) * Sin (Im (X)));
   end Cosh;

   ---------
   -- Cot --
   ---------

   function Cot (X : Complex) return Complex is
   begin
      if abs Re (X) < Square_Root_Epsilon and then
         abs Im (X) < Square_Root_Epsilon
      then
         return Complex_One  /  X;

      elsif Im (X) > Log_Inverse_Epsilon_2 then
         return -Complex_I;

      elsif Im (X) < -Log_Inverse_Epsilon_2 then
         return Complex_I;
      end if;

      return Cos (X) / Sin (X);
   end Cot;

   ----------
   -- Coth --
   ----------

   function Coth (X : Complex) return Complex is
   begin
      if abs Re (X) < Square_Root_Epsilon and then
         abs Im (X) < Square_Root_Epsilon
      then
         return Complex_One  /  X;

      elsif Re (X) > Log_Inverse_Epsilon_2 then
         return Complex_One;

      elsif Re (X) < -Log_Inverse_Epsilon_2 then
         return -Complex_One;

      else
         return Cosh (X) / Sinh (X);
      end if;
   end Coth;

   ---------
   -- Exp --
   ---------

   function Exp (X : Complex) return Complex is
      EXP_RE_X : constant Real'Base := Exp (Re (X));

   begin
      return Compose_From_Cartesian (EXP_RE_X * Cos (Im (X)),
                                     EXP_RE_X * Sin (Im (X)));
   end Exp;

   function Exp (X : Imaginary) return Complex is
      ImX : constant Real'Base := Im (X);

   begin
      return Compose_From_Cartesian (Cos (ImX), Sin (ImX));
   end Exp;

   ---------
   -- Log --
   ---------

   function Log (X : Complex) return Complex is
      ReX : Real'Base;
      ImX : Real'Base;
      Z   : Complex;

   begin
      if Re (X) = 0.0 and then Im (X) = 0.0 then
         raise Constraint_Error;

      elsif abs (1.0 - Re (X)) < Root_Root_Epsilon
        and then abs Im (X) < Root_Root_Epsilon
      then
         Z := X;
         Set_Re (Z, Re (Z) - 1.0);

         return (1.0 - (1.0 / 2.0 -
                       (1.0 / 3.0 - (1.0 / 4.0) * Z) * Z) * Z) * Z;
      end if;

      begin
         ReX := Log (Modulus (X));

      exception
         when Constraint_Error =>
            ReX := Log (Modulus (X / 2.0)) - Log_Two;
      end;

      ImX := Arctan (Im (X), Re (X));

      if ImX > PI then
         ImX := ImX - 2.0 * PI;
      end if;

      return Compose_From_Cartesian (ReX, ImX);
   end Log;

   ---------
   -- Sin --
   ---------

   function Sin (X : Complex) return Complex is
   begin
      if abs Re (X) < Square_Root_Epsilon
           and then
         abs Im (X) < Square_Root_Epsilon
      then
         return X;
      end if;

      return
        Compose_From_Cartesian
          (Sin (Re (X)) * Cosh (Im (X)),
           Cos (Re (X)) * Sinh (Im (X)));
   end Sin;

   ----------
   -- Sinh --
   ----------

   function Sinh (X : Complex) return Complex is
   begin
      if abs Re (X) < Square_Root_Epsilon and then
         abs Im (X) < Square_Root_Epsilon
      then
         return X;

      else
         return Compose_From_Cartesian (Sinh (Re (X)) * Cos (Im (X)),
                                        Cosh (Re (X)) * Sin (Im (X)));
      end if;
   end Sinh;

   ----------
   -- Sqrt --
   ----------

   function Sqrt (X : Complex) return Complex is
      ReX : constant Real'Base := Re (X);
      ImX : constant Real'Base := Im (X);
      XR  : constant Real'Base := abs Re (X);
      YR  : constant Real'Base := abs Im (X);
      R   : Real'Base;
      R_X : Real'Base;
      R_Y : Real'Base;

   begin
      --  Deal with pure real case, see (RM G.1.2(39))

      if ImX = 0.0 then
         if ReX > 0.0 then
            return
              Compose_From_Cartesian
                (Sqrt (ReX), 0.0);

         elsif ReX = 0.0 then
            return X;

         else
            return
              Compose_From_Cartesian
                (0.0, Real'Copy_Sign (Sqrt (-ReX), ImX));
         end if;

      elsif ReX = 0.0 then
         R_X := Sqrt (YR / 2.0);

         if ImX > 0.0 then
            return Compose_From_Cartesian (R_X, R_X);
         else
            return Compose_From_Cartesian (R_X, -R_X);
         end if;

      else
         R  := Sqrt (XR ** 2 + YR ** 2);

         --  If the square of the modulus overflows, try rescaling the
         --  real and imaginary parts. We cannot depend on an exception
         --  being raised on all targets.

         if R > Real'Base'Last then
            raise Constraint_Error;
         end if;

         --  We are solving the system

         --  XR = R_X ** 2 - Y_R ** 2      (1)
         --  YR = 2.0 * R_X * R_Y          (2)
         --
         --  The symmetric solution involves square roots for both R_X and
         --  R_Y, but it is more accurate to use the square root with the
         --  larger argument for either R_X or R_Y, and equation (2) for the
         --  other.

         if ReX < 0.0 then
            R_Y := Sqrt (0.5 * (R - ReX));
            R_X := YR / (2.0 * R_Y);

         else
            R_X := Sqrt (0.5 * (R + ReX));
            R_Y := YR / (2.0 * R_X);
         end if;
      end if;

      if Im (X) < 0.0 then -- halve angle, Sqrt of magnitude
         R_Y := -R_Y;
      end if;
      return Compose_From_Cartesian (R_X, R_Y);

   exception
      when Constraint_Error =>

         --  Rescale and try again

         R := Modulus (Compose_From_Cartesian (Re (X / 4.0), Im (X / 4.0)));
         R_X := 2.0 * Sqrt (0.5 * R + 0.5 * Re (X / 4.0));
         R_Y := 2.0 * Sqrt (0.5 * R - 0.5 * Re (X / 4.0));

         if Im (X) < 0.0 then -- halve angle, Sqrt of magnitude
            R_Y := -R_Y;
         end if;

         return Compose_From_Cartesian (R_X, R_Y);
   end Sqrt;

   ---------
   -- Tan --
   ---------

   function Tan (X : Complex) return Complex is
   begin
      if abs Re (X) < Square_Root_Epsilon and then
         abs Im (X) < Square_Root_Epsilon
      then
         return X;

      elsif Im (X) > Log_Inverse_Epsilon_2 then
         return Complex_I;

      elsif Im (X) < -Log_Inverse_Epsilon_2 then
         return -Complex_I;

      else
         return Sin (X) / Cos (X);
      end if;
   end Tan;

   ----------
   -- Tanh --
   ----------

   function Tanh (X : Complex) return Complex is
   begin
      if abs Re (X) < Square_Root_Epsilon and then
         abs Im (X) < Square_Root_Epsilon
      then
         return X;

      elsif Re (X) > Log_Inverse_Epsilon_2 then
         return Complex_One;

      elsif Re (X) < -Log_Inverse_Epsilon_2 then
         return -Complex_One;

      else
         return Sinh (X) / Cosh (X);
      end if;
   end Tanh;

end Ada.Numerics.Generic_Complex_Elementary_Functions;