separate ( Generic_Elementary_Functions )

function KF_AtnSm( Y : Common_Float ) return Common_Float is

-- On input, Y is a floating point value in Common_Float, |Y| <= 1/16;
-- On output, a Common_Float value is returned which represents the 
-- approximation of atan(Y)-Y.

   Q                    : Common_Float;

begin


-- 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, S, Poly : Working_Float;
         begin
            R := Working_Float( Y );
            S := R * R;
            Poly := -R*S*( 0.33333_23001 -
                       S*( 0.19918_59719 ));
            Q := 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, S, Poly : Working_Float;
         begin
            R := Working_Float( Y );
            S := R * R;
            Poly := -R*S*( 0.33333_33333_33331_70500 -
                       S*( 0.19999_99999_93932_23405 -
                       S*( 0.14285_71356_18071_69030 -
                       S*( 0.11110_73628_35145_25407 -
                       S*( 0.90029_81028_54497_84439E-01 ))))); 
            Q := Common_Float( Poly );
         end;

      when 16 =>

         declare
            type Working_Float is digits 
               (16+System.Max_Digits - abs(16-System.Max_Digits))/2;
            R, S, Poly : Working_Float; 
         begin
            R :=  Working_Float( Y );
            S := R * R;
            Poly := -R*S*( 0.33333_33333_33331_70500 -
                       S*( 0.19999_99999_93932_23405 -
                       S*( 0.14285_71356_18071_69030 -
                       S*( 0.11110_73628_35145_25407 -
                       S*( 0.90029_81028_54497_84439E-01 ))))); 
            Q := Common_Float( Poly );

         end;

      when 17..18 =>

         declare
            type Working_Float is digits
               (18+System.Max_Digits - abs(18-System.Max_Digits))/2;
            R, S, Poly : Working_Float;
         begin
            R := Working_Float( Y );
            S := R * R;
            Poly := -R*S*( 0.33333_33333_33333_33130_94892 -
                       S*( 0.19999_99999_99989_99552_86683 -
                       S*( 0.14285_71428_40674_81920_91353 -
                       S*( 0.11111_10986_45979_04500_50686 -
                       S*( 0.90904_31858_09744_75086_35684E-01 - 
                       S*( 0.76023_22269_40038_42302_02382E-01 ))))));
            Q := Common_Float( Poly );
         end;

      when 19..27 =>

         declare
            type Working_Float is digits
               (27+System.Max_Digits - abs(27-System.Max_Digits))/2;
            R, S, Poly : Working_Float; 
         begin
            R := Working_Float( Y );
            S := R * R;
            Poly := -R*S*( 0.33333_33333_33333_33333_33333_30640_688 -
                       S*( 0.19999_99999_99999_99999_99725_07169_093 -
                       S*( 0.14285_71428_57142_85704_61906_79692_957 -
                       S*( 0.11111_11111_11110_94447_44685_62289_640 -
                       S*( 0.90909_09090_89284_08030_74591_43510_263E-01 -
                       S*( 0.76923_07682_73618_54039_14363_83712_288E-01 -
                       S*( 0.66666_63192_87231_75829_95036_25183_458E-01 -
                       S*( 0.58815_92881_06830_95028_64129_28595_505E-01 -
                       S*( 0.51711_22330_76482_53731_52208_70654_932E-01 
                              )))))))));
            Q := Common_Float( Poly );
         end;

      when 28..33 =>

         declare
            type Working_Float is digits
               (33+System.Max_Digits - abs(33-System.Max_Digits))/2;
            R, S, Poly : Working_Float; 
         begin
            R := Working_Float( Y );
            S := R * R;
            Poly := -R*S*( 0.33333_33333_33333_33333_33333_33333_3304 -
                       S*( 0.19999_99999_99999_99999_99999_99956_0187 -
                       S*( 0.14285_71428_57142_85714_28569_15903_2129 -
                       S*( 0.11111_11111_11111_11111_05278_66512_9097 -
                       S*( 0.90909_09090_90909_08220_63173_05946_2434E-01 -
                       S*( 0.76923_07692_30761_06801_67642_90322_8423E-01 -
                       S*( 0.66666_66666_61651_60744_77042_68877_1500E-01 -
                       S*( 0.58823_52920_68496_03403_53270_56134_3084E-01 -
                       S*( 0.52631_52374_60356_91863_62352_54421_7989E-01 -
                       S*( 0.47609_62465_86213_53088_72992_19065_8779E-01 -
                       S*( 0.42553_42964_37218_11922_99573_45059_6737E-01 
                             )))))))))));
            Q := Common_Float( Poly );
         end;

      when others =>

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

   end case;

-- This completes the core approximation.

   return( Q );


end KF_AtnSm;