separate ( Generic_Elementary_Functions )

function KF_Atanh( Y : Common_Float ) return Common_Float is

-- On input, |Y| <= 2(exp(1/16)-1) / (exp(1/16)+1).
-- On output, the value of [log(1 + Y/2) - log(1 - Y/2)]/Y  -  1 
-- is returned.

-- The core approximation calculates 
--     Poly = [log(1 + Y/2) - log(1 - Y/2)]/Y  -  1 
-- in the type Working_Float, which is a type chosen to
-- have accuracy comparable to the base type of Float_Type.

   Result : Common_Float;

begin

   
-- Approximation.


-- The following is the core approximation. We approximate
--         [log(1 + Y/2) - log(1 - Y/2)]/Y  -  1 
-- by a polynomial Poly. The case analysis finds both a suitable 
-- floating-point type (less expensive to use than LONGEST_FLOAT)
-- and an appropriate polynomial approximation that will deliver 
-- a result accurate enough with respect to Float_Type'Base'Digits.
-- Note that the upper bounds of the cases below (6, 15, 16, 18, 
-- 27, and 33) are attributes of predefined floating types of
-- common systems.

   case Float_Type'Base'Digits is

      when 1..6 =>

         declare
            type Working_Float is digits 6;
            R, Poly : Working_Float;
         begin
            R := Working_Float( Y * Y );
            Poly :=  R * 8.33340_08285_51364E-02;
            Result := Common_Float( Poly );
         end;

      when 7..15 =>

         declare
            type Working_Float is digits 
               (15+System.Max_Digits - abs(15-System.Max_Digits))/2;
            -- this is min( 15, System.Max_Digits ) 
            R, Poly : Working_Float;
         begin
            R := Working_Float( Y * Y );
            Poly := R*( 8.33333_33333_33335_93622E-02 +
                    R*( 1.24999_99997_81386_68903E-02 +
                    R*( 2.23219_81075_85598_51206E-03 )));

            Result := Common_Float( Poly );
         end;

      when 16 =>

         declare
            type Working_Float is digits 
               (16+System.Max_Digits - abs(16-System.Max_Digits))/2;
            R, Poly : Working_Float; 
         begin
            R :=  Working_Float( Y * Y );
            Poly := R*( 8.33333_33333_33335_93622E-02 +
                    R*( 1.24999_99997_81386_68903E-02 +
                    R*( 2.23219_81075_85598_51206E-03 )));

            Result := Common_Float( Poly );
         end;

      when 17..18 =>

         declare
            type Working_Float is digits
               (18+System.Max_Digits - abs(18-System.Max_Digits))/2;
            R, Poly : Working_Float;
         begin
            R := Working_Float( Y * Y );
            Poly := R*( 8.33333_33333_33335_93622E-02 +
                    R*( 1.24999_99997_81386_68903E-02 +
                    R*( 2.23219_81075_85598_51206E-03 )));

            Result := Common_Float( Poly );
         end;

      when 19..27 =>

         declare
            type Working_Float is digits
               (27+System.Max_Digits - abs(27-System.Max_Digits))/2;
            R, Poly : Working_Float; 
         begin
            R := Working_Float( Y * Y );
            Poly := R*( 8.33333_33333_33333_33333_33334_07301_529E-02 +
                    R*( 1.24999_99999_99999_99998_61732_74718_869E-02 +
                    R*( 2.23214_28571_42866_13712_34336_23012_985E-03 +
                    R*( 4.34027_77751_26439_67391_35491_00214_979E-04 +
                    R*( 8.87820_39767_24501_02052_39367_49695_054E-05
                         )))));
            Result := Common_Float( Poly );
         end;

      when 28..33 =>

         declare
            type Working_Float is digits
               (33+System.Max_Digits - abs(33-System.Max_Digits))/2;
            R, Poly : Working_Float; 
         begin
            R := Working_Float( Y * Y );
            Poly := 
             R*( 8.33333_33333_33333_33333_33333_33332_96298_39318E-02 +
             R*( 1.25000_00000_00000_00000_00000_93488_19499_40702E-02 +
             R*( 2.23214_28571_42857_14277_26598_59261_40273_30694E-03 +
             R*( 4.34027_77777_77814_30973_20354_95180_362E-04 +
             R*( 8.87784_09009_03777_78533_78449_15942_610E-05 +
             R*( 1.87809_65740_24066_11924_19609_24471_232E-05
                 ))))));
            Result := Common_Float( Poly );
         end;

      when others =>

         raise PROGRAM_ERROR;  -- assumption (1) is violated.

   end case;

-- This completes the core approximation.

   return( Result );

end KF_Atanh;