From:	SCOVCB::GALVIN_FX    31-JAN-1989 09:00
To:	ARISIA::EVERHART
Subj:	

package body MATH is

   FLOAT_MAX_INT : constant float := float (integer'last);
   EPSLON : constant := 1.0E-4;
   PERCENT_ERROR : float;

   --    EPSLON        = relative amt. of error permissible btwn.
   --                    two successive terms in series expansion.
   --    PERCENT_ERROR = % error btwn. two successive terms
   --                    in series expansion.

   function LN (Z : in float) return float is

      --  {Natural Log of Z} <--  rx + (s*rz**2)/2 + (s**2)(rz**3)/3
      --                             +     ...     + (s**(n-1)(rz**n)/n +   ...
      --
      --       where  rz  <=>  MULTIPLIER,  n indicates the term no.,
      --              and
      --              if  (Z is greater than or equal to 2)  then
      --                 rz = (Z-1)/Z  and   s = +1
      --              else if  (Z is greater than 0, but less than 2)  then
      --                 rz = (Z-1)    and   s = -1
      --              else if  (Z is non-positive, erroneous)  then
      --                 let  {Natural Log of Z}  <--  float'first
      --              end if;


      MULTIPLIER : float := Z - 1.0;  -- Init. multiplier to (Z-1).
      LOG_Z_BASE_E : float := 0.0;    -- Init. {Natural log of Z} to 0.
      TERM_NUMBER : integer := 1;     -- Init. term # in series expansion.
      NUMERATOR, TERM : float;        -- NUMERATOR is the part of the 
                                      -- term which is the MULTIPLIER raised
                                      -- to a power determined by the its
                                      -- term number.
                                      -- TERM  <--  NUMERATOR / TERM_NUMBER.

   begin
      if  Z <= 0.0  then
         LOG_Z_BASE_E := float'first;
      elsif Z = 1.0 then
         LOG_Z_BASE_E := 0.0;
      else
         PERCENT_ERROR := 100.0;
         if  Z >= 2.0  then
            MULTIPLIER := MULTIPLIER / Z;  -- MULTIPLIER  <--  (Z-1)/Z.
            NUMERATOR := MULTIPLIER;       -- Init. numerator to (Z-1)/Z.
            LOG_Z_BASE_E := NUMERATOR;     -- Initialize series expansion.
         elsif  (0.0 < Z)  and  (Z < 2.0)  then
            NUMERATOR := MULTIPLIER;       -- Init. 1st numerator to (Z-1).
            LOG_Z_BASE_E := NUMERATOR;     -- Initialize series expansion.
            MULTIPLIER := - MULTIPLIER;    -- MULTIPLIER  <--  -(Z-1).
         end if;

         -- Until % of error is sufficiently small :
         while  ((PERCENT_ERROR > EPSLON) and then
                 (TERM_NUMBER < 13))  loop
            TERM_NUMBER := TERM_NUMBER + 1;
            NUMERATOR := NUMERATOR * MULTIPLIER;
            TERM := NUMERATOR / float(TERM_NUMBER);
            LOG_Z_BASE_E := LOG_Z_BASE_E + TERM;
            PERCENT_ERROR := TERM / LOG_Z_BASE_E;
         end loop;
      end if;
      return LOG_Z_BASE_E;
   end LN;

   function LN (Z : in integer) return float is
   begin 
      return  LN (float (Z));
   end LN;


   function EXP (Z : in float) return float is

      --  {EXP (Z) equals (e ** Z)}  <--  1  +  Z  +  (Z**2)/2!  +  (Z**3)/3!
      --                                     + ... +  (Z**n)/n!  +  ...
      --       where  n indicates the term no.,  and  Z is the input parameter.


      NUMERATOR, TERM : float := 1.0;    -- Init. next term in series.
      TERM_NUMBER : integer := 1;        -- Init. term # in series expansion.
      EXPONENTIAL_VALUE : float := TERM; -- Init. val. of exponential function.
      MAX_EXPONENTIAL : constant float := LN (float'last);

   begin
      if  Z >= MAX_EXPONENTIAL  then
         EXPONENTIAL_VALUE := float'last;
      elsif Z = 0.0 then
         EXPONENTIAL_VALUE := 1.0;
      else
         PERCENT_ERROR := 100.0;
         -- Until % of error btwn. terms is sufficiently small :
         while  ((PERCENT_ERROR > EPSLON)  and then
                 (TERM_NUMBER < 13))  loop
            NUMERATOR := NUMERATOR * Z;
            TERM := NUMERATOR / float(FACTORIAL(TERM_NUMBER));
            EXPONENTIAL_VALUE := EXPONENTIAL_VALUE + TERM;
            TERM_NUMBER := TERM_NUMBER + 1;
            PERCENT_ERROR := TERM / EXPONENTIAL_VALUE;
         end loop;
      end if;
      return EXPONENTIAL_VALUE;
   end EXP;

   function EXP (Z : in integer) return float is
   begin
      return EXP (float (Z));
   end EXP;


   function "**" (X, Y : in float) return float is

      --   {x ** y}  <=>  {e  **  [y  *  (log of x base e) ] }

      INT_Y : integer := integer (Y);

   begin
--      if  Y = 0.0  then
--         return 1.0;
      if float (INT_Y) = Y then
         return X ** INT_Y;
      elsif X = 0.0 then
         if Y > 0.0 then
            return 0.0;
         else
            raise DIVISION_BY_ZERO;
         end if;
      elsif X < 0.0 then
         raise FLOATING_ROOT_OF_NEGATIVE_BASE;
      else
         return EXP (Y * LN (X));
      end if;
   end "**";


   function MIN (X, Y : in float) return float is
   begin
      if X < Y then
         return X;
      else
         return Y;
      end if;
   end MIN;


   function MAX (X, Y : in float) return float is
   begin
      if X > Y then
         return X;
      else
         return Y;
      end if;
   end MAX;


   function MIN (X, Y : in integer) return integer is
   begin
      if X < Y then
         return X;
      else
         return Y;
      end if;
   end MIN;


   function MAX (X, Y : in integer) return integer is
   begin
      if X > Y then
         return X;
      else
         return Y;
      end if;
   end MAX;

   function SQRT (Z : in float) return float is
   begin 
      if Z > 0.0 then
         return  Z ** 0.5;
      elsif Z = 0.0 then
         return 0.0;
      else
         raise FLOATING_ROOT_OF_NEGATIVE_BASE;
      end if;
   end SQRT;

   function SQRT (Z : in integer) return float is
   begin 
      return  SQRT (float (Z));
   end SQRT;

   function FACTORIAL (Z : in integer) return integer is
      FACT : integer :=1;
   begin
      if Z = 0 then
         FACT := 0;
      else
         for i in 0..(Z-1) loop
            FACT := FACT * (Z-i);
         end loop;
      end if;
      return FACT;
   end FACTORIAL;

end MATH;
