⟦fa7f22591⟧

E3 Source Code



separate (Generic_Elementary_Functions)

-- On input, Y is a Common_Float number such that
--       |Y| <= 2**( Common_Float'Mantissa / 2 ).
-- On return, the sum R1 and R2 represents the reduced
-- argument (with respect to pi/2) to at least 3/2 of the precision
-- offered by Common_Float. Precisely
--
--  | R2 | <= Common_Float'Epsilon * | R1 |,
--
-- there is an integer value k such that
--
--  Y - k( pi/2 ) = (1 + delta) * ( R1 + R2 ),
--  | delta | <= (Common_Float'Epsilon) ** (3/2),
--
-- and
--
--  | R1 + R2 | <= (pi/4) * ( 1 + sqrt(Common_Float'Epsilon) ).
--
--
-- There are two major steps in the algorithm.
--
-- First, the value k is obtained by
--      k := round-to-nearest-integer( Y * Two_by_Pi ),
-- where Two_by_Pi is the working-precise value of (2/pi).
-- Because the magnitude of Y is never too large, the value k
-- obtained is very accurate. Moreover, note also that the value
-- k need not be identical to the integer closest to to
-- Y*(2/pi), provided that
--                    Y - k * (pi/2)
-- can be computed accurately and that value is small in magnitude
-- ( <= pi/4 roughly ).
--
-- Second, the difference Y - k*(pi/2) is calculated. In this
-- program, (pi/2) is represented by
--   pi/2 = Piby2_Piece(1) + Piby2_Piece(2) + ... +
--
--              Piby2_Piece( N_Pieces_of_Piby2 )
--
-- in a way such that  k * Piby2_Piece(J)  is exact for all
-- possible J's. Thus, to compute the reduced argument accurately,
-- we need only subtract correctly:
--
--     Y - k*Piby2_Piece(1) - k*Piby2_Piece(2) - ...
--             ... - k*Piby2_Piece( N_Pieces_of_Piby2 ).
--
-- Three steps are needed to perform the subtraction.
--
-- Step 1. Let X := Y. For J = 1, 2, ....,  perform
--            X := X - k * Piby2_Piece( J )
-- as long as the subtraction is exact. That is the case if
--    (k*Piby2_Piece(J))/2 <= X <= 2*(k*Piby2_Piece(J)).
--
-- Step 2. Suppose Step 1 terminates at J = N. Then,
--    X = Y - k*(Piby2_Piece(1) + ... + Piby2_Piece(N))
-- exactly. It can be shown that
--             X - k*Piby2_Piece(N+1)
-- can be represented exactly as
--       X - k*Piby2_Piece = X_Lead  +  X_Trail,
-- where
--     | X_Trail | <= Common_Float'Epislon * | X_Lead |.
--
-- Step 3. At this stage, X_Lead is the reduced argument
-- to the precision of Common_Float. We must now calculate
-- correction to X_Lead+X_Trail so that the reduced argument
-- can be (3/2) the accuracy of Common_Float. Therefore, we
-- modify X_Trail by
--    X_Trail := X_Trail - k*( Piby2_Piece(N+2) + ...
--                              + Piby2_Piece( good enough ) ).
--

procedure Kp_Pi2rd (Y      : in  Common_Float;
                    R1, R2 : out Common_Float;
                    I_Flag : out Quadrant_Flag) is

   X, X_Lead, X_Trail, X_Extra, Temp, F, K, K_Mod16, P, Two_P, Half_P :
      Common_Float;
   J, L1, L2 : Integer;
   Sign_Y : Quadrant_Flag;
   Same_Sign : Boolean;
   Tmp_R1 : Common_Float;
   Two_By_Pi : constant := 16#A.2F983_6E4E4_41529_FC275_7D1F5_34DDC_0DB63#E-1;

begin

-- Initialization

   X := Y;

   Sign_Y := 0;
   if X < 0.0 then
      X      := -X;
      Sign_Y := 1;
   end if;

-- Calculation of integer quotient K

   K := Round (X * Two_By_Pi);

   K_Mod16 := K * 16#0.1#;
   K_Mod16 := Truncate (K_Mod16);
   K_Mod16 := K - 16.0 * K_Mod16;
   I_Flag  := Quadrant_Flag (Integer (K_Mod16) mod 4);

   if K = 0.0 then
      R1     := Y;
      R2     := 0.0;
      I_Flag := 0;
      return;
   end if;


-- Step 1. Calcuate X := X - k*Piby2_Piece(J) whenever exact
   -- Initialization of Step 1
   J      := 1;
   P      := K * Piby2_Piece (J);
   Two_P  := P + P;
   Half_P := P * 2#1.0#E-1;

   -- The loop
   while X >= Half_P and X <= Two_P and J <= N_Pieces_Of_Piby2 - 2 loop
      if P > X then
         X := (X - Half_P) - Half_P;
      else
         X := X - P;
      end if;
      J      := J + 1;
      P      := K * Piby2_Piece (J);
      Two_P  := P + P;
      Half_P := P * 2#1.0#E-1;
   end loop;
-- End of Step 1

-- Step 2. Represent X - k*Piby2_Piece( J ) in two pieces
   P      := -K * Piby2_Piece (J);
   X_Lead := X + P;
   if X < 0.0 then
      Same_Sign := True;
   else
      Same_Sign := False;
   end if;
   if abs (X) < abs (P) then
      Temp := X;
      X    := P;
      P    := Temp;
   end if;

-- Ideally, if the arithmetic of the underlying machine is good,
-- X_Trail := (X - X_Lead) + P suffices for our purpose because
-- X - X_Lead must be exact. Unforturnately, there are machines
-- without a guard bit for subtraction. The trick below ensures
-- that  (X-F), (X_Lead-F), and (X-F) - (X_Lead-F) are all exact.
-- For a detailed explanation of why this trick work, refer to
-- pages 10-3 to 10-9 of the notes IMPLEMENTATION OF ALGORITHMS
-- by W. Kahan. This set of notes can be obtained through NTIS
-- in Virginia. Order number AD769 124 / LP, 1986.
   F := 0.0;
   if Same_Sign then
      F := (0.46 * X_Lead - X_Lead) + X_Lead;
   end if;
   X_Trail := ((X - F) - (X_Lead - F)) + P;

-- Step 3. Calcuate to higher precision.
   J       := J + 1;
   L1      := 2 * 4 * Common_Float'Base'Digits - N_Bits_Per_Piece;
   L2      := L1 / N_Bits_Per_Piece + 1;
   X_Extra := 0.0;
   if (J + L2 > N_Pieces_Of_Piby2) then
      L2 := N_Pieces_Of_Piby2 - J;
   end if;
   for I in reverse J .. J + L2 loop
      X_Extra := X_Extra + Piby2_Piece (I);
   end loop;
   X_Trail := X_Trail - K * X_Extra;


-- Finally, recast the answer in Common_Float.

   Tmp_R1 := X_Lead;
   R2     := ((X_Lead - F) - (Tmp_R1 - F)) + X_Trail;
   R1     := Tmp_R1;

   return;

end Kp_Pi2rd;

------------------------------------------------------------------------
DataMuseum.dk

Rational R1000/400

⟦fa7f22591⟧ Ada Source

Derivation

E3 Source Code

E3 Meta Data
