gcc/ada/libgnat/a-numaux__darwin.adb

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