separate ( Generic_Elementary_Functions )

function KF_Em1Sm( Y1, Y2 : Common_Float ) return Common_Float is

-- On input, Y1 and Y2 are floating point values in Common_Float;
--   These two variables represent the remainder of the reduced argument
--   X = N * log2/32 + remainder, where |remainder| <= log2/64.
-- On output, a Common_Float value is returned which represents the 
--   approximation of exp(Y1+Y2)-1 for small numbers.

   R1, R2, Q                    : Common_Float;

begin

   R1 := Y1;
   R2 := Y2;

-- The following is the core approximation. We approximate
-- exp(R1+R2)-1 by a polynomial. The case analysis finds both
-- a suitable floating-point type (less expensive to use than 
-- Common_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( R1 + R2 );
            Poly := R*R*( 5.00004_0481E-01 + R * 1.66667_6443E-01 );
            Q := R1 + ( R2 + 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( R1 + R2 );
            Poly := R*R*( 5.00000_00000_00000_08883E-01 + 
                      R*( 1.66666_66666_52608_78863E-01 + 
                      R*( 4.16666_66666_22607_95726E-02 + 
                      R*( 8.33336_79843_42196_16221E-03 + 
                      R*( 1.38889_49086_37771_99667E-03 )))));
            Q := R1 + ( R2 + 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( R1 + R2 );
            Poly := R*R*( 5.00000_00000_00000_08883E-01 + 
                      R*( 1.66666_66666_52608_78863E-01 + 
                      R*( 4.16666_66666_22607_95726E-02 + 
                      R*( 8.33336_79843_42196_16221E-03 + 
                      R*( 1.38889_49086_37771_99667E-03 )))));
            Q := R1 + ( R2 + 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( R1 + R2 );
            Poly := R*R*( 5.00000_00000_00000_07339E-01 + 
                      R*( 1.66666_66666_66666_69177E-01 + 
                      R*( 4.16666_66666_28680_32559E-02 + 
                      R*( 8.33333_33332_52083_91118E-03 + 
                      R*( 1.38889_44766_51246_30293E-03 + 
                      R*( 1.98413_53190_32208_33704E-04 ))))));
            Q := R1 + ( R2 + 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( R1 + R2 );
            Poly := R*R*( 4.99999_99999_99999_99999_99636_21075E-01 +
                      R*( 1.66666_66666_66666_66666_66512_04136E-01 +
                      R*( 4.16666_66666_66666_69681_59325_03184E-02 +
                      R*( 8.33333_33333_33333_40906_33326_46233E-03 +
                      R*( 1.38888_88888_81124_92492_26093_01620E-03 +
                      R*( 1.98412_69841_13983_54303_59568_15543E-04 +
                      R*( 2.48016_66086_20855_39725_92760_56125E-05 +
                      R*( 2.75574_13983_51388_82843_29291_74995E-06
                          ))))))));
            Q := R1 + ( R2 + 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( R1 + R2 );
            Poly := R*R*( 5.0E-01 +
          R*( 1.66666_66666_66666_66666_66666_66668_18891E-01 +
          R*( 4.16666_66666_66666_66666_66666_66671_98062E-02 +
          R*( 8.33333_33333_33333_33333_33182_72433_96473E-03 +
          R*( 1.38888_88888_88888_88888_88860_77788_96115E-03 +
          R*( 1.98412_69841_26984_13216_98830_39302_820E-04 +
          R*( 2.48015_87301_58730_16549_32617_44006_810E-05 +
          R*( 2.75573_19223_90497_50521_23337_44713_411E-06 +
          R*( 2.75573_19223_90383_09381_24531_22474_208E-07 +
          R*( 2.50521_67036_89710_14700_24557_88635_351E-08 +
          R*( 2.08768_06002_87469_73970_46716_40247_597E-09
                     )))))))))));
            Q := R1 + ( R2 + Common_Float( Poly ) );
         end;

      when others =>

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

   end case;

-- This completes the core approximation.

   return( Q );


end KF_Em1Sm;