math_lib.adb   [plain text]

------------------------------------------------------------------------------
--                                                                          --
--                         GNAT RUN-TIME COMPONENTS                         --
--                                                                          --
--                             M A T H _ L I B                              --
--                                                                          --
--                                 B o d y                                  --
--                                                                          --
--          Copyright (C) 1992-2005 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 2,  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.  See the GNU General Public License --
-- for  more details.  You should have  received  a copy of the GNU General --
-- Public License  distributed with GNAT;  see file COPYING.  If not, write --
-- to  the  Free Software Foundation,  51  Franklin  Street,  Fifth  Floor, --
-- Boston, MA 02110-1301, USA.                                              --
--                                                                          --
-- As a special exception,  if other files  instantiate  generics from this --
-- unit, or you link  this unit with other files  to produce an executable, --
-- this  unit  does not  by itself cause  the resulting  executable  to  be --
-- covered  by the  GNU  General  Public  License.  This exception does not --
-- however invalidate  any other reasons why  the executable file  might be --
-- covered by the  GNU Public License.                                      --
--                                                                          --
-- 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. This is not a "strict" implementation, but takes full
--  advantage of the C functions, e.g. in providing interface to hardware
--  provided versions of the elementary functions.

--  A known weakness is that on the x86, all computation is done in Double,
--  which means that a lot of accuracy is lost for the Long_Long_Float case.

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

--  This is an adaptation of Ada.Numerics.Generic_Elementary_Functions that
--  provides a compatible body for the DEC Math_Lib package.

with Ada.Numerics.Aux;
use type Ada.Numerics.Aux.Double;
with Ada.Numerics; use Ada.Numerics;

package body Math_Lib is

   Log_Two : constant := 0.69314_71805_59945_30941_72321_21458_17656_80755;

   Two_Pi     : constant Real'Base := 2.0 * Pi;
   Half_Pi    : constant Real'Base := Pi / 2.0;
   Fourth_Pi  : constant Real'Base := Pi / 4.0;
   Epsilon    : constant Real'Base := Real'Base'Epsilon;
   IEpsilon   : constant Real'Base := 1.0 / Epsilon;

   subtype Double is Aux.Double;

   DEpsilon    : constant Double := Double (Epsilon);
   DIEpsilon   : constant Double := Double (IEpsilon);

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

   function Arctan
     (Y    : Real;
      A    : Real := 1.0)
      return Real;

   function Arctan
     (Y     : Real;
      A     : Real := 1.0;
      Cycle : Real)
      return  Real;

   function Exact_Remainder
     (A    : Real;
      Y    : Real)
      return Real;
   --  Computes exact remainder of A divided by Y

   function Half_Log_Epsilon return Real;
   --  Function to provide constant: 0.5 * Log (Epsilon)

   function Local_Atan
     (Y    : Real;
      A    : Real := 1.0)
      return Real;
   --  Common code for arc tangent after cyele reduction

   function Log_Inverse_Epsilon return Real;
   --  Function to provide constant: Log (1.0 / Epsilon)

   function Square_Root_Epsilon return Real;
   --  Function to provide constant: Sqrt (Epsilon)

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

   function "**" (A1, A2 : Real) return Real is

   begin
      if A1 = 0.0
        and then A2 = 0.0
      then
         raise Argument_Error;

      elsif A1 < 0.0 then
         raise Argument_Error;

      elsif A2 = 0.0 then
         return 1.0;

      elsif A1 = 0.0 then
         if A2 < 0.0 then
            raise Constraint_Error;
         else
            return 0.0;
         end if;

      elsif A1 = 1.0 then
         return 1.0;

      elsif A2 = 1.0 then
         return A1;

      else
         begin
            if A2 = 2.0 then
               return A1 * A1;
            else
               return
                 Real (Aux.pow (Double (A1), Double (A2)));
            end if;

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

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

   --  Natural cycle

   function Arccos (A : Real) return Real is
      Temp : Real'Base;

   begin
      if abs A > 1.0 then
         raise Argument_Error;

      elsif abs A < Square_Root_Epsilon then
         return Pi / 2.0 - A;

      elsif A = 1.0 then
         return 0.0;

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

      Temp := Real (Aux.acos (Double (A)));

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

      return Temp;
   end Arccos;

   --  Arbitrary cycle

   function Arccos (A, Cycle : Real) return Real is
      Temp : Real'Base;

   begin
      if Cycle <= 0.0 then
         raise Argument_Error;

      elsif abs A > 1.0 then
         raise Argument_Error;

      elsif abs A < Square_Root_Epsilon then
         return Cycle / 4.0;

      elsif A = 1.0 then
         return 0.0;

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

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

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

      return Temp;
   end Arccos;

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

   function Arccosh (A : Real) return Real is
   begin
      --  Return Log (A - Sqrt (A * A - 1.0));  double valued,
      --    only positive value returned
      --  What is this comment ???

      if A < 1.0 then
         raise Argument_Error;

      elsif A < 1.0 + Square_Root_Epsilon then
         return A - 1.0;

      elsif abs A > 1.0 / Square_Root_Epsilon then
         return Log (A) + Log_Two;

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

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

   --  Natural cycle

   function Arccot
     (A    : Real;
      Y    : Real := 1.0)
      return Real
   is
   begin
      --  Just reverse arguments

      return Arctan (Y, A);
   end Arccot;

   --  Arbitrary cycle

   function Arccot
     (A     : Real;
      Y     : Real := 1.0;
      Cycle : Real)
      return  Real
   is
   begin
      --  Just reverse arguments

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

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

   function Arccoth (A : Real) return Real is
   begin
      if abs A = 1.0 then
         raise Constraint_Error;

      elsif abs A < 1.0 then
         raise Argument_Error;

      elsif abs A > 1.0 / Epsilon then
         return 0.0;

      else
         return 0.5 * Log ((1.0 + A) / (A - 1.0));
      end if;
   end Arccoth;

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

   --  Natural cycle

   function Arcsin (A : Real) return Real is
   begin
      if abs A > 1.0 then
         raise Argument_Error;

      elsif abs A < Square_Root_Epsilon then
         return A;

      elsif A = 1.0 then
         return Pi / 2.0;

      elsif A = -1.0 then
         return -Pi / 2.0;
      end if;

      return Real (Aux.asin (Double (A)));
   end Arcsin;

   --  Arbitrary cycle

   function Arcsin (A, Cycle : Real) return Real is
   begin
      if Cycle <= 0.0 then
         raise Argument_Error;

      elsif abs A > 1.0 then
         raise Argument_Error;

      elsif A = 0.0 then
         return A;

      elsif A = 1.0 then
         return Cycle / 4.0;

      elsif A = -1.0 then
         return -Cycle / 4.0;
      end if;

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

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

   function Arcsinh (A : Real) return Real is
   begin
      if abs A < Square_Root_Epsilon then
         return A;

      elsif A > 1.0 / Square_Root_Epsilon then
         return Log (A) + Log_Two;

      elsif A < -1.0 / Square_Root_Epsilon then
         return -(Log (-A) + Log_Two);

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

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

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

   --  Natural cycle

   function Arctan
     (Y    : Real;
      A    : Real := 1.0)
      return Real
   is
   begin
      if A = 0.0
        and then Y = 0.0
      then
         raise Argument_Error;

      elsif Y = 0.0 then
         if A > 0.0 then
            return 0.0;
         else -- A < 0.0
            return Pi;
         end if;

      elsif A = 0.0 then
         if Y > 0.0 then
            return Half_Pi;
         else -- Y < 0.0
            return -Half_Pi;
         end if;

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

   --  Arbitrary cycle

   function Arctan
     (Y     : Real;
      A     : Real := 1.0;
      Cycle : Real)
      return  Real
   is
   begin
      if Cycle <= 0.0 then
         raise Argument_Error;

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

      elsif Y = 0.0 then
         if A > 0.0 then
            return 0.0;
         else -- A < 0.0
            return Cycle / 2.0;
         end if;

      elsif A = 0.0 then
         if Y > 0.0 then
            return Cycle / 4.0;
         else -- Y < 0.0
            return -Cycle / 4.0;
         end if;

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

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

   function Arctanh (A : Real) return Real is
   begin
      if abs A = 1.0 then
         raise Constraint_Error;

      elsif abs A > 1.0 then
         raise Argument_Error;

      elsif abs A < Square_Root_Epsilon then
         return A;

      else
         return 0.5 * Log ((1.0 + A) / (1.0 - A));
      end if;
   end Arctanh;

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

   --  Natural cycle

   function Cos (A : Real) return Real is
   begin
      if A = 0.0 then
         return 1.0;

      elsif abs A < Square_Root_Epsilon then
         return 1.0;

      end if;

      return Real (Aux.Cos (Double (A)));
   end Cos;

   --  Arbitrary cycle

   function Cos (A, Cycle : Real) return Real is
      T : Real'Base;

   begin
      if Cycle <= 0.0 then
         raise Argument_Error;

      elsif A = 0.0 then
         return 1.0;
      end if;

      T := Exact_Remainder (abs (A), Cycle) / Cycle;

      if T = 0.25
        or else T = 0.75
        or else T = -0.25
        or else T = -0.75
      then
         return 0.0;

      elsif T = 0.5 or T = -0.5 then
         return -1.0;
      end if;

      return Real (Aux.Cos (Double (T * Two_Pi)));
   end Cos;

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

   function Cosh (A : Real) return Real is
   begin
      if abs A < Square_Root_Epsilon then
         return 1.0;

      elsif abs A > Log_Inverse_Epsilon then
         return Exp ((abs A) - Log_Two);
      end if;

      return Real (Aux.cosh (Double (A)));

   exception
      when others =>
         raise Constraint_Error;
   end Cosh;

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

   --  Natural cycle

   function Cot (A : Real) return Real is
   begin
      if A = 0.0 then
         raise Constraint_Error;

      elsif abs A < Square_Root_Epsilon then
         return 1.0 / A;
      end if;

      return Real (1.0 / Real'Base (Aux.tan (Double (A))));
   end Cot;

   --  Arbitrary cycle

   function Cot (A, Cycle : Real) return Real is
      T : Real'Base;

   begin
      if Cycle <= 0.0 then
         raise Argument_Error;

      elsif A = 0.0 then
         raise Constraint_Error;

      elsif abs A < Square_Root_Epsilon then
         return 1.0 / A;
      end if;

      T := Exact_Remainder (A, Cycle) / Cycle;

      if T = 0.0 or T = 0.5 or T = -0.5 then
         raise Constraint_Error;
      else
         return  Cos (T * Two_Pi) / Sin (T * Two_Pi);
      end if;
   end Cot;

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

   function Coth (A : Real) return Real is
   begin
      if A = 0.0 then
         raise Constraint_Error;

      elsif A < Half_Log_Epsilon then
         return -1.0;

      elsif A > -Half_Log_Epsilon then
         return 1.0;

      elsif abs A < Square_Root_Epsilon then
         return 1.0 / A;
      end if;

      return Real (1.0 / Real'Base (Aux.tanh (Double (A))));
   end Coth;

   ---------------------
   -- Exact_Remainder --
   ---------------------

   function Exact_Remainder
     (A    : Real;
      Y    : Real)
      return Real
   is
      Denominator : Real'Base := abs A;
      Divisor     : Real'Base := abs Y;
      Reducer     : Real'Base;
      Sign        : Real'Base := 1.0;

   begin
      if Y = 0.0 then
         raise Constraint_Error;

      elsif A = 0.0 then
         return 0.0;

      elsif A = Y then
         return 0.0;

      elsif Denominator < Divisor then
         return A;
      end if;

      while Denominator >= Divisor loop

         --  Put divisors mantissa with denominators exponent to make reducer

         Reducer := Divisor;

         begin
            while Reducer * 1_048_576.0 < Denominator loop
               Reducer := Reducer * 1_048_576.0;
            end loop;

         exception
            when others => null;
         end;

         begin
            while Reducer * 1_024.0 < Denominator loop
               Reducer := Reducer * 1_024.0;
            end loop;

         exception
            when others => null;
         end;

         begin
            while Reducer * 2.0 < Denominator loop
               Reducer := Reducer * 2.0;
            end loop;

         exception
            when others => null;
         end;

         Denominator := Denominator - Reducer;
      end loop;

      if A < 0.0 then
         return -Denominator;
      else
         return Denominator;
      end if;
   end Exact_Remainder;

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

   function Exp (A : Real) return Real is
      Result : Real'Base;

   begin
      if A = 0.0 then
         return 1.0;

      else
         Result := Real (Aux.Exp (Double (A)));

         --  The check here catches the case of Exp returning IEEE infinity

         if Result > Real'Last then
            raise Constraint_Error;
         else
            return Result;
         end if;
      end if;
   end Exp;

   ----------------------
   -- Half_Log_Epsilon --
   ----------------------

   --  Cannot precompute this constant, because this is required to be a
   --  pure package, which allows no state. A pity, but no way around it!

   function Half_Log_Epsilon return Real is
   begin
      return Real (0.5 * Real'Base (Aux.Log (DEpsilon)));
   end Half_Log_Epsilon;

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

   function Local_Atan
     (Y    : Real;
      A    : Real := 1.0)
      return Real
   is
      Z        : Real'Base;
      Raw_Atan : Real'Base;

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

      if Z < Square_Root_Epsilon then
         Raw_Atan := Z;

      elsif Z = 1.0 then
         Raw_Atan := Pi / 4.0;

      elsif Z < Square_Root_Epsilon then
         Raw_Atan := Z;

      else
         Raw_Atan := Real'Base (Aux.Atan (Double (Z)));
      end if;

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

      if A > 0.0 then
         if Y > 0.0 then
            return Raw_Atan;
         else                 --  Y < 0.0
            return -Raw_Atan;
         end if;

      else                    --  A < 0.0
         if Y > 0.0 then
            return Pi - Raw_Atan;
         else                  --  Y < 0.0
            return -(Pi - Raw_Atan);
         end if;
      end if;
   end Local_Atan;

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

   --  Natural base

   function Log (A : Real) return Real is
   begin
      if A < 0.0 then
         raise Argument_Error;

      elsif A = 0.0 then
         raise Constraint_Error;

      elsif A = 1.0 then
         return 0.0;
      end if;

      return Real (Aux.Log (Double (A)));
   end Log;

   --  Arbitrary base

   function Log (A, Base : Real) return Real is
   begin
      if A < 0.0 then
         raise Argument_Error;

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

      elsif A = 0.0 then
         raise Constraint_Error;

      elsif A = 1.0 then
         return 0.0;
      end if;

      return Real (Aux.Log (Double (A)) / Aux.Log (Double (Base)));
   end Log;

   -------------------------
   -- Log_Inverse_Epsilon --
   -------------------------

   --  Cannot precompute this constant, because this is required to be a
   --  pure package, which allows no state. A pity, but no way around it!

   function Log_Inverse_Epsilon return Real is
   begin
      return Real (Aux.Log (DIEpsilon));
   end Log_Inverse_Epsilon;

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

   --  Natural cycle

   function Sin (A : Real) return Real is
   begin
      if abs A < Square_Root_Epsilon then
         return A;
      end if;

      return Real (Aux.Sin (Double (A)));
   end Sin;

   --  Arbitrary cycle

   function Sin (A, Cycle : Real) return Real is
      T : Real'Base;

   begin
      if Cycle <= 0.0 then
         raise Argument_Error;

      elsif A = 0.0 then
         return A;
      end if;

      T := Exact_Remainder (A, Cycle) / Cycle;

      if T = 0.0 or T = 0.5 or T = -0.5 then
         return 0.0;

      elsif T = 0.25 or T = -0.75 then
         return 1.0;

      elsif T = -0.25 or T = 0.75 then
         return -1.0;

      end if;

      return Real (Aux.Sin (Double (T * Two_Pi)));
   end Sin;

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

   function Sinh (A : Real) return Real is
   begin
      if abs A < Square_Root_Epsilon then
         return A;

      elsif  A > Log_Inverse_Epsilon then
         return Exp (A - Log_Two);

      elsif A < -Log_Inverse_Epsilon then
         return -Exp ((-A) - Log_Two);
      end if;

      return Real (Aux.Sinh (Double (A)));

   exception
      when others =>
         raise Constraint_Error;
   end Sinh;

   -------------------------
   -- Square_Root_Epsilon --
   -------------------------

   --  Cannot precompute this constant, because this is required to be a
   --  pure package, which allows no state. A pity, but no way around it!

   function Square_Root_Epsilon return Real is
   begin
      return Real (Aux.Sqrt (DEpsilon));
   end Square_Root_Epsilon;

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

   function Sqrt (A : Real) return Real is
   begin
      if A < 0.0 then
         raise Argument_Error;

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

      elsif A = 0.0 then
         return A;

      --  Sqrt (1.0) must be exact for good complex accuracy

      elsif A = 1.0 then
         return 1.0;

      end if;

      return Real (Aux.Sqrt (Double (A)));
   end Sqrt;

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

   --  Natural cycle

   function Tan (A : Real) return Real is
   begin
      if abs A < Square_Root_Epsilon then
         return A;

      elsif abs A = Pi / 2.0 then
         raise Constraint_Error;
      end if;

      return Real (Aux.tan (Double (A)));
   end Tan;

   --  Arbitrary cycle

   function Tan (A, Cycle : Real) return Real is
      T : Real'Base;

   begin
      if Cycle <= 0.0 then
         raise Argument_Error;

      elsif A = 0.0 then
         return A;
      end if;

      T := Exact_Remainder (A, Cycle) / Cycle;

      if T = 0.25
        or else T = 0.75
        or else T = -0.25
        or else T = -0.75
      then
         raise Constraint_Error;

      else
         return  Sin (T * Two_Pi) / Cos (T * Two_Pi);
      end if;
   end Tan;

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

   function Tanh (A : Real) return Real is
   begin
      if A < Half_Log_Epsilon then
         return -1.0;

      elsif A > -Half_Log_Epsilon then
         return 1.0;

      elsif abs A < Square_Root_Epsilon then
         return A;
      end if;

      return Real (Aux.tanh (Double (A)));
   end Tanh;

   ----------------------------
   -- DEC-Specific functions --
   ----------------------------

   function LOG10  (A : REAL) return REAL is
   begin
      return Log (A, 10.0);
   end LOG10;

   function LOG2   (A : REAL) return REAL is
   begin
      return Log (A, 2.0);
   end LOG2;

   function ASIN (A : REAL) return REAL renames Arcsin;
   function ACOS (A : REAL) return REAL renames Arccos;

   function ATAN (A : REAL) return REAL is
   begin
      return Arctan (A, 1.0);
   end ATAN;

   function ATAN2 (A1, A2 : REAL) return REAL renames Arctan;

   function SIND   (A : REAL) return REAL is
   begin
      return Sin (A, 360.0);
   end SIND;

   function COSD   (A : REAL) return REAL is
   begin
      return  Cos (A, 360.0);
   end COSD;

   function TAND   (A : REAL) return REAL is
   begin
      return  Tan (A, 360.0);
   end TAND;

   function ASIND  (A : REAL) return REAL is
   begin
      return  Arcsin (A, 360.0);
   end ASIND;

   function ACOSD  (A : REAL) return REAL is
   begin
      return  Arccos (A, 360.0);
   end ACOSD;

   function Arctan  (A : REAL) return REAL is
   begin
      return  Arctan (A, 1.0, 360.0);
   end Arctan;

   function ATAND (A : REAL) return REAL is
   begin
      return Arctan (A, 1.0, 360.0);
   end ATAND;

   function ATAN2D (A1, A2 : REAL) return REAL is
   begin
      return Arctan (A1, A2, 360.0);
   end ATAN2D;

end Math_Lib;