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-2009, 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 -- File a-numaux.adb <- a-numaux-darwin.adb
35 package body Ada.Numerics.Aux is
37 -----------------------
38 -- Local subprograms --
39 -----------------------
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.
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.
49 function Sine_Approx (X : Double) return Double;
50 function Cosine_Approx (X : Double) return Double;
52 pragma Inline (Reduce);
53 pragma Inline (Sine_Approx);
54 pragma Inline (Cosine_Approx);
56 function Cosine_Approx (X : Double) return Double is
57 XX : constant Double := X * X;
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;
68 function Sine_Approx (X : Double) return Double is
69 XX : constant Double := X * X;
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;
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;
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
95 P6 : constant Double := Double'Model (Half_Pi - P1 - P2 - P3 - P4 - P5);
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).
105 K := X * Two_Over_Pi;
106 while abs K >= 2.0 ** HM loop
107 K := K * M - (K * M - K);
109 (((((X - K * P1) - K * P2) - K * P3) - K * P4) - K * P5) - K * P6;
110 K := X * Two_Over_Pi;
113 -- If K is not a number (because X was not finite) raise exception
116 raise Constraint_Error;
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;
129 function Cos (X : Double) return Double is
130 Reduced_X : Double := abs X;
131 Quadrant : Natural range 0 .. 3;
134 if Reduced_X > Pi / 4.0 then
135 Reduce (Reduced_X, Quadrant);
139 return Cosine_Approx (Reduced_X);
142 return Sine_Approx (-Reduced_X);
145 return -Cosine_Approx (Reduced_X);
148 return Sine_Approx (Reduced_X);
152 return Cosine_Approx (Reduced_X);
159 function Sin (X : Double) return Double is
160 Reduced_X : Double := X;
161 Quadrant : Natural range 0 .. 3;
164 if abs X > Pi / 4.0 then
165 Reduce (Reduced_X, Quadrant);
169 return Sine_Approx (Reduced_X);
172 return Cosine_Approx (Reduced_X);
175 return Sine_Approx (-Reduced_X);
178 return -Cosine_Approx (Reduced_X);
182 return Sine_Approx (Reduced_X);
185 end Ada.Numerics.Aux;