mathpack.vhd

fast fourior transform

		variable newval : real ;

    begin 
		
		-- compute root for special cases
		if X = 0.0 then
			return 0.0;
		elsif ( X = 1.0 ) then
			return 1.0;
		else
			if X = -1.0 then
				return -1.0;
			end if;
		end if;

		-- compute root for general cases
		if negative then
			xlocal := -X;
		end if;
		
		oldval := inival;
		newval := (xlocal/(oldval*oldval) + 2.0*oldval)/3.0;

		while ( abs(newval -oldval) > relative_err ) loop
			oldval := newval;
			newval :=(xlocal/(oldval*oldval) + 2.0*oldval)/3.0;
		end loop;

		if negative then
			newval := -newval;
		end if;

		return newval;

    end CBRT;

    function "**" (X : integer; Y : real) return real is
    	-- returns Y power of X ==>  X**Y;
    	-- error if X = 0 and Y <= 0.0
    	-- error if X < 0 and Y does not have an integer value
    begin
	-- check validity of argument
	if ( X = 0  ) and ( Y <= 0.0 ) then
		assert false report "X = 0 and Y <= 0.0 in X**Y" 
			severity ERROR;
		return (0.0);
	end if;

	if ( X < 0  ) and ( Y /= REAL(INTEGER(Y)) ) then
		assert false report "X < 0 and Y \= integer in X**Y" 
			severity ERROR;
		return (0.0);
	end if;

	-- compute the result
	return EXP (Y * LOG (REAL(X)));
    end "**";

    function "**" (X : real; Y : real) return real is
    	-- returns Y power of X ==>  X**Y;
    	-- error if X = 0.0 and Y <= 0.0
    	-- error if X < 0.0 and Y does not have an integer value
    begin
	-- check validity of argument
	if ( X = 0.0  ) and ( Y <= 0.0 ) then
		assert false report "X = 0.0 and Y <= 0.0 in X**Y" 
			severity ERROR;
		return (0.0);
	end if;

	if ( X < 0.0  ) and ( Y /= REAL(INTEGER(Y)) ) then
		assert false report "X < 0.0 and Y \= integer in X**Y" 
			severity ERROR;
		return (0.0);
	end if;

	-- compute the result
	return EXP (Y * LOG (X));
    end "**";

    function EXP  (X : real ) return real is
    	-- returns e**X; where e = MATH_E
  	--
  	-- This function computes the exponential using the following series:
  	--    exp(x) = 1 + x + x**2/2! + x**3/3! + ... ; x > 0
  	--
	constant eps : real := 0.000001;	-- precision criteria

    	variable reciprocal: boolean := x < 0.0;-- check sign of argument
    	variable xlocal : real := abs(x);       -- use positive value
    	variable oldval: real ;			-- following variables are
    	variable num: real ;			-- used for series evaluation
    	variable count: integer ;
    	variable denom: real ;
    	variable newval: real ;

     begin
    	-- compute value for special cases
	if X = 0.0 then
		return 1.0;
	else
		if  X = 1.0  then
			return MATH_E;
		end if;
	end if;

    	-- compute value for general cases
    	oldval := 1.0;
    	num := xlocal;
    	count := 1;
    	denom := 1.0;
    	newval:= oldval + num/denom;

    	while ( abs(newval - oldval) > eps ) loop
		oldval := newval;
		num := num*xlocal;
        	count := count +1;
		denom := denom*(real(count));
        	newval := oldval + num/denom;
    	end loop;

    	if reciprocal then
        	newval := 1.0/newval;
    	end if;

    	return newval;
     end EXP;


    function LOG (X : real ) return real is
    	-- returns natural logarithm of X; X > 0
	--
  	-- This function computes the exponential using the following series:
  	--    log(x) = 2[ (x-1)/(x+1) + (((x-1)/(x+1))**3)/3.0 + ...] ; x > 0
	--
    	
	constant eps : real := 0.000001;	-- precision criteria

    	variable xlocal: real ;			-- following variables are
    	variable oldval: real ;			-- used to evaluate the series
    	variable xlocalsqr: real ;
	variable factor : real ;
    	variable count: integer ;
    	variable newval: real ;
    	
     begin
    	-- check validity of argument
    	if ( x <= 0.0 ) then
       	  assert false report "X <= 0 in LOG(X)" 
			severity ERROR;
            return(REAL'LOW);
    	end if;

    	-- compute value for special cases
    	if ( X = 1.0 ) then
		return 0.0;
	else
		if ( X = MATH_E ) then
			return 1.0;
		end if;
	end if;

    	-- compute value for general cases
    	xlocal := (X - 1.0)/(X + 1.0);
    	oldval := xlocal;
    	xlocalsqr := xlocal*xlocal;
	factor := xlocal*xlocalsqr; 
    	count := 3;
    	newval := oldval + (factor/real(count));

	while ( abs(newval - oldval) > eps ) loop
		oldval := newval;
        	count := count +2;
		factor := factor * xlocalsqr;
		newval := oldval + factor/real(count);
    	end loop;

	newval := newval * 2.0;
    	return newval;
    end LOG;

    function LOG (BASE: positive; X : real) return real is
    	-- returns logarithm base BASE of X; X > 0
    begin
    	-- check validity of argument
    	if ( BASE <= 0 ) or ( x <= 0.0 ) then
       	  assert false report "BASE <= 0 or X <= 0.0 in LOG(BASE, X)" 
				severity ERROR;
            return(REAL'LOW);
    	end if;

	-- compute the value
	return ( LOG(X)/LOG(REAL(BASE)));
    end LOG;

    function  SIN (X : real ) return real is
    	-- returns sin X; X in radians
        variable n : INTEGER;
    begin 
      if (x < 1.6 ) and (x > -1.6) then
         return    CORDIC( KC, 0.0, x, 27, ROTATION)(1);
      end if; 

      n := INTEGER( x / HALF_PI );
      case QUADRANT( n mod 4 ) is 
         when 0 => 
            return  CORDIC( KC, 0.0, x - REAL(n) * HALF_PI, 27, ROTATION)(1);
         when 1 => 
            return  CORDIC( KC, 0.0, x - REAL(n) * HALF_PI, 27, ROTATION)(0);
         when 2 =>
            return -CORDIC( KC, 0.0, x - REAL(n) * HALF_PI, 27, ROTATION)(1);
         when 3 => 
            return -CORDIC( KC, 0.0, x - REAL(n) * HALF_PI, 27, ROTATION)(0);
       end case;
    end SIN;

   
   function COS (x : REAL) return REAL is 
    	-- returns cos X; X in radians
      variable n : INTEGER;
   begin 
      if (x < 1.6 ) and (x > -1.6) then
         return CORDIC( KC, 0.0, x, 27, ROTATION)(0);
      end if; 

      n := INTEGER( x / HALF_PI );
      case QUADRANT( n mod 4 ) is 
         when 0 => 
            return  CORDIC( KC, 0.0, x - REAL(n) * HALF_PI, 27, ROTATION)(0);
         when 1 => 
            return -CORDIC( KC, 0.0, x - REAL(n) * HALF_PI, 27, ROTATION)(1);
         when 2 =>
            return -CORDIC( KC, 0.0, x - REAL(n) * HALF_PI, 27, ROTATION)(0);
         when 3 => 
            return  CORDIC( KC, 0.0, x - REAL(n) * HALF_PI, 27, ROTATION)(1);
       end case;
   end COS;
   
   function TAN (x : REAL) return REAL is 
    	-- returns tan X; X in radians
    	-- X /= ((2k+1) * PI/2), where k is an integer
      variable n : INTEGER := INTEGER( x / HALF_PI );
      variable v : REAL_ARR_3 :=
                       CORDIC( KC, 0.0, x - REAL(n) * HALF_PI, 27, ROTATION);
   begin
      if n mod 2 = 0 then
         return v(1)/v(0);
      else
         return -v(0)/v(1);
      end if;
   end TAN; 
    
   function ASIN (x : real ) return real is
    	-- returns  -PI/2 < asin X < PI/2; | X | <= 1
   begin   
      if abs x > 1.0 then 
         assert false report "Out of range parameter passed to ASIN" 
			severity ERROR;
         return x;
      elsif abs x < 0.9 then 
         return atan(x/(sqrt(1.0 - x*x)));
      elsif x > 0.0 then 
         return HALF_PI - atan(sqrt(1.0 - x*x)/x);
      else 
         return - HALF_PI + atan((sqrt(1.0 - x*x))/x);
      end if;
   end ASIN; 
   
   function ACOS (x : REAL) return REAL is
    	-- returns  0 < acos X < PI; | X | <= 1
   begin  
      if abs x > 1.0 then 
         assert false report "Out of range parameter passed to ACOS" 
			severity ERROR; 
         return x;
      elsif abs x > 0.9 then
         if x > 0.0 then  
            return atan(sqrt(1.0 - x*x)/x);
         else
            return MATH_PI - atan(sqrt(1.0 - x*x)/x); 
         end if; 
      else 
         return HALF_PI - atan(x/sqrt(1.0 - x*x));
      end if;
   end ACOS; 
   
   function ATAN (x : REAL) return REAL is
    	-- returns  -PI/2 < atan X < PI/2
   begin
      return  CORDIC( 1.0, x, 0.0, 27, VECTORING )(2);      
   end ATAN; 

   function ATAN2 (x : REAL; y : REAL) return REAL is 
    	-- returns  atan (X/Y); -PI < atan2(X,Y) < PI; Y /= 0.0
   begin   
     if y = 0.0 then 
        if x = 0.0 then 
           assert false report "atan2(0.0, 0.0) is undetermined, returned 0,0" 
           	severity NOTE;
           return 0.0; 
        elsif x > 0.0 then 
           return 0.0;
        else 
           return MATH_PI;
        end if;
     elsif x > 0.0 then  
        return  CORDIC( x, y, 0.0, 27, VECTORING )(2); 
     else 
        return MATH_PI + CORDIC( x, y, 0.0, 27, VECTORING )(2); 
     end if;     
   end ATAN2; 


    function SINH (X : real) return real is
    	-- hyperbolic sine; returns (e**X - e**(-X))/2
    begin
		return ( (EXP(X) - EXP(-X))/2.0 );
    end SINH;

    function  COSH (X : real) return real is
    	-- hyperbolic cosine; returns (e**X + e**(-X))/2
    begin
		return ( (EXP(X) + EXP(-X))/2.0 );
    end COSH;

    function  TANH (X : real) return real is
    	-- hyperbolic tangent; -- returns (e**X - e**(-X))/(e**X + e**(-X))
    begin
		return ( (EXP(X) - EXP(-X))/(EXP(X) + EXP(-X)) );
    end TANH;
    
    function ASINH (X : real) return real is
    	-- returns ln( X + sqrt( X**2 + 1))
    begin
		return ( LOG( X + SQRT( X**2 + 1.0)) );
    end ASINH;

    function ACOSH (X : real) return real is
    	-- returns ln( X + sqrt( X**2 - 1));   X >= 1
    begin
      	if abs x >= 1.0 then 
         	assert false report "Out of range parameter passed to ACOSH" 
			severity ERROR; 
         	return x;
      	end if;

	return ( LOG( X + SQRT( X**2 - 1.0)) );
    end ACOSH;

    function ATANH (X : real) return real is
    	-- returns (ln( (1 + X)/(1 - X)))/2 ; | X | < 1
    begin
      	if abs x < 1.0 then 
        	assert false report "Out of range parameter passed to ATANH" 
			severity ERROR; 
        	return x;
      	end if;

	return( LOG( (1.0+X)/(1.0-X) )/2.0 );
    end ATANH;

end  MATH_REAL;



--------------------------------------------------------------- 
--
-- This source file may be used and distributed without restriction.
-- No declarations or definitions shall be included in this package.
-- This package cannot be sold or distributed for profit. 
--
--   ****************************************************************
--   *                                                              *
--   *                      W A R N I N G		 	    *
--   *								    *
--   *   This DRAFT version IS NOT endorsed or approved by IEEE     *
--   *								    *
--   ****************************************************************
--
-- Title:    PACKAGE BODY MATH_COMPLEX
--
-- Purpose:  VHDL declarations for mathematical package MATH_COMPLEX