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-2019, 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 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 211end Ada.Numerics.Aux; 212