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