1 ------------------------------------------------------------------------------
3 -- GNAT RUN-TIME COMPONENTS --
5 -- A D A . N U M E R I C S . A U X --
8 -- (Apple OS X Version) --
10 -- Copyright (C) 1998-2016, Free Software Foundation, Inc. --
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. --
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. --
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/>. --
28 -- GNAT was originally developed by the GNAT team at New York University. --
29 -- Extensive contributions were provided by Ada Core Technologies Inc. --
31 ------------------------------------------------------------------------------
33 package body Ada
.Numerics
.Aux
is
35 -----------------------
36 -- Local subprograms --
37 -----------------------
39 function Is_Nan
(X
: Double
) return Boolean;
40 -- Return True iff X is a IEEE NaN value
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.
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.
53 function Sine_Approx
(X
: Double
) return Double
;
54 function Cosine_Approx
(X
: Double
) return Double
;
56 pragma Inline
(Reduce
);
57 pragma Inline
(Sine_Approx
);
58 pragma Inline
(Cosine_Approx
);
64 function Cosine_Approx
(X
: Double
) return Double
is
65 XX
: constant Double
:= X
* X
;
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;
80 function Sine_Approx
(X
: Double
) return Double
is
81 XX
: constant Double
:= X
* X
;
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
;
95 function Is_Nan
(X
: Double
) return Boolean is
97 -- The IEEE NaN values are the only ones that do not equal themselves
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
;
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
118 P6
: constant Double
:= Double
'Model (Half_Pi
- P1
- P2
- P3
- P4
- P5
);
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).
129 K
:= X
* Two_Over_Pi
;
130 while abs K
>= 2.0**HM
loop
131 K
:= K
* M
- (K
* M
- K
);
133 (((((X
- K
* P1
) - K
* P2
) - K
* P3
) - K
* P4
) - K
* P5
) - K
* P6
;
134 K
:= X
* Two_Over_Pi
;
137 -- If K is not a number (because X was not finite) raise exception
140 raise Constraint_Error
;
143 -- Go through an integer temporary so as to use machine instructions
145 R
:= Integer (Double
'Rounding (K
));
148 X
:= (((((X
- K
* P1
) - K
* P2
) - K
* P3
) - K
* P4
) - K
* P5
) - K
* P6
;
155 function Cos
(X
: Double
) return Double
is
156 Reduced_X
: Double
:= abs X
;
157 Quadrant
: Natural range 0 .. 3;
160 if Reduced_X
> Pi
/ 4.0 then
161 Reduce
(Reduced_X
, Quadrant
);
165 return Cosine_Approx
(Reduced_X
);
168 return Sine_Approx
(-Reduced_X
);
171 return -Cosine_Approx
(Reduced_X
);
174 return Sine_Approx
(Reduced_X
);
178 return Cosine_Approx
(Reduced_X
);
185 function Sin
(X
: Double
) return Double
is
186 Reduced_X
: Double
:= X
;
187 Quadrant
: Natural range 0 .. 3;
190 if abs X
> Pi
/ 4.0 then
191 Reduce
(Reduced_X
, Quadrant
);
195 return Sine_Approx
(Reduced_X
);
198 return Cosine_Approx
(Reduced_X
);
201 return Sine_Approx
(-Reduced_X
);
204 return -Cosine_Approx
(Reduced_X
);
208 return Sine_Approx
(Reduced_X
);
211 end Ada
.Numerics
.Aux
;