1------------------------------------------------------------------------------
2--                                                                          --
3--                         GNAT RUN-TIME COMPONENTS                         --
4--                                                                          --
5--                     A D A . N U M E R I C S . A U X                      --
6--                                                                          --
7--                                 B o d y                                  --
8--                          (Apple OS X Version)                            --
9--                                                                          --
10--          Copyright (C) 1998-2009, Free Software Foundation, Inc.         --
11--                                                                          --
12-- GNAT is free software;  you can  redistribute it  and/or modify it under --
13-- terms of the  GNU General Public License as published  by the Free Soft- --
14-- ware  Foundation;  either version 3,  or (at your option) any later ver- --
15-- sion.  GNAT is distributed in the hope that it will be useful, but WITH- --
16-- OUT ANY WARRANTY;  without even the  implied warranty of MERCHANTABILITY --
17-- or FITNESS FOR A PARTICULAR PURPOSE.                                     --
18--                                                                          --
19-- As a special exception under Section 7 of GPL version 3, you are granted --
20-- additional permissions described in the GCC Runtime Library Exception,   --
21-- version 3.1, as published by the Free Software Foundation.               --
22--                                                                          --
23-- You should have received a copy of the GNU General Public License and    --
24-- a copy of the GCC Runtime Library Exception along with this program;     --
25-- see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see    --
26-- <http://www.gnu.org/licenses/>.                                          --
27--                                                                          --
28-- GNAT was originally developed  by the GNAT team at  New York University. --
29-- Extensive contributions were provided by Ada Core Technologies Inc.      --
30--                                                                          --
31------------------------------------------------------------------------------
32
33--  File a-numaux.adb <- a-numaux-darwin.adb
34
35package body Ada.Numerics.Aux is
36
37   -----------------------
38   -- Local subprograms --
39   -----------------------
40
41   procedure Reduce (X : in out Double; Q : out Natural);
42   --  Implements reduction of X by Pi/2. Q is the quadrant of the final
43   --  result in the range 0 .. 3. The absolute value of X is at most Pi/4.
44
45   --  The following three functions implement Chebishev approximations
46   --  of the trigonometric functions in their reduced domain.
47   --  These approximations have been computed using Maple.
48
49   function Sine_Approx (X : Double) return Double;
50   function Cosine_Approx (X : Double) return Double;
51
52   pragma Inline (Reduce);
53   pragma Inline (Sine_Approx);
54   pragma Inline (Cosine_Approx);
55
56   function Cosine_Approx (X : Double) return Double is
57      XX : constant Double := X * X;
58   begin
59      return (((((16#8.DC57FBD05F640#E-08 * XX
60              - 16#4.9F7D00BF25D80#E-06) * XX
61              + 16#1.A019F7FDEFCC2#E-04) * XX
62              - 16#5.B05B058F18B20#E-03) * XX
63              + 16#A.AAAAAAAA73FA8#E-02) * XX
64              - 16#7.FFFFFFFFFFDE4#E-01) * XX
65              - 16#3.655E64869ECCE#E-14 + 1.0;
66   end Cosine_Approx;
67
68   function Sine_Approx (X : Double) return Double is
69      XX : constant Double := X * X;
70   begin
71      return (((((16#A.EA2D4ABE41808#E-09 * XX
72              - 16#6.B974C10F9D078#E-07) * XX
73              + 16#2.E3BC673425B0E#E-05) * XX
74              - 16#D.00D00CCA7AF00#E-04) * XX
75              + 16#2.222222221B190#E-02) * XX
76              - 16#2.AAAAAAAAAAA44#E-01) * (XX * X) + X;
77   end Sine_Approx;
78
79   ------------
80   -- Reduce --
81   ------------
82
83   procedure Reduce (X : in out Double; Q : out Natural) is
84      Half_Pi     : constant := Pi / 2.0;
85      Two_Over_Pi : constant := 2.0 / Pi;
86
87      HM : constant := Integer'Min (Double'Machine_Mantissa / 2, Natural'Size);
88      M  : constant Double := 0.5 + 2.0**(1 - HM); -- Splitting constant
89      P1 : constant Double := Double'Leading_Part (Half_Pi, HM);
90      P2 : constant Double := Double'Leading_Part (Half_Pi - P1, HM);
91      P3 : constant Double := Double'Leading_Part (Half_Pi - P1 - P2, HM);
92      P4 : constant Double := Double'Leading_Part (Half_Pi - P1 - P2 - P3, HM);
93      P5 : constant Double := Double'Leading_Part (Half_Pi - P1 - P2 - P3
94                                                                 - P4, HM);
95      P6 : constant Double := Double'Model (Half_Pi - P1 - P2 - P3 - P4 - P5);
96      K  : Double;
97
98   begin
99      --  For X < 2.0**HM, all products below are computed exactly.
100      --  Due to cancellation effects all subtractions are exact as well.
101      --  As no double extended floating-point number has more than 75
102      --  zeros after the binary point, the result will be the correctly
103      --  rounded result of X - K * (Pi / 2.0).
104
105      K := X * Two_Over_Pi;
106      while abs K >= 2.0 ** HM loop
107         K := K * M - (K * M - K);
108         X :=
109           (((((X - K * P1) - K * P2) - K * P3) - K * P4) - K * P5) - K * P6;
110         K := X * Two_Over_Pi;
111      end loop;
112
113      --  If K is not a number (because X was not finite) raise exception
114
115      if K /= K then
116         raise Constraint_Error;
117      end if;
118
119      K := Double'Rounding (K);
120      Q := Integer (K) mod 4;
121      X := (((((X - K * P1) - K * P2) - K * P3)
122                  - K * P4) - K * P5) - K * P6;
123   end Reduce;
124
125   ---------
126   -- Cos --
127   ---------
128
129   function Cos (X : Double) return Double is
130      Reduced_X : Double := abs X;
131      Quadrant  : Natural range 0 .. 3;
132
133   begin
134      if Reduced_X > Pi / 4.0 then
135         Reduce (Reduced_X, Quadrant);
136
137         case Quadrant is
138            when 0 =>
139               return Cosine_Approx (Reduced_X);
140
141            when 1 =>
142               return Sine_Approx (-Reduced_X);
143
144            when 2 =>
145               return -Cosine_Approx (Reduced_X);
146
147            when 3 =>
148               return Sine_Approx (Reduced_X);
149         end case;
150      end if;
151
152      return Cosine_Approx (Reduced_X);
153   end Cos;
154
155   ---------
156   -- Sin --
157   ---------
158
159   function Sin (X : Double) return Double is
160      Reduced_X : Double := X;
161      Quadrant  : Natural range 0 .. 3;
162
163   begin
164      if abs X > Pi / 4.0 then
165         Reduce (Reduced_X, Quadrant);
166
167         case Quadrant is
168            when 0 =>
169               return Sine_Approx (Reduced_X);
170
171            when 1 =>
172               return Cosine_Approx (Reduced_X);
173
174            when 2 =>
175               return Sine_Approx (-Reduced_X);
176
177            when 3 =>
178               return -Cosine_Approx (Reduced_X);
179         end case;
180      end if;
181
182      return Sine_Approx (Reduced_X);
183   end Sin;
184
185end Ada.Numerics.Aux;
186