separate ( Generic_Elementary_Functions )

function KF_R1pSm( Y : Common_Float ) return Common_Float is

-- On input, Y lies in the interval [-1/64, 1/64]
-- On output, the value sqrt( 1 + Y ) - 1 is returned.
-- The approximation employs both minmax polynomial and Newton Iteration.
--     
-- First, sqrt(1+Y)-1 is approximated by a polynomial of the form
--    0.5 Y - Y^2*(A1 - Y*(A2 - Y*(A3 ... ))).
-- Depending of the number of digits required, this polynomial may or may
-- not be the value return. Should the polynomial fall short of the
-- accuracy requirement, one step of the Newton iteration is applied to
-- the equation
--      (P + 1)^2   =   1 + Y, where P is the unknown.
-- Thus,
--      new P  :=   P   -  (P^2 + 2P - Y)/(2P + 2)
-- where the initial guess P is the value of the polynomial.

   Poly : Common_Float;

begin


   case Float_Type'Base'Digits is

      when 1..6 =>

         declare
            type Working_Float is digits 6;
            R, P : Working_Float;
         begin
            R := Working_Float( Y );
            P := R*0.5 -
		 R*R*( 0.12500_79870 -
		   R*( 0.62505_22999E-01 ));
            Poly := Common_Float( P );
         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, P, Two_P, Correction : Working_Float;
         begin
            R := Working_Float( Y );
            P := R*0.5 -
		 R*R*( 0.12500_79870 -
		   R*( 0.62505_22999E-01 ));
            Two_P := P + P;
            Correction := (P*P + (Two_P-R)) / (Two_P + 2.0);
            Poly := Common_Float(P) - Common_Float(Correction);
         end;

      when 16 =>

         declare
            type Working_Float is digits 
               (16+System.Max_Digits - abs(16-System.Max_Digits))/2;
            R, P, Two_P, Correction : Working_Float; 
         begin
            R := Working_Float( Y );
            P := R*0.5 -
		 R*R*( 0.12500_79870 -
		   R*( 0.62505_22999E-01 ));
            Two_P := P + P;
            Correction := (P*P + (Two_P-R)) / (Two_P + 2.0);
            Poly := Common_Float(P) - Common_Float(Correction);
         end;

      when 17..18 =>

         declare
            type Working_Float is digits
               (18+System.Max_Digits - abs(18-System.Max_Digits))/2;
            R, P, Two_P, Correction : Working_Float; 
         begin
            R := Working_Float( Y );
            P := R*0.5 -
		 R*R*( 0.125 -
		   R*( 0.62505_85095_32379_40606E-01 -
                   R*( 0.39067_55576_62793_58660E-01 )));
            Two_P := P + P;
            Correction := (P*P + (Two_P-R)) / (Two_P + 2.0);
            Poly := Common_Float(P) - Common_Float(Correction);
         end;

      when 19..27 =>

         declare
            type Working_Float is digits
               (27+System.Max_Digits - abs(27-System.Max_Digits))/2;
            R, P, Two_P, Correction : Working_Float; 
         begin
            R := Working_Float( Y );
            P := R*0.5 -
		 R*R*( 0.12500_00000_00043_29869_79603     -
		   R*( 0.62499_99949_58313_28657_64839E-01 -
                   R*( 0.39062_49897_81141_25807_82311E-01 -
                   R*( 0.27349_66395_87080_51587_13279E-01 -
                   R*( 0.20514_47322_06561_78273_82000E-01 )))));
            Two_P := P + P;
            Correction := (P*P + (Two_P-R)) / (Two_P + 2.0);
            Poly := Common_Float(P) - Common_Float(Correction);
         end;

      when 28..33 =>

         declare
            type Working_Float is digits
               (33+System.Max_Digits - abs(33-System.Max_Digits))/2;
            R, P, Two_P, Correction : Working_Float; 
         begin
            R := Working_Float( Y );
            P := R*0.5 -
		 R*R*( 0.12500_00000_00035_79677_13859_89202_617 -
		   R*( 0.62500_00000_00664_48622_98045_60336E-01 -
                   R*( 0.39062_49911_61750_24807_32583_41508E-01 -
                   R*( 0.27343_74885_03551_93697_91233_69204E-01 -
                   R*( 0.20514_00699_96599_52762_98109_70261E-01 -
                   R*( 0.16119_54134_09746_99174_73737_89068E-01 ))))));
            Two_P := P + P;
            Correction := (P*P + (Two_P-R)) / (Two_P + 2.0);
            Poly := Common_Float(P) - Common_Float(Correction);
         end;

      when others =>

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

   end case;

   return( Poly );

end KF_R1pSm;