old-autovect-branch/gcc/ada/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-2005, 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 2,  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.  See the GNU General Public License --
  18 -- for  more details.  You should have  received  a copy of the GNU General --
  19 -- Public License  distributed with GNAT;  see file COPYING.  If not, write --
  20 -- to  the  Free Software Foundation,  51  Franklin  Street,  Fifth  Floor, --
  21 -- Boston, MA 02110-1301, USA.                                              --
  22 --                                                                          --
  23 -- As a special exception,  if other files  instantiate  generics from this --
  24 -- unit, or you link  this unit with other files  to produce an executable, --
  25 -- this  unit  does not  by itself cause  the resulting  executable  to  be --
  26 -- covered  by the  GNU  General  Public  License.  This exception does not --
  27 -- however invalidate  any other reasons why  the executable file  might be --
  28 -- covered by the  GNU Public License.                                      --
  29 --                                                                          --
  30 -- GNAT was originally developed  by the GNAT team at  New York University. --
  31 -- Extensive contributions were provided by Ada Core Technologies Inc.      --
  32 --                                                                          --
  33 ------------------------------------------------------------------------------
  34
  35 --  File a-numaux.adb <- a-numaux-d arwin.adb
  36
  37 package body Ada.Numerics.Aux is
  38
  39    -----------------------
  40    -- Local subprograms --
  41    -----------------------
  42
  43    procedure Reduce (X : in out Double; Q : out Natural);
  44    --  Implements reduction of X by Pi/2. Q is the quadrant of the final
  45    --  result in the range 0 .. 3. The absolute value of X is at most Pi/4.
  46
  47    --  The following three functions implement Chebishev approximations
  48    --  of the trigoniometric functions in their reduced domain.
  49    --  These approximations have been computed using Maple.
  50
  51    function Sine_Approx (X : Double) return Double;
  52    function Cosine_Approx (X : Double) return Double;
  53
  54    pragma Inline (Reduce);
  55    pragma Inline (Sine_Approx);
  56    pragma Inline (Cosine_Approx);
  57
  58    function Cosine_Approx (X : Double) return Double is
  59       XX : constant Double := X * X;
  60    begin
  61       return (((((16#8.DC57FBD05F640#E-08 * XX
  62               - 16#4.9F7D00BF25D80#E-06) * XX
  63               + 16#1.A019F7FDEFCC2#E-04) * XX
  64               - 16#5.B05B058F18B20#E-03) * XX
  65               + 16#A.AAAAAAAA73FA8#E-02) * XX
  66               - 16#7.FFFFFFFFFFDE4#E-01) * XX
  67               - 16#3.655E64869ECCE#E-14 + 1.0;
  68    end Cosine_Approx;
  69
  70    function Sine_Approx (X : Double) return Double is
  71       XX : constant Double := X * X;
  72    begin
  73       return (((((16#A.EA2D4ABE41808#E-09 * XX
  74               - 16#6.B974C10F9D078#E-07) * XX
  75               + 16#2.E3BC673425B0E#E-05) * XX
  76               - 16#D.00D00CCA7AF00#E-04) * XX
  77               + 16#2.222222221B190#E-02) * XX
  78               - 16#2.AAAAAAAAAAA44#E-01) * (XX * X) + X;
  79    end Sine_Approx;
  80
  81    ------------
  82    -- Reduce --
  83    ------------
  84
  85    procedure Reduce (X : in out Double; Q : out Natural) is
  86       Half_Pi     : constant := Pi / 2.0;
  87       Two_Over_Pi : constant := 2.0 / Pi;
  88
  89       HM : constant := Integer'Min (Double'Machine_Mantissa / 2, Natural'Size);
  90       M  : constant Double := 0.5 + 2.0**(1 - HM); -- Splitting constant
  91       P1 : constant Double := Double'Leading_Part (Half_Pi, HM);
  92       P2 : constant Double := Double'Leading_Part (Half_Pi - P1, HM);
  93       P3 : constant Double := Double'Leading_Part (Half_Pi - P1 - P2, HM);
  94       P4 : constant Double := Double'Leading_Part (Half_Pi - P1 - P2 - P3, HM);
  95       P5 : constant Double := Double'Leading_Part (Half_Pi - P1 - P2 - P3
  96                                                                  - P4, HM);
  97       P6 : constant Double := Double'Model (Half_Pi - P1 - P2 - P3 - P4 - P5);
  98       K  : Double;
  99
 100    begin
 101       --  For X < 2.0**HM, all products below are computed exactly.
 102       --  Due to cancellation effects all subtractions are exact as well.
 103       --  As no double extended floating-point number has more than 75
 104       --  zeros after the binary point, the result will be the correctly
 105       --  rounded result of X - K * (Pi / 2.0).
 106
 107       K := X * Two_Over_Pi;
 108       while abs K >= 2.0 ** HM loop
 109          K := K * M - (K * M - K);
 110          X :=
 111            (((((X - K * P1) - K * P2) - K * P3) - K * P4) - K * P5) - K * P6;
 112          K := X * Two_Over_Pi;
 113       end loop;
 114
 115       --  If K is not a number (because X was not finite) raise exception
 116
 117       if K /= K then
 118          raise Constraint_Error;
 119       end if;
 120
 121       K := Double'Rounding (K);
 122       Q := Integer (K) mod 4;
 123       X := (((((X - K * P1) - K * P2) - K * P3)
 124                   - K * P4) - K * P5) - K * P6;
 125    end Reduce;
 126
 127    ---------
 128    -- Cos --
 129    ---------
 130
 131    function Cos (X : Double) return Double is
 132       Reduced_X : Double := abs X;
 133       Quadrant  : Natural range 0 .. 3;
 134
 135    begin
 136       if Reduced_X > Pi / 4.0 then
 137          Reduce (Reduced_X, Quadrant);
 138
 139          case Quadrant is
 140             when 0 =>
 141                return Cosine_Approx (Reduced_X);
 142
 143             when 1 =>
 144                return Sine_Approx (-Reduced_X);
 145
 146             when 2 =>
 147                return -Cosine_Approx (Reduced_X);
 148
 149             when 3 =>
 150                return Sine_Approx (Reduced_X);
 151          end case;
 152       end if;
 153
 154       return Cosine_Approx (Reduced_X);
 155    end Cos;
 156
 157    ---------
 158    -- Sin --
 159    ---------
 160
 161    function Sin (X : Double) return Double is
 162       Reduced_X : Double := X;
 163       Quadrant  : Natural range 0 .. 3;
 164
 165    begin
 166       if abs X > Pi / 4.0 then
 167          Reduce (Reduced_X, Quadrant);
 168
 169          case Quadrant is
 170             when 0 =>
 171                return Sine_Approx (Reduced_X);
 172
 173             when 1 =>
 174                return Cosine_Approx (Reduced_X);
 175
 176             when 2 =>
 177                return Sine_Approx (-Reduced_X);
 178
 179             when 3 =>
 180                return -Cosine_Approx (Reduced_X);
 181          end case;
 182       end if;
 183
 184       return Sine_Approx (Reduced_X);
 185    end Sin;
 186
 187 end Ada.Numerics.Aux;