------------------------------------------------------------------------------
-- --
--- GNAT RUNTIME COMPONENTS --
+-- GNAT RUN-TIME COMPONENTS --
-- --
-- ADA.NUMERICS.GENERIC_ELEMENTARY_FUNCTIONS --
-- --
-- B o d y --
-- --
--- --
--- Copyright (C) 1992-2001, Free Software Foundation, Inc. --
+-- Copyright (C) 1992-2009, 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- --
+-- 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. 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, 59 Temple Place - Suite 330, Boston, --
--- MA 02111-1307, USA. --
+-- 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. --
-- --
--- 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. --
+-- 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. --
--- It is now maintained by Ada Core Technologies Inc (http://www.gnat.com). --
+-- Extensive contributions were provided by Ada Core Technologies Inc. --
-- --
------------------------------------------------------------------------------
-- advantage of the C functions, e.g. in providing interface to hardware
-- provided versions of the elementary functions.
--- Uses functions sqrt, exp, log, pow, sin, asin, cos, acos, tan, atan,
--- sinh, cosh, tanh from C library via math.h
+-- 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;
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;
- Fourth_Pi : constant T := Pi / 4.0;
+ Two_Pi : constant T := 2.0 * Pi;
+ Half_Pi : constant T := Pi / 2.0;
- Epsilon : constant T := 2.0 ** (1 - T'Model_Mantissa);
- IEpsilon : constant T := 2.0 ** (T'Model_Mantissa - 1);
- Log_Epsilon : constant T := T (1 - T'Model_Mantissa) * Log_Two;
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);
- DEpsilon : constant Double := Double (Epsilon);
- DIEpsilon : constant Double := Double (IEpsilon);
-
-----------------------
-- 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.
+ -- 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 cyele reduction
+ (Y : Float_Type'Base;
+ X : Float_Type'Base := 1.0) return Float_Type'Base;
+ -- Common code for arc tangent after cycle reduction
----------
-- "**" --
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.
+ -- 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)
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.
+ -- 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;
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.
+ -- 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;
return Pi / 2.0;
elsif X = -1.0 then
- return -Pi / 2.0;
+ return -(Pi / 2.0);
end if;
return Float_Type'Base (Aux.Asin (Double (X)));
return Cycle / 4.0;
elsif X = -1.0 then
- return -Cycle / 4.0;
+ return -(Cycle / 4.0);
end if;
return Arctan (X / Sqrt ((1.0 - X) * (1.0 + X)), 1.0, Cycle);
elsif X > 1.0 / Sqrt_Epsilon then
return Log (X) + Log_Two;
- elsif X < -1.0 / Sqrt_Epsilon then
+ elsif X < -(1.0 / Sqrt_Epsilon) then
return -(Log (-X) + Log_Two);
elsif X < 0.0 then
return Float_Type'Base
is
begin
- if X = 0.0
- and then Y = 0.0
- then
+ if X = 0.0 and then Y = 0.0 then
raise Argument_Error;
elsif Y = 0.0 then
end if;
elsif X = 0.0 then
- if Y > 0.0 then
- return Half_Pi;
- else -- Y < 0.0
- return -Half_Pi;
- end if;
+ return Float_Type'Copy_Sign (Half_Pi, Y);
else
return Local_Atan (Y, X);
if Cycle <= 0.0 then
raise Argument_Error;
- elsif X = 0.0
- and then Y = 0.0
- then
+ elsif X = 0.0 and then Y = 0.0 then
raise Argument_Error;
elsif Y = 0.0 then
end if;
elsif X = 0.0 then
- if Y > 0.0 then
- return Cycle / 4.0;
- else -- Y < 0.0
- return -Cycle / 4.0;
- end if;
+ return Float_Type'Copy_Sign (Cycle / 4.0, Y);
else
return Local_Atan (Y, X) * Cycle / Two_Pi;
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
-- 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.
+ -- 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
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) + ..
+
+ -- 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;
-
- -- else
- -- return 0.5 * Log ((X + 1.0) / (1.0 - X));
- -- why are above lines commented out ???
end if;
end Arctanh;
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.
+ -- 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;
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 : Float_Type'Base := abs X;
+ Y : constant Float_Type'Base := abs X;
Z : Float_Type'Base;
begin
T := Float_Type'Base'Remainder (X, Cycle);
- if T = 0.0 or abs T = 0.5 * Cycle then
+ if T = 0.0 or else abs T = 0.5 * Cycle then
raise Constraint_Error;
elsif abs T < Sqrt_Epsilon then
else
T := T / Cycle * Two_Pi;
- return Cos (T) / Sin (T);
+ return Cos (T) / Sin (T);
end if;
end Cot;
-- 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.
+ -- 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;
----------------
function Local_Atan
- (Y : Float_Type'Base;
- X : Float_Type'Base := 1.0)
- return Float_Type'Base
+ (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
- if abs Y > abs X then
- Z := abs (X / Y);
- else
- Z := abs (Y / X);
- end if;
-
- if Z < Sqrt_Epsilon then
- Raw_Atan := Z;
+ Z := (if abs Y > abs X then abs (X / Y) else abs (Y / X));
- elsif Z = 1.0 then
- Raw_Atan := Pi / 4.0;
-
- else
- Raw_Atan := Float_Type'Base (Aux.Atan (Double (Z)));
- end if;
+ 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
- if Y > 0.0 then
- return Raw_Atan;
- else -- Y < 0.0
- return -Raw_Atan;
- end if;
-
- else -- X < 0.0
- if Y > 0.0 then
- return Pi - Raw_Atan;
- else -- Y < 0.0
- return -(Pi - Raw_Atan);
- end if;
+ return Float_Type'Copy_Sign (Raw_Atan, Y);
+ else
+ return Float_Type'Copy_Sign (Pi - Raw_Atan, Y);
end if;
end Local_Atan;
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
- -- Is this test really needed on any machine ???
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 sinus function
- -- a different (more accurate) value of Pi in its
- -- reduction than is used in the multiplication with Two_Pi.
+ -- 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 that
- -- this is worth the extra test though.
+ -- 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)));
+ return Float_Type'Base (Aux.Sin (Double (T / Cycle * Two_Pi)));
end Sin;
----------
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 : Float_Type'Base := abs X;
+ Y : constant Float_Type'Base := abs X;
F : constant Float_Type'Base := Y * Y;
Z : Float_Type'Base;
elsif X = 0.0 then
return X;
-
end if;
return Float_Type'Base (Aux.Sqrt (Double (X)));
----------
function Tanh (X : Float_Type'Base) return Float_Type'Base is
- P0 : constant Float_Type'Base := -0.16134_11902E4;
- P1 : constant Float_Type'Base := -0.99225_92967E2;
- P2 : constant Float_Type'Base := -0.96437_49299E0;
+ 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_35707E4;
- Q1 : constant Float_Type'Base := 0.22337_72071E4;
- Q2 : constant Float_Type'Base := 0.11274_47438E3;
- Q3 : constant Float_Type'Base := 0.10000000000E1;
+ 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;
+ Half_Ln3 : constant Float_Type'Base := 0.54930_61443_34054_84570;
P, Q, R : Float_Type'Base;
- Y : Float_Type'Base := abs X;
- G : Float_Type'Base := Y * Y;
+ 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;
OSDN Git Service
