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-- 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 33package 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 183end Ada.Numerics.Aux; 184