------------------------------------------------------------------------------
--                                                                          --
--                         GNAT RUN-TIME COMPONENTS                         --
--                                                                          --
--                ADA.NUMERICS.GENERIC_ELEMENTARY_FUNCTIONS                 --
--                                                                          --
--                                 B o d y                                  --
--                                                                          --
--          Copyright (C) 1992-2015, 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.      --
--                                                                          --
------------------------------------------------------------------------------

--  This body is specifically for using an Ada interface to C math.h to get
--  the computation engine. Many special cases are handled locally to avoid
--  unnecessary calls or to meet Annex G strict mode requirements.

--  Uses functions sqrt, exp, log, pow, sin, asin, cos, acos, tan, atan, sinh,
--  cosh, tanh from C library via math.h

with Ada.Numerics.Aux;

package body Ada.Numerics.Generic_Elementary_Functions is

   use type Ada.Numerics.Aux.Double;

   Sqrt_Two : constant := 1.41421_35623_73095_04880_16887_24209_69807_85696;
   Log_Two  : constant := 0.69314_71805_59945_30941_72321_21458_17656_80755;

   Half_Log_Two : constant := Log_Two / 2;

   subtype T is Float_Type'Base;
   subtype Double is Aux.Double;

   Two_Pi  : constant T := 2.0 * Pi;
   Half_Pi : constant T := Pi / 2.0;

   Half_Log_Epsilon    : constant T := T (1 - T'Model_Mantissa) * Half_Log_Two;
   Log_Inverse_Epsilon : constant T := T (T'Model_Mantissa - 1) * Log_Two;
   Sqrt_Epsilon        : constant T := Sqrt_Two ** (1 - T'Model_Mantissa);

   -----------------------
   -- Local Subprograms --
   -----------------------

   function Exp_Strict (X : Float_Type'Base) return Float_Type'Base;
   --  Cody/Waite routine, supposedly more precise than the library version.
   --  Currently only needed for Sinh/Cosh on X86 with the largest FP type.

   function Local_Atan
     (Y : Float_Type'Base;
      X : Float_Type'Base := 1.0) return Float_Type'Base;
   --  Common code for arc tangent after cycle reduction

   ----------
   -- "**" --
   ----------

   function "**" (Left, Right : Float_Type'Base) return Float_Type'Base is
      A_Right  : Float_Type'Base;
      Int_Part : Integer;
      Result   : Float_Type'Base;
      R1       : Float_Type'Base;
      Rest     : Float_Type'Base;

   begin
      if Left = 0.0
        and then Right = 0.0
      then
         raise Argument_Error;

      elsif Left < 0.0 then
         raise Argument_Error;

      elsif Right = 0.0 then
         return 1.0;

      elsif Left = 0.0 then
         if Right < 0.0 then
            raise Constraint_Error;
         else
            return 0.0;
         end if;

      elsif Left = 1.0 then
         return 1.0;

      elsif Right = 1.0 then
         return Left;

      else
         begin
            if Right = 2.0 then
               return Left * Left;

            elsif Right = 0.5 then
               return Sqrt (Left);

            else
               A_Right := abs (Right);

               --  If exponent is larger than one, compute integer exponen-
               --  tiation if possible, and evaluate fractional part with more
               --  precision. The relative error is now proportional to the
               --  fractional part of the exponent only.

               if A_Right > 1.0
                 and then A_Right < Float_Type'Base (Integer'Last)
               then
                  Int_Part := Integer (Float_Type'Base'Truncation (A_Right));
                  Result := Left ** Int_Part;
                  Rest := A_Right - Float_Type'Base (Int_Part);

                  --  Compute with two leading bits of the mantissa using
                  --  square roots. Bound  to be better than logarithms, and
                  --  easily extended to greater precision.

                  if Rest >= 0.5 then
                     R1 := Sqrt (Left);
                     Result := Result * R1;
                     Rest := Rest - 0.5;

                     if Rest >= 0.25 then
                        Result := Result * Sqrt (R1);
                        Rest := Rest - 0.25;
                     end if;

                  elsif Rest >= 0.25 then
                     Result := Result * Sqrt (Sqrt (Left));
                     Rest := Rest - 0.25;
                  end if;

                  Result := Result *
                    Float_Type'Base (Aux.Pow (Double (Left), Double (Rest)));

                  if Right >= 0.0 then
                     return Result;
                  else
                     return (1.0 / Result);
                  end if;
               else
                  return
                    Float_Type'Base (Aux.Pow (Double (Left), Double (Right)));
               end if;
            end if;

         exception
            when others =>
               raise Constraint_Error;
         end;
      end if;
   end "**";

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

   --  Natural cycle

   function Arccos (X : Float_Type'Base) return Float_Type'Base is
      Temp : Float_Type'Base;

   begin
      if abs X > 1.0 then
         raise Argument_Error;

      elsif abs X < Sqrt_Epsilon then
         return Pi / 2.0 - X;

      elsif X = 1.0 then
         return 0.0;

      elsif X = -1.0 then
         return Pi;
      end if;

      Temp := Float_Type'Base (Aux.Acos (Double (X)));

      if Temp < 0.0 then
         Temp := Pi + Temp;
      end if;

      return Temp;
   end Arccos;

   --  Arbitrary cycle

   function Arccos (X, Cycle : Float_Type'Base) return Float_Type'Base is
      Temp : Float_Type'Base;

   begin
      if Cycle <= 0.0 then
         raise Argument_Error;

      elsif abs X > 1.0 then
         raise Argument_Error;

      elsif abs X < Sqrt_Epsilon then
         return Cycle / 4.0;

      elsif X = 1.0 then
         return 0.0;

      elsif X = -1.0 then
         return Cycle / 2.0;
      end if;

      Temp := Arctan (Sqrt ((1.0 - X) * (1.0 + X)) / X, 1.0, Cycle);

      if Temp < 0.0 then
         Temp := Cycle / 2.0 + Temp;
      end if;

      return Temp;
   end Arccos;

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

   function Arccosh (X : Float_Type'Base) return Float_Type'Base is
   begin
      --  Return positive branch of Log (X - Sqrt (X * X - 1.0)), or the proper
      --  approximation for X close to 1 or >> 1.

      if X < 1.0 then
         raise Argument_Error;

      elsif X < 1.0 + Sqrt_Epsilon then
         return Sqrt (2.0 * (X - 1.0));

      elsif X > 1.0 / Sqrt_Epsilon then
         return Log (X) + Log_Two;

      else
         return Log (X + Sqrt ((X - 1.0) * (X + 1.0)));
      end if;
   end Arccosh;

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

   --  Natural cycle

   function Arccot
     (X    : Float_Type'Base;
      Y    : Float_Type'Base := 1.0)
      return Float_Type'Base
   is
   begin
      --  Just reverse arguments

      return Arctan (Y, X);
   end Arccot;

   --  Arbitrary cycle

   function Arccot
     (X     : Float_Type'Base;
      Y     : Float_Type'Base := 1.0;
      Cycle : Float_Type'Base)
      return  Float_Type'Base
   is
   begin
      --  Just reverse arguments

      return Arctan (Y, X, Cycle);
   end Arccot;

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

   function Arccoth (X : Float_Type'Base) return Float_Type'Base is
   begin
      if abs X > 2.0 then
         return Arctanh (1.0 / X);

      elsif abs X = 1.0 then
         raise Constraint_Error;

      elsif abs X < 1.0 then
         raise Argument_Error;

      else
         --  1.0 < abs X <= 2.0. One of X + 1.0 and X - 1.0 is exact, the other
         --  has error 0 or Epsilon.

         return 0.5 * (Log (abs (X + 1.0)) - Log (abs (X - 1.0)));
      end if;
   end Arccoth;

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

   --  Natural cycle

   function Arcsin (X : Float_Type'Base) return Float_Type'Base is
   begin
      if abs X > 1.0 then
         raise Argument_Error;

      elsif abs X < Sqrt_Epsilon then
         return X;

      elsif X = 1.0 then
         return Pi / 2.0;

      elsif X = -1.0 then
         return -(Pi / 2.0);
      end if;

      return Float_Type'Base (Aux.Asin (Double (X)));
   end Arcsin;

   --  Arbitrary cycle

   function Arcsin (X, Cycle : Float_Type'Base) return Float_Type'Base is
   begin
      if Cycle <= 0.0 then
         raise Argument_Error;

      elsif abs X > 1.0 then
         raise Argument_Error;

      elsif X = 0.0 then
         return X;

      elsif X = 1.0 then
         return Cycle / 4.0;

      elsif X = -1.0 then
         return -(Cycle / 4.0);
      end if;

      return Arctan (X / Sqrt ((1.0 - X) * (1.0 + X)), 1.0, Cycle);
   end Arcsin;

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

   function Arcsinh (X : Float_Type'Base) return Float_Type'Base is
   begin
      if abs X < Sqrt_Epsilon then
         return X;

      elsif X > 1.0 / Sqrt_Epsilon then
         return Log (X) + Log_Two;

      elsif X < -(1.0 / Sqrt_Epsilon) then
         return -(Log (-X) + Log_Two);

      elsif X < 0.0 then
         return -Log (abs X + Sqrt (X * X + 1.0));

      else
         return Log (X + Sqrt (X * X + 1.0));
      end if;
   end Arcsinh;

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

   --  Natural cycle

   function Arctan
     (Y    : Float_Type'Base;
      X    : Float_Type'Base := 1.0)
      return Float_Type'Base
   is
   begin
      if X = 0.0 and then Y = 0.0 then
         raise Argument_Error;

      elsif Y = 0.0 then
         if X > 0.0 then
            return 0.0;
         else -- X < 0.0
            return Pi * Float_Type'Copy_Sign (1.0, Y);
         end if;

      elsif X = 0.0 then
         return Float_Type'Copy_Sign (Half_Pi, Y);

      else
         return Local_Atan (Y, X);
      end if;
   end Arctan;

   --  Arbitrary cycle

   function Arctan
     (Y     : Float_Type'Base;
      X     : Float_Type'Base := 1.0;
      Cycle : Float_Type'Base)
      return  Float_Type'Base
   is
   begin
      if Cycle <= 0.0 then
         raise Argument_Error;

      elsif X = 0.0 and then Y = 0.0 then
         raise Argument_Error;

      elsif Y = 0.0 then
         if X > 0.0 then
            return 0.0;
         else -- X < 0.0
            return Cycle / 2.0 * Float_Type'Copy_Sign (1.0, Y);
         end if;

      elsif X = 0.0 then
         return Float_Type'Copy_Sign (Cycle / 4.0, Y);

      else
         return Local_Atan (Y, X) *  Cycle / Two_Pi;
      end if;
   end Arctan;

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

   function Arctanh (X : Float_Type'Base) return Float_Type'Base is
      A, B, D, A_Plus_1, A_From_1 : Float_Type'Base;

      Mantissa : constant Integer := Float_Type'Base'Machine_Mantissa;

   begin
      --  The naive formula:

      --     Arctanh (X) := (1/2) * Log  (1 + X) / (1 - X)

      --   is not well-behaved numerically when X < 0.5 and when X is close
      --   to one. The following is accurate but probably not optimal.

      if abs X = 1.0 then
         raise Constraint_Error;

      elsif abs X >= 1.0 - 2.0 ** (-Mantissa) then

         if abs X >= 1.0 then
            raise Argument_Error;
         else

            --  The one case that overflows if put through the method below:
            --  abs X = 1.0 - Epsilon.  In this case (1/2) log (2/Epsilon) is
            --  accurate. This simplifies to:

            return Float_Type'Copy_Sign (
               Half_Log_Two * Float_Type'Base (Mantissa + 1), X);
         end if;

      --  elsif abs X <= 0.5 then
      --  why is above line commented out ???

      else
         --  Use several piecewise linear approximations. A is close to X,
         --  chosen so 1.0 + A, 1.0 - A, and X - A are exact. The two scalings
         --  remove the low-order bits of X.

         A := Float_Type'Base'Scaling (
             Float_Type'Base (Long_Long_Integer
               (Float_Type'Base'Scaling (X, Mantissa - 1))), 1 - Mantissa);

         B := X - A;                --  This is exact; abs B <= 2**(-Mantissa).
         A_Plus_1 := 1.0 + A;       --  This is exact.
         A_From_1 := 1.0 - A;       --  Ditto.
         D := A_Plus_1 * A_From_1;  --  1 - A*A.

         --  use one term of the series expansion:

         --    f (x + e) = f(x) + e * f'(x) + ..

         --  The derivative of Arctanh at A is 1/(1-A*A). Next term is
         --  A*(B/D)**2 (if a quadratic approximation is ever needed).

         return 0.5 * (Log (A_Plus_1) - Log (A_From_1)) + B / D;
      end if;
   end Arctanh;

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

   --  Natural cycle

   function Cos (X : Float_Type'Base) return Float_Type'Base is
   begin
      if abs X < Sqrt_Epsilon then
         return 1.0;
      end if;

      return Float_Type'Base (Aux.Cos (Double (X)));
   end Cos;

   --  Arbitrary cycle

   function Cos (X, Cycle : Float_Type'Base) return Float_Type'Base is
   begin
      --  Just reuse the code for Sin. The potential small loss of speed is
      --  negligible with proper (front-end) inlining.

      return -Sin (abs X - Cycle * 0.25, Cycle);
   end Cos;

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

   function Cosh (X : Float_Type'Base) return Float_Type'Base is
      Lnv      : constant Float_Type'Base := 8#0.542714#;
      V2minus1 : constant Float_Type'Base := 0.13830_27787_96019_02638E-4;
      Y        : constant Float_Type'Base := abs X;
      Z        : Float_Type'Base;

   begin
      if Y < Sqrt_Epsilon then
         return 1.0;

      elsif Y > Log_Inverse_Epsilon then
         Z := Exp_Strict (Y - Lnv);
         return (Z + V2minus1 * Z);

      else
         Z := Exp_Strict (Y);
         return 0.5 * (Z + 1.0 / Z);
      end if;

   end Cosh;

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

   --  Natural cycle

   function Cot (X : Float_Type'Base) return Float_Type'Base is
   begin
      if X = 0.0 then
         raise Constraint_Error;

      elsif abs X < Sqrt_Epsilon then
         return 1.0 / X;
      end if;

      return 1.0 / Float_Type'Base (Aux.Tan (Double (X)));
   end Cot;

   --  Arbitrary cycle

   function Cot (X, Cycle : Float_Type'Base) return Float_Type'Base is
      T : Float_Type'Base;

   begin
      if Cycle <= 0.0 then
         raise Argument_Error;
      end if;

      T := Float_Type'Base'Remainder (X, Cycle);

      if T = 0.0 or else abs T = 0.5 * Cycle then
         raise Constraint_Error;

      elsif abs T < Sqrt_Epsilon then
         return 1.0 / T;

      elsif abs T = 0.25 * Cycle then
         return 0.0;

      else
         T := T / Cycle * Two_Pi;
         return Cos (T) / Sin (T);
      end if;
   end Cot;

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

   function Coth (X : Float_Type'Base) return Float_Type'Base is
   begin
      if X = 0.0 then
         raise Constraint_Error;

      elsif X < Half_Log_Epsilon then
         return -1.0;

      elsif X > -Half_Log_Epsilon then
         return 1.0;

      elsif abs X < Sqrt_Epsilon then
         return 1.0 / X;
      end if;

      return 1.0 / Float_Type'Base (Aux.Tanh (Double (X)));
   end Coth;

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

   function Exp (X : Float_Type'Base) return Float_Type'Base is
      Result : Float_Type'Base;

   begin
      if X = 0.0 then
         return 1.0;
      end if;

      Result := Float_Type'Base (Aux.Exp (Double (X)));

      --  Deal with case of Exp returning IEEE infinity. If Machine_Overflows
      --  is False, then we can just leave it as an infinity (and indeed we
      --  prefer to do so). But if Machine_Overflows is True, then we have
      --  to raise a Constraint_Error exception as required by the RM.

      if Float_Type'Machine_Overflows and then not Result'Valid then
         raise Constraint_Error;
      end if;

      return Result;
   end Exp;

   ----------------
   -- Exp_Strict --
   ----------------

   function Exp_Strict (X : Float_Type'Base) return Float_Type'Base is
      G : Float_Type'Base;
      Z : Float_Type'Base;

      P0 : constant := 0.25000_00000_00000_00000;
      P1 : constant := 0.75753_18015_94227_76666E-2;
      P2 : constant := 0.31555_19276_56846_46356E-4;

      Q0 : constant := 0.5;
      Q1 : constant := 0.56817_30269_85512_21787E-1;
      Q2 : constant := 0.63121_89437_43985_02557E-3;
      Q3 : constant := 0.75104_02839_98700_46114E-6;

      C1 : constant := 8#0.543#;
      C2 : constant := -2.1219_44400_54690_58277E-4;
      Le : constant := 1.4426_95040_88896_34074;

      XN : Float_Type'Base;
      P, Q, R : Float_Type'Base;

   begin
      if X = 0.0 then
         return 1.0;
      end if;

      XN := Float_Type'Base'Rounding (X * Le);
      G := (X - XN * C1) - XN * C2;
      Z := G * G;
      P := G * ((P2 * Z + P1) * Z + P0);
      Q := ((Q3 * Z + Q2) * Z + Q1) * Z + Q0;
      R := 0.5 + P / (Q - P);

      R := Float_Type'Base'Scaling (R, Integer (XN) + 1);

      --  Deal with case of Exp returning IEEE infinity. If Machine_Overflows
      --  is False, then we can just leave it as an infinity (and indeed we
      --  prefer to do so). But if Machine_Overflows is True, then we have to
      --  raise a Constraint_Error exception as required by the RM.

      if Float_Type'Machine_Overflows and then not R'Valid then
         raise Constraint_Error;
      else
         return R;
      end if;

   end Exp_Strict;

   ----------------
   -- Local_Atan --
   ----------------

   function Local_Atan
     (Y : Float_Type'Base;
      X : Float_Type'Base := 1.0) return Float_Type'Base
   is
      Z        : Float_Type'Base;
      Raw_Atan : Float_Type'Base;

   begin
      Z := (if abs Y > abs X then abs (X / Y) else abs (Y / X));

      Raw_Atan :=
        (if Z < Sqrt_Epsilon then Z
         elsif Z = 1.0 then Pi / 4.0
         else Float_Type'Base (Aux.Atan (Double (Z))));

      if abs Y > abs X then
         Raw_Atan := Half_Pi - Raw_Atan;
      end if;

      if X > 0.0 then
         return Float_Type'Copy_Sign (Raw_Atan, Y);
      else
         return Float_Type'Copy_Sign (Pi - Raw_Atan, Y);
      end if;
   end Local_Atan;

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

   --  Natural base

   function Log (X : Float_Type'Base) return Float_Type'Base is
   begin
      if X < 0.0 then
         raise Argument_Error;

      elsif X = 0.0 then
         raise Constraint_Error;

      elsif X = 1.0 then
         return 0.0;
      end if;

      return Float_Type'Base (Aux.Log (Double (X)));
   end Log;

   --  Arbitrary base

   function Log (X, Base : Float_Type'Base) return Float_Type'Base is
   begin
      if X < 0.0 then
         raise Argument_Error;

      elsif Base <= 0.0 or else Base = 1.0 then
         raise Argument_Error;

      elsif X = 0.0 then
         raise Constraint_Error;

      elsif X = 1.0 then
         return 0.0;
      end if;

      return Float_Type'Base (Aux.Log (Double (X)) / Aux.Log (Double (Base)));
   end Log;

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

   --  Natural cycle

   function Sin (X : Float_Type'Base) return Float_Type'Base is
   begin
      if abs X < Sqrt_Epsilon then
         return X;
      end if;

      return Float_Type'Base (Aux.Sin (Double (X)));
   end Sin;

   --  Arbitrary cycle

   function Sin (X, Cycle : Float_Type'Base) return Float_Type'Base is
      T : Float_Type'Base;

   begin
      if Cycle <= 0.0 then
         raise Argument_Error;

      --  If X is zero, return it as the result, preserving the argument sign.
      --  Is this test really needed on any machine ???

      elsif X = 0.0 then
         return X;
      end if;

      T := Float_Type'Base'Remainder (X, Cycle);

      --  The following two reductions reduce the argument to the interval
      --  [-0.25 * Cycle, 0.25 * Cycle]. This reduction is exact and is needed
      --  to prevent inaccuracy that may result if the sine function uses a
      --  different (more accurate) value of Pi in its reduction than is used
      --  in the multiplication with Two_Pi.

      if abs T > 0.25 * Cycle then
         T := 0.5 * Float_Type'Copy_Sign (Cycle, T) - T;
      end if;

      --  Could test for 12.0 * abs T = Cycle, and return an exact value in
      --  those cases. It is not clear this is worth the extra test though.

      return Float_Type'Base (Aux.Sin (Double (T / Cycle * Two_Pi)));
   end Sin;

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

   function Sinh (X : Float_Type'Base) return Float_Type'Base is
      Lnv      : constant Float_Type'Base := 8#0.542714#;
      V2minus1 : constant Float_Type'Base := 0.13830_27787_96019_02638E-4;
      Y        : constant Float_Type'Base := abs X;
      F        : constant Float_Type'Base := Y * Y;
      Z        : Float_Type'Base;

      Float_Digits_1_6 : constant Boolean := Float_Type'Digits < 7;

   begin
      if Y < Sqrt_Epsilon then
         return X;

      elsif Y > Log_Inverse_Epsilon then
         Z := Exp_Strict (Y - Lnv);
         Z := Z + V2minus1 * Z;

      elsif Y < 1.0 then

         if Float_Digits_1_6 then

            --  Use expansion provided by Cody and Waite, p. 226. Note that
            --  leading term of the polynomial in Q is exactly 1.0.

            declare
               P0 : constant := -0.71379_3159E+1;
               P1 : constant := -0.19033_3399E+0;
               Q0 : constant := -0.42827_7109E+2;

            begin
               Z := Y + Y * F * (P1 * F + P0) / (F + Q0);
            end;

         else
            declare
               P0 : constant := -0.35181_28343_01771_17881E+6;
               P1 : constant := -0.11563_52119_68517_68270E+5;
               P2 : constant := -0.16375_79820_26307_51372E+3;
               P3 : constant := -0.78966_12741_73570_99479E+0;
               Q0 : constant := -0.21108_77005_81062_71242E+7;
               Q1 : constant :=  0.36162_72310_94218_36460E+5;
               Q2 : constant := -0.27773_52311_96507_01667E+3;

            begin
               Z := Y + Y * F * (((P3 * F + P2) * F + P1) * F + P0)
                              / (((F + Q2) * F + Q1) * F + Q0);
            end;
         end if;

      else
         Z := Exp_Strict (Y);
         Z := 0.5 * (Z - 1.0 / Z);
      end if;

      if X > 0.0 then
         return Z;
      else
         return -Z;
      end if;
   end Sinh;

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

   function Sqrt (X : Float_Type'Base) return Float_Type'Base is
   begin
      if X < 0.0 then
         raise Argument_Error;

      --  Special case Sqrt (0.0) to preserve possible minus sign per IEEE

      elsif X = 0.0 then
         return X;
      end if;

      return Float_Type'Base (Aux.Sqrt (Double (X)));
   end Sqrt;

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

   --  Natural cycle

   function Tan (X : Float_Type'Base) return Float_Type'Base is
   begin
      if abs X < Sqrt_Epsilon then
         return X;
      end if;

      --  Note: if X is exactly pi/2, then we should raise an exception, since
      --  the result would overflow. But for all floating-point formats we deal
      --  with, it is impossible for X to be exactly pi/2, and the result is
      --  always in range.

      return Float_Type'Base (Aux.Tan (Double (X)));
   end Tan;

   --  Arbitrary cycle

   function Tan (X, Cycle : Float_Type'Base) return Float_Type'Base is
      T : Float_Type'Base;

   begin
      if Cycle <= 0.0 then
         raise Argument_Error;

      elsif X = 0.0 then
         return X;
      end if;

      T := Float_Type'Base'Remainder (X, Cycle);

      if abs T = 0.25 * Cycle then
         raise Constraint_Error;

      elsif abs T = 0.5 * Cycle then
         return 0.0;

      else
         T := T / Cycle * Two_Pi;
         return Sin (T) / Cos (T);
      end if;

   end Tan;

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

   function Tanh (X : Float_Type'Base) return Float_Type'Base is
      P0 : constant Float_Type'Base := -0.16134_11902_39962_28053E+4;
      P1 : constant Float_Type'Base := -0.99225_92967_22360_83313E+2;
      P2 : constant Float_Type'Base := -0.96437_49277_72254_69787E+0;

      Q0 : constant Float_Type'Base :=  0.48402_35707_19886_88686E+4;
      Q1 : constant Float_Type'Base :=  0.22337_72071_89623_12926E+4;
      Q2 : constant Float_Type'Base :=  0.11274_47438_05349_49335E+3;
      Q3 : constant Float_Type'Base :=  0.10000_00000_00000_00000E+1;

      Half_Ln3 : constant Float_Type'Base := 0.54930_61443_34054_84570;

      P, Q, R : Float_Type'Base;
      Y : constant Float_Type'Base := abs X;
      G : constant Float_Type'Base := Y * Y;

      Float_Type_Digits_15_Or_More : constant Boolean :=
        Float_Type'Digits > 14;

   begin
      if X < Half_Log_Epsilon then
         return -1.0;

      elsif X > -Half_Log_Epsilon then
         return 1.0;

      elsif Y < Sqrt_Epsilon then
         return X;

      elsif Y < Half_Ln3
        and then Float_Type_Digits_15_Or_More
      then
         P := (P2 * G + P1) * G + P0;
         Q := ((Q3 * G + Q2) * G + Q1) * G + Q0;
         R := G * (P / Q);
         return X + X * R;

      else
         return Float_Type'Base (Aux.Tanh (Double (X)));
      end if;
   end Tanh;

end Ada.Numerics.Generic_Elementary_Functions;