DataMuseum.dk

top - metrics - download

    Length: 5022 (0x139e)
    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_Atanh (Y : Common_Float) return Common_Float is

-- On input, |Y| <= 2(exp(1/16)-1) / (exp(1/16)+1).
-- On output, the value of [log(1 + Y/2) - log(1 - Y/2)]/Y  -  1
-- is returned.

-- The core approximation calculates
--     Poly = [log(1 + Y/2) - log(1 - Y/2)]/Y  -  1
-- in the type Working_Float, which is a type chosen to
-- have accuracy comparable to the base type of Float_Type.

   Result : Common_Float;

begin


-- Approximation.


-- The following is the core approximation. We approximate
--         [log(1 + Y/2) - log(1 - Y/2)]/Y  -  1
-- by a polynomial Poly. The case analysis finds both a suitable
-- floating-point type (less expensive to use than LONGEST_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 (Y * Y);
            Poly   := R * 8.33340_08285_51364E-02;
            Result := 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 (Y * Y);
            Poly := R * (8.33333_33333_33335_93622E-02 +
                         R * (1.24999_99997_81386_68903E-02 +
                              R * (2.23219_81075_85598_51206E-03)));

            Result := 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 (Y * Y);
            Poly := R * (8.33333_33333_33335_93622E-02 +
                         R * (1.24999_99997_81386_68903E-02 +
                              R * (2.23219_81075_85598_51206E-03)));

            Result := 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 (Y * Y);
            Poly := R * (8.33333_33333_33335_93622E-02 +
                         R * (1.24999_99997_81386_68903E-02 +
                              R * (2.23219_81075_85598_51206E-03)));

            Result := 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 (Y * Y);
            Poly   :=
               R *
                  (8.33333_33333_33333_33333_33334_07301_529E-02 +
                   R *
                      (1.24999_99999_99999_99998_61732_74718_869E-02 +
                       R *
                          (2.23214_28571_42866_13712_34336_23012_985E-03 +
                           R *
                              (4.34027_77751_26439_67391_35491_00214_979E-04 +
                               R *
                                  (8.87820_39767_24501_02052_39367_49695_054E-05)))));
            Result := 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 (Y * Y);
            Poly   :=
               R *
                  (8.33333_33333_33333_33333_33333_33332_96298_39318E-02 +
                   R *
                      (1.25000_00000_00000_00000_00000_93488_19499_40702E-02 +
                       R *
                          (2.23214_28571_42857_14277_26598_59261_40273_30694E-03 +
                           R *
                              (4.34027_77777_77814_30973_20354_95180_362E-04 +
                               R *
                                  (8.87784_09009_03777_78533_78449_15942_610E-05 +
                                   R *
                                      (1.87809_65740_24066_11924_19609_24471_232E-05))))));
            Result := Common_Float (Poly);
         end;

      when others =>

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

   end case;

-- This completes the core approximation.

   return (Result);

end Kf_Atanh;