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