File : a-numaux-darwin.adb
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-2014, 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 -- --
20 -- --
21 -- --
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 package body Ada.Numerics.Aux is
34
35 -----------------------
36 -- Local subprograms --
37 -----------------------
38
39 procedure Reduce (X : in out Double; Q : out Natural);
40 -- Implements reduction of X by Pi/2. Q is the quadrant of the final
41 -- result in the range 0 .. 3. The absolute value of X is at most Pi/4.
42
43 -- The following three functions implement Chebishev approximations
44 -- of the trigonometric functions in their reduced domain.
45 -- These approximations have been computed using Maple.
46
47 function Sine_Approx (X : Double) return Double;
48 function Cosine_Approx (X : Double) return Double;
49
50 pragma Inline (Reduce);
51 pragma Inline (Sine_Approx);
52 pragma Inline (Cosine_Approx);
53
54 function Cosine_Approx (X : Double) return Double is
55 XX : constant Double := X * X;
56 begin
57 return (((((16#8.DC57FBD05F640#E-08 * XX
58 - 16#4.9F7D00BF25D80#E-06) * XX
59 + 16#1.A019F7FDEFCC2#E-04) * XX
60 - 16#5.B05B058F18B20#E-03) * XX
61 + 16#A.AAAAAAAA73FA8#E-02) * XX
62 - 16#7.FFFFFFFFFFDE4#E-01) * XX
63 - 16#3.655E64869ECCE#E-14 + 1.0;
64 end Cosine_Approx;
65
66 function Sine_Approx (X : Double) return Double is
67 XX : constant Double := X * X;
68 begin
69 return (((((16#A.EA2D4ABE41808#E-09 * XX
70 - 16#6.B974C10F9D078#E-07) * XX
71 + 16#2.E3BC673425B0E#E-05) * XX
72 - 16#D.00D00CCA7AF00#E-04) * XX
73 + 16#2.222222221B190#E-02) * XX
74 - 16#2.AAAAAAAAAAA44#E-01) * (XX * X) + X;
75 end Sine_Approx;
76
77 ------------
78 -- Reduce --
79 ------------
80
81 procedure Reduce (X : in out Double; Q : out Natural) is
82 Half_Pi : constant := Pi / 2.0;
83 Two_Over_Pi : constant := 2.0 / Pi;
84
85 HM : constant := Integer'Min (Double'Machine_Mantissa / 2, Natural'Size);
86 M : constant Double := 0.5 + 2.0**(1 - HM); -- Splitting constant
87 P1 : constant Double := Double'Leading_Part (Half_Pi, HM);
88 P2 : constant Double := Double'Leading_Part (Half_Pi - P1, HM);
89 P3 : constant Double := Double'Leading_Part (Half_Pi - P1 - P2, HM);
90 P4 : constant Double := Double'Leading_Part (Half_Pi - P1 - P2 - P3, HM);
91 P5 : constant Double := Double'Leading_Part (Half_Pi - P1 - P2 - P3
92 - P4, HM);
93 P6 : constant Double := Double'Model (Half_Pi - P1 - P2 - P3 - P4 - P5);
94 K : Double;
95
96 begin
97 -- For X < 2.0**HM, all products below are computed exactly.
98 -- Due to cancellation effects all subtractions are exact as well.
99 -- As no double extended floating-point number has more than 75
100 -- zeros after the binary point, the result will be the correctly
101 -- rounded result of X - K * (Pi / 2.0).
102
103 K := X * Two_Over_Pi;
104 while abs K >= 2.0 ** HM loop
105 K := K * M - (K * M - K);
106 X :=
107 (((((X - K * P1) - K * P2) - K * P3) - K * P4) - K * P5) - K * P6;
108 K := X * Two_Over_Pi;
109 end loop;
110
111 -- If K is not a number (because X was not finite) raise exception
112
113 if K /= K then
114 raise Constraint_Error;
115 end if;
116
117 K := Double'Rounding (K);
118 Q := Integer (K) mod 4;
119 X := (((((X - K * P1) - K * P2) - K * P3)
120 - K * P4) - K * P5) - K * P6;
121 end Reduce;
122
123 ---------
124 -- Cos --
125 ---------
126
127 function Cos (X : Double) return Double is
128 Reduced_X : Double := abs X;
129 Quadrant : Natural range 0 .. 3;
130
131 begin
132 if Reduced_X > Pi / 4.0 then
133 Reduce (Reduced_X, Quadrant);
134
135 case Quadrant is
136 when 0 =>
137 return Cosine_Approx (Reduced_X);
138
139 when 1 =>
140 return Sine_Approx (-Reduced_X);
141
142 when 2 =>
143 return -Cosine_Approx (Reduced_X);
144
145 when 3 =>
146 return Sine_Approx (Reduced_X);
147 end case;
148 end if;
149
150 return Cosine_Approx (Reduced_X);
151 end Cos;
152
153 ---------
154 -- Sin --
155 ---------
156
157 function Sin (X : Double) return Double is
158 Reduced_X : Double := X;
159 Quadrant : Natural range 0 .. 3;
160
161 begin
162 if abs X > Pi / 4.0 then
163 Reduce (Reduced_X, Quadrant);
164
165 case Quadrant is
166 when 0 =>
167 return Sine_Approx (Reduced_X);
168
169 when 1 =>
170 return Cosine_Approx (Reduced_X);
171
172 when 2 =>
173 return Sine_Approx (-Reduced_X);
174
175 when 3 =>
176 return -Cosine_Approx (Reduced_X);
177 end case;
178 end if;
179
180 return Sine_Approx (Reduced_X);
181 end Sin;
182
183 end Ada.Numerics.Aux;