aboutsummaryrefslogtreecommitdiffstats
path: root/gcc-4.4.3/gcc/ada/a-numaux-darwin.adb
diff options
context:
space:
mode:
Diffstat (limited to 'gcc-4.4.3/gcc/ada/a-numaux-darwin.adb')
-rw-r--r--gcc-4.4.3/gcc/ada/a-numaux-darwin.adb185
1 files changed, 185 insertions, 0 deletions
diff --git a/gcc-4.4.3/gcc/ada/a-numaux-darwin.adb b/gcc-4.4.3/gcc/ada/a-numaux-darwin.adb
new file mode 100644
index 000000000..1444603d6
--- /dev/null
+++ b/gcc-4.4.3/gcc/ada/a-numaux-darwin.adb
@@ -0,0 +1,185 @@
+------------------------------------------------------------------------------
+-- --
+-- GNAT RUN-TIME COMPONENTS --
+-- --
+-- A D A . N U M E R I C S . A U X --
+-- --
+-- B o d y --
+-- (Apple OS X Version) --
+-- --
+-- Copyright (C) 1998-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 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. --
+-- --
+------------------------------------------------------------------------------
+
+-- File a-numaux.adb <- a-numaux-darwin.adb
+
+package body Ada.Numerics.Aux is
+
+ -----------------------
+ -- Local subprograms --
+ -----------------------
+
+ procedure Reduce (X : in out Double; Q : out Natural);
+ -- Implements reduction of X by Pi/2. Q is the quadrant of the final
+ -- result in the range 0 .. 3. The absolute value of X is at most Pi/4.
+
+ -- The following three functions implement Chebishev approximations
+ -- of the trigonometric functions in their reduced domain.
+ -- These approximations have been computed using Maple.
+
+ function Sine_Approx (X : Double) return Double;
+ function Cosine_Approx (X : Double) return Double;
+
+ pragma Inline (Reduce);
+ pragma Inline (Sine_Approx);
+ pragma Inline (Cosine_Approx);
+
+ function Cosine_Approx (X : Double) return Double is
+ XX : constant Double := X * X;
+ begin
+ return (((((16#8.DC57FBD05F640#E-08 * XX
+ - 16#4.9F7D00BF25D80#E-06) * XX
+ + 16#1.A019F7FDEFCC2#E-04) * XX
+ - 16#5.B05B058F18B20#E-03) * XX
+ + 16#A.AAAAAAAA73FA8#E-02) * XX
+ - 16#7.FFFFFFFFFFDE4#E-01) * XX
+ - 16#3.655E64869ECCE#E-14 + 1.0;
+ end Cosine_Approx;
+
+ function Sine_Approx (X : Double) return Double is
+ XX : constant Double := X * X;
+ begin
+ return (((((16#A.EA2D4ABE41808#E-09 * XX
+ - 16#6.B974C10F9D078#E-07) * XX
+ + 16#2.E3BC673425B0E#E-05) * XX
+ - 16#D.00D00CCA7AF00#E-04) * XX
+ + 16#2.222222221B190#E-02) * XX
+ - 16#2.AAAAAAAAAAA44#E-01) * (XX * X) + X;
+ end Sine_Approx;
+
+ ------------
+ -- Reduce --
+ ------------
+
+ procedure Reduce (X : in out Double; Q : out Natural) is
+ Half_Pi : constant := Pi / 2.0;
+ Two_Over_Pi : constant := 2.0 / Pi;
+
+ HM : constant := Integer'Min (Double'Machine_Mantissa / 2, Natural'Size);
+ M : constant Double := 0.5 + 2.0**(1 - HM); -- Splitting constant
+ P1 : constant Double := Double'Leading_Part (Half_Pi, HM);
+ P2 : constant Double := Double'Leading_Part (Half_Pi - P1, HM);
+ P3 : constant Double := Double'Leading_Part (Half_Pi - P1 - P2, HM);
+ P4 : constant Double := Double'Leading_Part (Half_Pi - P1 - P2 - P3, HM);
+ P5 : constant Double := Double'Leading_Part (Half_Pi - P1 - P2 - P3
+ - P4, HM);
+ P6 : constant Double := Double'Model (Half_Pi - P1 - P2 - P3 - P4 - P5);
+ K : Double;
+
+ begin
+ -- For X < 2.0**HM, all products below are computed exactly.
+ -- Due to cancellation effects all subtractions are exact as well.
+ -- As no double extended floating-point number has more than 75
+ -- zeros after the binary point, the result will be the correctly
+ -- rounded result of X - K * (Pi / 2.0).
+
+ K := X * Two_Over_Pi;
+ while abs K >= 2.0 ** HM loop
+ K := K * M - (K * M - K);
+ X :=
+ (((((X - K * P1) - K * P2) - K * P3) - K * P4) - K * P5) - K * P6;
+ K := X * Two_Over_Pi;
+ end loop;
+
+ -- If K is not a number (because X was not finite) raise exception
+
+ if K /= K then
+ raise Constraint_Error;
+ end if;
+
+ K := Double'Rounding (K);
+ Q := Integer (K) mod 4;
+ X := (((((X - K * P1) - K * P2) - K * P3)
+ - K * P4) - K * P5) - K * P6;
+ end Reduce;
+
+ ---------
+ -- Cos --
+ ---------
+
+ function Cos (X : Double) return Double is
+ Reduced_X : Double := abs X;
+ Quadrant : Natural range 0 .. 3;
+
+ begin
+ if Reduced_X > Pi / 4.0 then
+ Reduce (Reduced_X, Quadrant);
+
+ case Quadrant is
+ when 0 =>
+ return Cosine_Approx (Reduced_X);
+
+ when 1 =>
+ return Sine_Approx (-Reduced_X);
+
+ when 2 =>
+ return -Cosine_Approx (Reduced_X);
+
+ when 3 =>
+ return Sine_Approx (Reduced_X);
+ end case;
+ end if;
+
+ return Cosine_Approx (Reduced_X);
+ end Cos;
+
+ ---------
+ -- Sin --
+ ---------
+
+ function Sin (X : Double) return Double is
+ Reduced_X : Double := X;
+ Quadrant : Natural range 0 .. 3;
+
+ begin
+ if abs X > Pi / 4.0 then
+ Reduce (Reduced_X, Quadrant);
+
+ case Quadrant is
+ when 0 =>
+ return Sine_Approx (Reduced_X);
+
+ when 1 =>
+ return Cosine_Approx (Reduced_X);
+
+ when 2 =>
+ return Sine_Approx (-Reduced_X);
+
+ when 3 =>
+ return -Cosine_Approx (Reduced_X);
+ end case;
+ end if;
+
+ return Sine_Approx (Reduced_X);
+ end Sin;
+
+end Ada.Numerics.Aux;