⟦53a092ec0⟧

E3 Source Code



separate (Generic_Elementary_Functions)

function Kf_Power (X, Y : Float_Type) return Float_Type is

-- On input, X and Y are floating-point value in Float_Type;
-- On output, the value of X**Y is returned.

   Result : Float_Type;

   Xx, Yy, R1, R2, Z1, Z2, W1, W2, Q, Y1, Y2, Tmp, F : Common_Float;
   N, M, J : Common_Int;
   L, Ll : Positive;

   Two_To : constant array (Common_Int range -3 .. 3) of Common_Float :=
      (0.125, 0.25, 0.5, 1.0, 2.0, 4.0, 8.0);
   Thirty_Two_By_Log2 : constant :=
      16#2E.2A8EC_A5705_FC2EE_FA1FF_B41A4_74FA2_3AD5E#;

   Large_Threshold : constant Common_Float :=
      4.0 * Common_Float (Float_Type'Machine_Emax) /
         Common_Float (Float_Type'Base'Epsilon);


begin

-- Filter out exceptional cases.

   Xx := Common_Float (X);
   Yy := Common_Float (Y);

   if (Xx = 0.0) then
      if (Y < 0.0) then
         raise Constraint_Error;
      elsif (Y > 0.0) then
         return (0.0);
      else
         raise Argument_Error;
      end if;
   end if;

   if (Xx < 0.0) then
      raise Argument_Error;
   end if;

   if (Xx = 1.0 or Y = 0.0) then
      return (1.0);
   end if;

   if (Yy = 1.0) then
      return (X);
   end if;

   if abs (Yy) >= Large_Threshold then
      if (Yy > 0.0) then
         -- Yy := Common_Float (Common_Float'Machine_Radix **
         --                     Common_Float'Machine_Emax);
         -- return (Float_Type (Yy * Yy * Yy));
         raise Constraint_Error; --pbk
      else
         return (Float_Type (0.0));
      end if;
   end if;

   L  := (2 * Common_Float'Machine_Mantissa) / 5;
   Ll := Common_Float'Machine_Mantissa;

-- Get log(X) to high accuracy

   Kp_Logpw (Xx, Z1, Z2);

   Y1 := Leading_Part (Yy, L);
   Y2 := Yy - Y1;

   W1 := Y1 * Z1;
   W2 := Y1 * Z2 + (Y2 * Z1 + Y2 * Z2);

   Z1 := Leading_Part (W1 + W2, Ll - 2);       -- force storage
   Z2 := W1 - Z1;
   Z2 := W2 + Z2;
-- Z1 + Z2 is Y*log(X) to high precision


-- Now get exp(Z1+Z2)

-- Perform argument reduction in Common_Float.

   N   := Common_Int (Z1 * Thirty_Two_By_Log2);
   Tmp := Common_Float (N) * Log2_By_32_Lead;
   if abs (Z1) >= abs (Tmp) then
      R1 := Z1 - Tmp;
   else
      Tmp := 0.5 * Tmp;
      R1  := (Z1 - Tmp) - Tmp;
   end if;
   R2 := Z2 - Common_Float (N) * Log2_By_32_Trail;
   J  := N mod 32;
   M  := (N - J) / 32;

-- Step 2. Approximation.

   Q := Kf_Em1sm (R1, R2);

-- Step 3. Function value reconstruction.

-- We now reconstruct the exponential of the input argument
-- so that Exp(Z1+Z2) = 2**M * (W1 + W2).
-- The order of the computation below must be strictly observed.

   F := Two_To_Jby32 (J, Lead) + Two_To_Jby32 (J, Trail);

   W1 := Two_To_Jby32 (J, Lead);
   W2 := Two_To_Jby32 (J, Trail) + F * Q;


   case Radix is

      when 2 =>
         Tmp := W1 + W2;

      when others =>
         J   := M rem 4;
         M   := (M - J) / 4;
         W1  := W1 * Two_To (J);
         W2  := W2 * Two_To (J);
         Tmp := W1 + W2;

   end case;

   Result := Float_Type (Scale (Tmp, M));

   return (Result);


end Kf_Power;
DataMuseum.dk

Rational R1000/400

⟦53a092ec0⟧ Ada Source

Derivation

E3 Source Code

E3 Meta Data
