DataMuseum.dk

top - download

    Length: 5240 (0x1478)
    Types: TextFile
    Names: »B«

└─⟦5f3412b64⟧ Bits:30000745 8mm tape, Rational 1000, ENVIRONMENT 12_6_5 TOOLS 
    └─ ⟦91c658230⟧ »DATA« 
        └─⟦458657fb6⟧ 
            └─⟦1472c4407⟧ 
                └─⟦this⟧ 
└─⟦d10a02448⟧ Bits:30000409 8mm tape, Rational 1000, ENVIRONMENT, D_12_7_3
    └─ ⟦fc9b38f02⟧ »DATA« 
        └─⟦9b46a407a⟧ 
            └─⟦2e03b931c⟧ 
                └─⟦this⟧ 

TextFile

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;