--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/gcc/ada/libgnat/a-ngcefu.adb	Fri Oct 27 22:46:09 2017 +0900
@@ -0,0 +1,710 @@
+------------------------------------------------------------------------------
+--                                                                          --
+--                         GNAT RUN-TIME COMPONENTS                         --
+--                                                                          --
+--            ADA.NUMERICS.GENERIC_COMPLEX_ELEMENTARY_FUNCTIONS             --
+--                                                                          --
+--                                 B o d y                                  --
+--                                                                          --
+--          Copyright (C) 1992-2017, 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;