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
	   raise Constraint_Error;
        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;