----------------------------------------------------------------------------
   --  Natural_Logarithm

   --  Implementation Notes:
   --    - Let f = (x - 1) / (x + 1).
   --    - Uses series
   --        ln (x) = 2 * (f + f**3 / 3 + f**5 / 5 + ...) when 1/2 < x < 2,
   --        ln (x) = 2 * ln (sqrt (x)) otherwise.
   function Natural_Logarithm (X : Decimal) return Decimal is

      A : Decimal := X;
      L : Decimal;
      P : Decimal := To_Number (2);

      Old_Model : constant Rounding_Model := Get_Rounding;

   begin  --  Natural_Logarithm

      if A <= To_Number (0) then
         raise Non_Existent_Logarithm;
      end if;

      Set_Rounding
        (Scale  => Get_Scale (Old_Model) + 6,
         Method => Truncate);

      Normalize : declare

         Two  : constant Decimal := To_Number (2);
         Half : constant Decimal := To_Number ("0.5");
         S    :          Natural := 0;

      begin

         while A >= Two or else A <= Half loop
            A := Root (A, Two);
            P := P * Two;
            S := S + 1;
         end loop;

         Set_Scale (Get_Scale + S);

      end Normalize;

      Computation : declare

         F :          Decimal := (A - To_Number (1)) / (A + To_Number (1));
         M : constant Decimal := F * F;
         N :          Decimal := To_Number (3);
         S :          Decimal;

      begin

         L := F;

         loop
            F := F * M;
            S := F / N;
            exit when S = To_Number (0);
            L := L + S;
            N := N + To_Number (2);
         end loop;

      end Computation;

      Denormalize : begin

         L := L * P;

      end Denormalize;

      Round
        (Self  => L,
         Model => New_Model
           (Scale  => Get_Scale (Old_Model),
            Method => Truncate));

      Set_Rounding (Old_Model);

      return L;

   end Natural_Logarithm;