DataMuseum.dk

top - metrics - download

    Length: 6020 (0x1784)
    Types: TextFile
    Names: »B«

└─⟦180fe333a⟧ Bits:30000405 8mm tape, Rational 1000, SW CATALOG, 10_20_0
└─⟦180fe333a⟧ Bits:30000537 8mm tape, Rational 1000, SW Catalog 10_20_0
    └─⟦5cb1d1d7f⟧ »DATA« 
        └─⟦3b1ee7bd8⟧ 
            └─⟦this⟧ 

TextFile

with System;
with Natural_Package;
package body Raised is
    --==================================================================
    -- This package encapsulates an overloading of the exponential
    -- operator that will raise a floating point number to a floating
    -- point power.
    --
    -- This package attempts to do some minimal numerical analysis in
    -- order to compute the result to the requested accuracy.  The
    -- expression for the termination of the loop comes from the book
    -- "Portability and Style in Ada" by Nissen and Wallis.
    --
    -- The method used to compute the exponentail is based on the
    -- fact that
    --   a ** x = e ** (x * ln(a))
    -- where ln is the natural log function and e is its base.
    -- This is then further refined by saying that x = y.z where
    -- y is the integer part of the number and z is the fractional
    -- part.  Therefore:
    --   a ** x = a ** y.z = a ** (y.0 + 0.z) = (a ** y.0) * (a ** 0.z)
    -- and then calculating a ** y.0 as a ** y using Ada's exponential
    -- function and expanding the other half to be:
    --   (a ** y) * [e ** (0.z * ln(a))]
    -- Now e ** (0.z * ln(a)) is calculated using the Taylor series
    -- expansion:
    --   e ** (0.z * ln(a)) = 1 + 0.z*ln(a) + (0.z*ln(a))**2/2! + ...
    --
    --
    -- Will raise the exception Negative_Base if the user tries to raise
    -- a number <= 0 to a fractional power.
    --
    -- Will raise the exception Value_Too_Large if the result is too
    -- large to compute on the particular machine
    --
    --
    -- Version 2.0 December 6, 1985
    --
    -- Written by Brad Balfour with help and suggestions from
    -- Ed Berard, Johan Margono and Gary Russell
    --==================================================================


    function "**" (Some_Base : in Base; To_The_Power : in Exponent)
                  return Return_Type is
        --===============================================================
        --
        -- calculates a ** x for a and x being floating point numbers.
        --
        -- Version 2.0 December 6, 1985
        --
        -- Written by Brad Balfour with help and suggestions from
        -- Ed Berard, Johan Margono and Gary Russell
        --===============================================================


        type Accurate_Float is digits System.Max_Digits;
        -- a very accurate floating point number used for all
        -- internal calculations

        package New_Natural is
           new Natural_Package (Argument => Base,
                                Return_Value => Accurate_Float);

        Numerator : Accurate_Float;
        -- holds the numerator for the Taylor series.  This
        -- will be the fractional part of the exponent, To_The_Power,
        -- times the log of the base.

        Current_Term : Accurate_Float := 0.0;
        -- the current term of the series

        Computed_Answer : Accurate_Float := 1.0;  -- takes care of the
        -- first term by
        -- initializing it to 1

        type Very_Large_Integer is range 0 .. System.Max_Int;

        Integer_Exponent : Very_Large_Integer := 0;
        --
        -- will hold the integer part of the Exponent
        --

        Term_Number : Integer := 1;
        --
        -- counts the number of terms
        --


        function Factorial (Of_Number : in Accurate_Float)
                           return Accurate_Float is
            --========================================================
            -- This calculates the factorial of a number by
            -- subtracting one and recursively calling itself.
            --
            -- a floating point number is used to improve the range of
            -- values of the (usual) integer argument.  This can be
            -- done because the integers are all model numbers in
            -- any floating point type.
            --
            --========================================================

        begin
            -- Factorial

            if Of_Number > 1.0 then

                return (Of_Number * Factorial (Of_Number => Of_Number - 1.0));

            else
                -- it's one, so don't recurse

                return 1.0;

            end if;

        end Factorial;


    begin
        -- "**"

        if Some_Base <= 0.0 then

            raise Negative_Base;

        else
            -- separate into integer and fractional part

            Integer_Exponent := Very_Large_Integer (To_The_Power);
            -- this will round it off

            if Exponent (Integer_Exponent) > To_The_Power then
                -- this must have rounded up, so subtract 1 to truncate
                Integer_Exponent := Integer_Exponent - 1;

            else
                -- rounded down which equals trucation

                null;

            end if;

            Numerator := Accurate_Float (To_The_Power -
                                         Exponent (Integer_Exponent));
            -- this is the fractional part.

            Numerator := Numerator * New_Natural.Log (Some_Base);

            -- the numerator is the fractional part * ln(base)


            Compute_Mantissa:
                loop

                    Current_Term := Numerator ** Term_Number /
                                       Factorial (Accurate_Float (Term_Number));

                    exit Compute_Mantissa when
                       Current_Term < Computed_Answer * Return_Type'Epsilon;

                    Computed_Answer := Computed_Answer + Current_Term;
                    Term_Number := Term_Number + 1;

                end loop Compute_Mantissa;

            return Return_Type ((Accurate_Float (Some_Base) **
                                 Integer (Integer_Exponent)) * Computed_Answer);
            -- the final answer is a**y.0 * e**(0.z*ln(a))

        end if; -- < 0 or not

    exception

        when Constraint_Error | Numeric_Error =>

            raise Value_Too_Large;

    end "**";

end Raised;