separate ( Generic_Elementary_Functions )

function KF_Cos( 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 * (pi/2) + remainder, where |remainder| <= pi/4.
-- On output, a Common_Float value is returned which represents the 
--   approximation of cos( Y1+Y2 ).

   R1, R2, Rsq, Q                  : Common_Float;

begin

   R1  := Y1;
   R2  := Y2;
   Rsq := R1 + R2;
   Rsq := Rsq * Rsq;

-- The following is the core approximation. We approximate
-- cos(Y1+Y2) by an even 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( R1 + R2 );
            S := R * R;
            Poly := 0.5 -
		    S*( 0.41666_65673E-01 - 
		    S*( 0.13887_77331E-02 -
		    S*( 0.24478_79547E-04 )));
            Q := 1.0 - Rsq*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( R1 + R2 );
            S := R * R;
            Poly := 0.5 -
	       S*( 0.41666_66666_66666_01903E-01 -
	       S*( 0.13888_88888_88744_48186E-02 -
	       S*( 0.24801_58728_96438_74032E-04 -
	       S*( 0.27557_31439_04201_86748E-06 -
	       S*( 0.20875_72754_71895_52619E-08 -
	       S*( 0.11359_83163_10688_30184E-10 ))))));
            Q := 1.0 - Rsq*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( R1 + R2 );
            S := R * R;
            Poly := 0.5 -
	       S*( 0.41666_66666_66666_01903E-01 -
	       S*( 0.13888_88888_88744_48186E-02 -
	       S*( 0.24801_58728_96438_74032E-04 -
	       S*( 0.27557_31439_04201_86748E-06 -
	       S*( 0.20875_72754_71895_52619E-08 -
	       S*( 0.11359_83163_10688_30184E-10 ))))));
            Q := 1.0 - Rsq*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( R1 + R2 );
            S := R * R;
            Poly := 0.5 -
		    S*( 0.41666_66666_66666_66603_42153E-01 - 
		    S*( 0.13888_88888_88888_71536_63177E-02 -
		    S*( 0.24801_58730_15691_85047_02943E-04 -
		    S*( 0.27557_31921_43827_45008_49677E-06 -
		    S*( 0.20876_75414_25307_03561_14757E-08 -
		    S*( 0.11470_26790_89435_73177_12619E-10 -
		    S*( 0.47369_66690_47318_75384_42068E-13 ))))))); 
            Q := 1.0 - Rsq*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( R1 + R2 );
            S := R * R;
            Poly := 0.5 -
		 S*( 0.41666_66666_66666_66666_66665_16383_4051E-01 -
		 S*( 0.13888_88888_88888_88888_88245_47342_9026E-02 -
		 S*( 0.24801_58730_15873_01576_57181_53770_2548E-04 -
		 S*( 0.27557_31922_39858_81219_61255_06316_8107E-06 -
		 S*( 0.20876_75698_78631_80489_36726_20113_1863E-08 -
		 S*( 0.11470_74559_61289_73694_98585_54814_7595E-10 -
		 S*( 0.47794_77003_63678_43164_38976_00318_0721E-13 -
		 S*( 0.15618_79265_92016_38912_62591_16219_0901E-15 -
		 S*( 0.40810_43877_26966_16816_80251_59059_3738E-18 ))))))))); 
            Q := 1.0 - Rsq*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( R1 + R2 );
            S := R * R;
            Poly := 0.5 -
		 S*( 0.41666_66666_66666_66666_66666_66666_6480E-01 -
		 S*( 0.13888_88888_88888_88888_88888_88877_6335E-02 -
		 S*( 0.24801_58730_15873_01587_30158_46180_0299E-04 -
		 S*( 0.27557_31922_39858_90652_55387_96436_4088E-06 -
		 S*( 0.20876_75698_78680_98976_51897_22915_8020E-08 -
		 S*( 0.11470_74559_77297_23343_14028_68248_0321E-10 -
		 S*( 0.47794_77332_38691_60941_21737_53668_2651E-13 -
		 S*( 0.15619_20696_74971_44814_63427_34672_9062E-15 -
		 S*( 0.41103_17454_13508_37011_88930_50799_6253E-18 -
		 S*( 0.88966_22926_27803_79129_98598_30659_9363E-21 -
		 S*( 0.16020_11193_30748_03892_30511_97354_2056E-23
                         ))))))))))); 
            Q := 1.0 - Rsq*Common_Float( Poly );

         end;

      when others =>

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

   end case;

-- This completes the core approximation.

   return( Q );


end KF_Cos;