------------------------------------------------------------------------
--
-- This source file may be used and distributed without restriction.
-- No declarations or definitions shall be added to 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 MATH_REAL
--
-- Library: This package shall be compiled into a library
-- symbolically named IEEE.
--
-- Purpose: VHDL declarations for mathematical package MATH_REAL
-- which contains common real constants, common real
-- functions, and real trascendental functions.
--
-- Author: IEEE VHDL Math Package Study Group
--
-- Notes:
-- The package body shall be considered the formal definition of
-- the semantics of this package. Tool developers may choose to implement
-- the package body in the most efficient manner available to them.
--
-- History:
-- Version 0.1 (Strawman) Jose A. Torres 6/22/92
-- Version 0.2 Jose A. Torres 1/15/93
-- Version 0.3 Jose A. Torres 4/13/93
-- Version 0.4 Jose A. Torres 4/19/93
-- Version 0.5 Jose A. Torres 4/20/93 Added RANDOM()
-- Version 0.6 Jose A. Torres 4/23/93 Renamed RANDOM as
-- UNIFORM. Modified
-- rights banner.
-- Version 0.7 Jose A. Torres 5/28/93 Rev up for compatibility
-- with package body.
-------------------------------------------------------------
Library IEEE;
Package MATH_REAL is
--
-- commonly used constants
--
constant MATH_E : real := 2.71828_18284_59045_23536;
-- value of e
constant MATH_1_E: real := 0.36787_94411_71442_32160;
-- value of 1/e
constant MATH_PI : real := 3.14159_26535_89793_23846;
-- value of pi
constant MATH_1_PI : real := 0.31830_98861_83790_67154;
-- value of 1/pi
constant MATH_LOG_OF_2: real := 0.69314_71805_59945_30942;
-- natural log of 2
constant MATH_LOG_OF_10: real := 2.30258_50929_94045_68402;
-- natural log of10
constant MATH_LOG2_OF_E: real := 1.44269_50408_88963_4074;
-- log base 2 of e
constant MATH_LOG10_OF_E: real := 0.43429_44819_03251_82765;
-- log base 10 of e
constant MATH_SQRT2: real := 1.41421_35623_73095_04880;
-- sqrt of 2
constant MATH_SQRT1_2: real := 0.70710_67811_86547_52440;
-- sqrt of 1/2
constant MATH_SQRT_PI: real := 1.77245_38509_05516_02730;
-- sqrt of pi
constant MATH_DEG_TO_RAD: real := 0.01745_32925_19943_29577;
-- conversion factor from degree to radian
constant MATH_RAD_TO_DEG: real := 57.29577_95130_82320_87685;
-- conversion factor from radian to degree
--
-- attribute for functions whose implementation is foreign (C native)
--
attribute FOREIGN : string; -- predefined attribute in VHDL-1992
--
-- function declarations
--
function SIGN (X: real ) return real;
-- returns 1.0 if X > 0.0; 0.0 if X == 0.0; -1.0 if X < 0.0
function CEIL (X : real ) return real;
-- returns smallest integer value (as real) not less than X
function FLOOR (X : real ) return real;
-- returns largest integer value (as real) not greater than X
function ROUND (X : real ) return real;
-- returns integer FLOOR(X + 0.5) if X > 0;
-- return integer CEIL(X - 0.5) if X < 0
function FMAX (X, Y : real ) return real;
-- returns the algebraically larger of X and Y
function FMIN (X, Y : real ) return real;
-- returns the algebraically smaller of X and Y
procedure UNIFORM (variable Seed1,Seed2:inout integer; variable X:out real);
-- returns a pseudo-random number with uniform distribution in the
-- interval (0.0, 1.0).
-- Before the first call to UNIFORM, the seed values (Seed1, Seed2) must
-- be initialized to values in the range [1, 2147483562] and
-- [1, 2147483398] respectively. The seed values are modified after
-- each call to UNIFORM.
-- This random number generator is portable for 32-bit computers, and
-- it has period ~2.30584*(10**18) for each set of seed values.
--
-- For VHDL-1992, the seeds will be global variables, functions to
-- initialize their values (INIT_SEED) will be provided, and the UNIFORM
-- procedure call will be modified accordingly.
function SRAND (seed: in integer ) return integer;
--
-- sets value of seed for sequence of
-- pseudo-random numbers.
-- It uses the foreign native C function srand().
attribute FOREIGN of SRAND : function is "C_NATIVE";
function RAND return integer;
--
-- returns an integer pseudo-random number with uniform distribution.
-- It uses the foreign native C function rand().
-- Seed for the sequence is initialized with the
-- SRAND() function and value of the seed is changed every
-- time SRAND() is called, but it is not visible.
-- The range of generated values is platform dependent.
attribute FOREIGN of RAND : function is "C_NATIVE";
function GET_RAND_MAX return integer;
--
-- returns the upper bound of the range of the
-- pseudo-random numbers generated by RAND().
-- The support for this function is platform dependent, and
-- it uses foreign native C functions or constants.
-- It may not be available in some platforms.
-- Note: the value of (RAND() / GET_RAND_MAX()) is a
-- pseudo-random number distributed between 0 & 1.
attribute FOREIGN of GET_RAND_MAX : function is "C_NATIVE";
function SQRT (X : real ) return real;
-- returns square root of X; X >= 0
function CBRT (X : real ) return real;
-- returns cube root of X
function "**" (X : integer; Y : real) return real;
-- 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
function "**" (X : real; Y : real) return real;
-- 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
function EXP (X : real ) return real;
-- returns e**X; where e = MATH_E
function LOG (X : real ) return real;
-- returns natural logarithm of X; X > 0
function LOG (BASE: positive; X : real) return real;
-- returns logarithm base BASE of X; X > 0
function SIN (X : real ) return real;
-- returns sin X; X in radians
function COS ( X : real ) return real;
-- returns cos X; X in radians
function TAN (X : real ) return real;
-- returns tan X; X in radians
-- X /= ((2k+1) * PI/2), where k is an integer
function ASIN (X : real ) return real;
-- returns -PI/2 < asin X < PI/2; | X | <= 1
function ACOS (X : real ) return real;
-- returns 0 < acos X < PI; | X | <= 1
function ATAN (X : real) return real;
-- returns -PI/2 < atan X < PI/2
function ATAN2 (X : real; Y : real) return real;
-- returns atan (X/Y); -PI < atan2(X,Y) < PI; Y /= 0.0
function SINH (X : real) return real;
-- hyperbolic sine; returns (e**X - e**(-X))/2
function COSH (X : real) return real;
-- hyperbolic cosine; returns (e**X + e**(-X))/2
function TANH (X : real) return real;
-- hyperbolic tangent; -- returns (e**X - e**(-X))/(e**X + e**(-X))
function ASINH (X : real) return real;
-- returns ln( X + sqrt( X**2 + 1))
function ACOSH (X : real) return real;
-- returns ln( X + sqrt( X**2 - 1)); X >= 1
function ATANH (X : real) return real;
-- returns (ln( (1 + X)/(1 - X)))/2 ; | X | < 1
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 MATH_COMPLEX
--
-- Purpose: VHDL declarations for mathematical package MATH_COMPLEX
-- which contains common complex constants and basic complex
-- functions and operations.
--
-- Author: IEEE VHDL Math Package Study Group
--
-- Notes:
-- The package body uses package IEEE.MATH_REAL
--
-- The package body shall be considered the formal definition of
-- the semantics of this package. Tool developers may choose to implement
-- the package body in the most efficient manner available to them.
--
-- History:
-- Version 0.1 (Strawman) Jose A. Torres 6/22/92
-- Version 0.2 Jose A. Torres 1/15/93
-- Version 0.3 Jose A. Torres 4/13/93
-- Version 0.4 Jose A. Torres 4/19/93
-- Version 0.5 Jose A. Torres 4/20/93
-- Version 0.6 Jose A. Torres 4/23/93 Added unary minus
-- and CONJ for polar
-- Version 0.7 Jose A. Torres 5/28/93 Rev up for compatibility
-- with package body.
-------------------------------------------------------------
Library IEEE;
Package MATH_COMPLEX is
type COMPLEX is record RE, IM: real; end record;
type COMPLEX_VECTOR is array (integer range <>) of COMPLEX;
type COMPLEX_POLAR is record MAG: real; ARG: real; end record;
constant CBASE_1: complex := COMPLEX'(1.0, 0.0);
constant CBASE_j: complex := COMPLEX'(0.0, 1.0);
constant CZERO: complex := COMPLEX'(0.0, 0.0);
function CABS(Z: in complex ) return real;
-- returns absolute value (magnitude) of Z
function CARG(Z: in complex ) return real;
-- returns argument (angle) in radians of a complex number
function CMPLX(X: in real; Y: in real:= 0.0 ) return complex;
-- returns complex number X + iY
function "-" (Z: in complex ) return complex;
-- unary minus
function "-" (Z: in complex_polar ) return complex_polar;
-- unary minus
function CONJ (Z: in complex) return complex;
-- returns complex conjugate
function CONJ (Z: in complex_polar) return complex_polar;
-- returns complex conjugate
function CSQRT(Z: in complex ) return complex_vector;
-- returns square root of Z; 2 values
function CEXP(Z: in complex ) return complex;
-- returns e**Z
function COMPLEX_TO_POLAR(Z: in complex ) return complex_polar;
-- converts complex to complex_polar
function POLAR_TO_COMPLEX(Z: in complex_polar ) return complex;
-- converts complex_polar to complex
-- arithmetic operators
function "+" ( L: in complex; R: in complex ) return complex;
function "+" ( L: in complex_polar; R: in complex_polar) return complex;
function "+" ( L: in complex_polar; R: in complex ) return complex;
function "+" ( L: in complex; R: in complex_polar) return complex;
function "+" ( L: in real; R: in complex ) return complex;
function "+" ( L: in complex; R: in real ) return complex;
function "+" ( L: in real; R: in complex_polar) return complex;
function "+" ( L: in complex_polar; R: in real) return complex;
function "-" ( L: in complex; R: in complex ) return complex;
function "-" ( L: in complex_polar; R: in complex_polar) return complex;
function "-" ( L: in complex_polar; R: in complex ) return complex;
function "-" ( L: in complex; R: in complex_polar) return complex;
function "-" ( L: in real; R: in complex ) return complex;
function "-" ( L: in complex; R: in real ) return complex;
function "-" ( L: in real; R: in complex_polar) return complex;
function "-" ( L: in complex_polar; R: in real) return complex;
function "*" ( L: in complex; R: in complex ) return complex;
function "*" ( L: in complex_polar; R: in complex_polar) return complex;
function "*" ( L: in complex_polar; R: in complex ) return complex;
function "*" ( L: in complex; R: in complex_polar) return complex;
function "*" ( L: in real; R: in complex ) return complex;
function "*" ( L: in complex; R: in real ) return complex;
function "*" ( L: in real; R: in complex_polar) return complex;
function "*" ( L: in complex_polar; R: in real) return complex;
function "/" ( L: in complex; R: in complex ) return complex;
function "/" ( L: in complex_polar; R: in complex_polar) return complex;
function "/" ( L: in complex_polar; R: in complex ) return complex;
function "/" ( L: in complex; R: in complex_polar) return complex;
function "/" ( L: in real; R: in complex ) return complex;
function "/" ( L: in complex; R: in real ) return complex;
function "/" ( L: in real; R: in complex_polar) return complex;
function "/" ( L: in complex_polar; R: in real) return complex;
end MATH_COMPLEX;
---------------------------------------------------------------
--
-- This source file may be used and distributed without restriction.
-- No declarations or definitions shall be added to 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_REAL
--
-- Library: This package shall be compiled into a library
-- symbolically named IEEE.
--
-- Purpose: VHDL declarations for mathematical package MATH_REAL
-- which contains common real constants, common real
-- functions, and real trascendental functions.
--
-- Author: IEEE VHDL Math Package Study Group
--
-- Notes:
-- The package body shall be considered the formal definition of
-- the semantics of this package. Tool developers may choose to implement
-- the package body in the most efficient manner available to them.
--
-- Source code and algorithms for this package body comes from the
-- following sources:
-- IEEE VHDL Math Package Study Group participants,
-- U. of Mississippi, Mentor Graphics, Synopsys,
-- Viewlogic/Vantage, Communications of the ACM (June 1988, Vol
-- 31, Number 6, pp. 747, Pierre L'Ecuyer, Efficient and Portable
-- Random Number Generators), Handbook of Mathematical Functions
-- by Milton Abramowitz and Irene A. Stegun (Dover).
--
-- History:
-- Version 0.1 Jose A. Torres 4/23/93 First draft
-- Version 0.2 Jose A. Torres 5/28/93 Fixed potentially illegal code
-------------------------------------------------------------
Library IEEE;
Package body MATH_REAL is
--
-- some constants for use in the package body only
--
constant Q_PI : real := MATH_PI/4.0;
constant HALF_PI : real := MATH_PI/2.0;
constant TWO_PI : real := MATH_PI*2.0;
constant MAX_ITER: integer := 27; -- max precision factor for cordic
--
-- some type declarations for cordic operations
--
constant KC : REAL := 6.0725293500888142e-01; -- constant for cordic
type REAL_VECTOR is array (NATURAL range <>) of REAL;
type NATURAL_VECTOR is array (NATURAL range <>) of NATURAL;
subtype REAL_VECTOR_N is REAL_VECTOR (0 to max_iter);
subtype REAL_ARR_2 is REAL_VECTOR (0 to 1);
subtype REAL_ARR_3 is REAL_VECTOR (0 to 2);
subtype QUADRANT is INTEGER range 0 to 3;
type CORDIC_MODE_TYPE is (ROTATION, VECTORING);
--
-- auxiliary functions for cordic algorithms
--
function POWER_OF_2_SERIES (d : NATURAL_VECTOR; initial_value : REAL;
number_of_values : NATURAL) return REAL_VECTOR is
variable v : REAL_VECTOR (0 to number_of_values);
variable temp : REAL := initial_value;
variable flag : boolean := true;
begin
for i in 0 to number_of_values loop
v(i) := temp;
for p in d'range loop
if i = d(p) then
flag := false;
end if;
end loop;
if flag then
temp := temp/2.0;
end if;
flag := true;
end loop;
return v;
end POWER_OF_2_SERIES;
constant two_at_minus : REAL_VECTOR := POWER_OF_2_SERIES(
NATURAL_VECTOR'(100, 90),1.0,
MAX_ITER);
constant epsilon : REAL_VECTOR_N := (
7.8539816339744827e-01,
4.6364760900080606e-01,
2.4497866312686413e-01,
1.2435499454676144e-01,
6.2418809995957351e-02,
3.1239833430268277e-02,
1.5623728620476830e-02,
7.8123410601011116e-03,
3.9062301319669717e-03,
1.9531225164788189e-03,
9.7656218955931937e-04,
4.8828121119489829e-04,
2.4414062014936175e-04,
1.2207031189367021e-04,
6.1035156174208768e-05,
3.0517578115526093e-05,
1.5258789061315760e-05,
7.6293945311019699e-06,
3.8146972656064960e-06,
1.9073486328101870e-06,
9.5367431640596080e-07,
4.7683715820308876e-07,
2.3841857910155801e-07,
1.1920928955078067e-07,
5.9604644775390553e-08,
2.9802322387695303e-08,
1.4901161193847654e-08,
7.4505805969238281e-09
);
function CORDIC ( x0 : REAL;
y0 : REAL;
z0 : REAL;
n : NATURAL; -- precision factor
CORDIC_MODE : CORDIC_MODE_TYPE -- rotation (z -> 0)
-- or vectoring (y -> 0)
) return REAL_ARR_3 is
variable x : REAL := x0;
variable y : REAL := y0;
variable z : REAL := z0;
variable x_temp : REAL;
begin
if CORDIC_MODE = ROTATION then
for k in 0 to n loop
x_temp := x;
if ( z >= 0.0) then
x := x - y * two_at_minus(k);
y := y + x_temp * two_at_minus(k);
z := z - epsilon(k);
else
x := x + y * two_at_minus(k);
y := y - x_temp * two_at_minus(k);
z := z + epsilon(k);
end if;
end loop;
else
for k in 0 to n loop
x_temp := x;
if ( y < 0.0) then
x := x - y * two_at_minus(k);
y := y + x_temp * two_at_minus(k);
z := z - epsilon(k);
else
x := x + y * two_at_minus(k);
y := y - x_temp * two_at_minus(k);
z := z + epsilon(k);
end if;
end loop;
end if;
return REAL_ARR_3'(x, y, z);
end CORDIC;
--
-- non-trascendental functions
--
function SIGN (X: real ) return real is
-- returns 1.0 if X > 0.0; 0.0 if X == 0.0; -1.0 if X < 0.0
begin
if ( X > 0.0 ) then
return 1.0;
elsif ( X < 0.0 ) then
return -1.0;
else
return 0.0;
end if;
end SIGN;
function CEIL (X : real ) return real is
-- returns smallest integer value (as real) not less than X
-- No conversion to an integer type is expected, so truncate cannot
-- overflow for large arguments.
variable large: real := 1073741824.0;
type long is range -1073741824 to 1073741824;
-- 2**30 is longer than any single-precision mantissa
variable rd: real;
begin
if abs( X) >= large then
return X;
else
rd := real ( long( X));
if X > 0.0 then
if rd >= X then
return rd;
else
return rd + 1.0;
end if;
elsif X = 0.0 then
return 0.0;
else
if rd <= X then
return rd;
else
return rd - 1.0;
end if;
end if;
end if;
end CEIL;
function FLOOR (X : real ) return real is
-- returns largest integer value (as real) not greater than X
-- No conversion to an integer type is expected, so truncate
-- cannot overflow for large arguments.
--
variable large: real := 1073741824.0;
type long is range -1073741824 to 1073741824;
-- 2**30 is longer than any single-precision mantissa
variable rd: real;
begin
if abs( X ) >= large then
return X;
else
rd := real ( long( X));
if X > 0.0 then
if rd <= X then
return rd;
else
return rd - 1.0;
end if;
elsif X = 0.0 then
return 0.0;
else
if rd >= X then
return rd;
else
return rd + 1.0;
end if;
end if;
end if;
end FLOOR;
function ROUND (X : real ) return real is
-- returns integer FLOOR(X + 0.5) if X > 0;
-- return integer CEIL(X - 0.5) if X < 0
begin
if X > 0.0 then
return FLOOR(X + 0.5);
elsif X < 0.0 then
return CEIL( X - 0.5);
else
return 0.0;
end if;
end ROUND;
function FMAX (X, Y : real ) return real is
-- returns the algebraically larger of X and Y
begin
if X > Y then
return X;
else
return Y;
end if;
end FMAX;
function FMIN (X, Y : real ) return real is
-- returns the algebraically smaller of X and Y
begin
if X < Y then
return X;
else
return Y;
end if;
end FMIN;
--
-- Pseudo-random number generators
--
procedure UNIFORM(variable Seed1,Seed2:inout integer;variable X:out real) is
-- returns a pseudo-random number with uniform distribution in the
-- interval (0.0, 1.0).
-- Before the first call to UNIFORM, the seed values (Seed1, Seed2) must
-- be initialized to values in the range [1, 2147483562] and
-- [1, 2147483398] respectively. The seed values are modified after
-- each call to UNIFORM.
-- This random number generator is portable for 32-bit computers, and
-- it has period ~2.30584*(10**18) for each set of seed values.
--
-- For VHDL-1992, the seeds will be global variables, functions to
-- initialize their values (INIT_SEED) will be provided, and the UNIFORM
-- procedure call will be modified accordingly.
variable z, k: integer;
begin
k := Seed1/53668;
Seed1 := 40014 * (Seed1 - k * 53668) - k * 12211;
if Seed1 < 0 then
Seed1 := Seed1 + 2147483563;
end if;
k := Seed2/52774;
Seed2 := 40692 * (Seed2 - k * 52774) - k * 3791;
if Seed2 < 0 then
Seed2 := Seed2 + 2147483399;
end if;
z := Seed1 - Seed2;
if z < 1 then
z := z + 2147483562;
end if;
X := REAL(Z)*4.656613e-10;
end UNIFORM;
function SRAND (seed: in integer ) return integer is
--
-- sets value of seed for sequence of
-- pseudo-random numbers.
-- Returns the value of the seed.
-- It uses the foreign native C function srand().
begin
end SRAND;
function RAND return integer is
--
-- returns an integer pseudo-random number with uniform distribution.
-- It uses the foreign native C function rand().
-- Seed for the sequence is initialized with the
-- SRAND() function and value of the seed is changed every
-- time SRAND() is called, but it is not visible.
-- The range of generated values is platform dependent.
begin
end RAND;
function GET_RAND_MAX return integer is
--
-- returns the upper bound of the range of the
-- pseudo-random numbers generated by RAND().
-- The support for this function is platform dependent, and
-- it uses foreign native C functions or constants.
-- It may not be available in some platforms.
-- Note: the value of (RAND / GET_RAND_MAX) is a
-- pseudo-random number distributed between 0 & 1.
begin
end GET_RAND_MAX;
--
-- trascendental and trigonometric functions
--
function SQRT (X : real ) return real is
-- returns square root of X; X >= 0
--
-- Computes square root using the Newton-Raphson approximation:
-- F(n+1) = 0.5*[F(n) + x/F(n)];
--
constant inival: real := 1.5;
constant eps : real := 0.000001;
constant relative_err : real := eps*X;
variable oldval : real ;
variable newval : real ;
begin
-- check validity of argument
if ( X < 0.0 ) then
assert false report "X < 0 in SQRT(X)"
severity ERROR;
return (0.0);
end if;
-- get the square root for special cases
if X = 0.0 then
return 0.0;
else
if ( X = 1.0 ) then
return 1.0; -- return exact value
end if;
end if;
-- get the square root for general cases
oldval := inival;
newval := (X/oldval + oldval)/2.0;
while ( abs(newval -oldval) > relative_err ) loop
oldval := newval;
newval := (X/oldval + oldval)/2.0;
end loop;
return newval;
end SQRT;
function CBRT (X : real ) return real is
-- returns cube root of X
-- Computes square root using the Newton-Raphson approximation:
-- F(n+1) = (1/3)*[2*F(n) + x/F(n)**2];
--
constant inival: real := 1.5;
constant eps : real := 0.000001;
constant relative_err : real := eps*abs(X);
variable xlocal : real := X;
variable negative : boolean := X < 0.0;
variable oldval : real ;
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
-- which contains common complex constants and basic complex
-- functions and operations.
--
-- Author: IEEE VHDL Math Package Study Group
--
-- Notes:
-- The package body uses package IEEE.MATH_REAL
--
-- The package body shall be considered the formal definition of
-- the semantics of this package. Tool developers may choose to implement
-- the package body in the most efficient manner available to them.
--
-- Source code for this package body comes from the following
-- following sources:
-- IEEE VHDL Math Package Study Group participants,
-- U. of Mississippi, Mentor Graphics, Synopsys,
-- Viewlogic/Vantage, Communications of the ACM (June 1988, Vol
-- 31, Number 6, pp. 747, Pierre L'Ecuyer, Efficient and Portable
-- Random Number Generators, Handbook of Mathematical Functions
-- by Milton Abramowitz and Irene A. Stegun (Dover).
--
-- History:
-- Version 0.1 Jose A. Torres 4/23/93 First draft
-- Version 0.2 Jose A. Torres 5/28/93 Fixed potentially illegal code
--
-------------------------------------------------------------
Library IEEE;
Use IEEE.MATH_REAL.all; -- real trascendental operations
Package body MATH_COMPLEX is
function CABS(Z: in complex ) return real is
-- returns absolute value (magnitude) of Z
variable ztemp : complex_polar;
begin
ztemp := COMPLEX_TO_POLAR(Z);
return ztemp.mag;
end CABS;
function CARG(Z: in complex ) return real is
-- returns argument (angle) in radians of a complex number
variable ztemp : complex_polar;
begin
ztemp := COMPLEX_TO_POLAR(Z);
return ztemp.arg;
end CARG;
function CMPLX(X: in real; Y: in real := 0.0 ) return complex is
-- returns complex number X + iY
begin
return COMPLEX'(X, Y);
end CMPLX;
function "-" (Z: in complex ) return complex is
-- unary minus; returns -x -jy for z= x + jy
begin
return COMPLEX'(-z.Re, -z.Im);
end "-";
function "-" (Z: in complex_polar ) return complex_polar is
-- unary minus; returns (z.mag, z.arg + MATH_PI)
begin
return COMPLEX_POLAR'(z.mag, z.arg + MATH_PI);
end "-";
function CONJ (Z: in complex) return complex is
-- returns complex conjugate (x-jy for z = x+ jy)
begin
return COMPLEX'(z.Re, -z.Im);
end CONJ;
function CONJ (Z: in complex_polar) return complex_polar is
-- returns complex conjugate (z.mag, -z.arg)
begin
return COMPLEX_POLAR'(z.mag, -z.arg);
end CONJ;
function CSQRT(Z: in complex ) return complex_vector is
-- returns square root of Z; 2 values
variable ztemp : complex_polar;
variable zout : complex_vector (0 to 1);
variable temp : real;
begin
ztemp := COMPLEX_TO_POLAR(Z);
temp := SQRT(ztemp.mag);
zout(0).re := temp*COS(ztemp.arg/2.0);
zout(0).im := temp*SIN(ztemp.arg/2.0);
zout(1).re := temp*COS(ztemp.arg/2.0 + MATH_PI);
zout(1).im := temp*SIN(ztemp.arg/2.0 + MATH_PI);
return zout;
end CSQRT;
function CEXP(Z: in complex ) return complex is
-- returns e**Z
begin
return COMPLEX'(EXP(Z.re)*COS(Z.im), EXP(Z.re)*SIN(Z.im));
end CEXP;
function COMPLEX_TO_POLAR(Z: in complex ) return complex_polar is
-- converts complex to complex_polar
begin
return COMPLEX_POLAR'(sqrt(z.re**2 + z.im**2),atan2(z.re,z.im));
end COMPLEX_TO_POLAR;
function POLAR_TO_COMPLEX(Z: in complex_polar ) return complex is
-- converts complex_polar to complex
begin
return COMPLEX'( z.mag*cos(z.arg), z.mag*sin(z.arg) );
end POLAR_TO_COMPLEX;
--
-- arithmetic operators
--
function "+" ( L: in complex; R: in complex ) return complex is
begin
return COMPLEX'(L.Re + R.Re, L.Im + R.Im);
end "+";
function "+" (L: in complex_polar; R: in complex_polar) return complex is
variable zL, zR : complex;
begin
zL := POLAR_TO_COMPLEX( L );
zR := POLAR_TO_COMPLEX( R );
return COMPLEX'(zL.Re + zR.Re, zL.Im + zR.Im);
end "+";
function "+" ( L: in complex_polar; R: in complex ) return complex is
variable zL : complex;
begin
zL := POLAR_TO_COMPLEX( L );
return COMPLEX'(zL.Re + R.Re, zL.Im + R.Im);
end "+";
function "+" ( L: in complex; R: in complex_polar) return complex is
variable zR : complex;
begin
zR := POLAR_TO_COMPLEX( R );
return COMPLEX'(L.Re + zR.Re, L.Im + zR.Im);
end "+";
function "+" ( L: in real; R: in complex ) return complex is
begin
return COMPLEX'(L + R.Re, R.Im);
end "+";
function "+" ( L: in complex; R: in real ) return complex is
begin
return COMPLEX'(L.Re + R, L.Im);
end "+";
function "+" ( L: in real; R: in complex_polar) return complex is
variable zR : complex;
begin
zR := POLAR_TO_COMPLEX( R );
return COMPLEX'(L + zR.Re, zR.Im);
end "+";
function "+" ( L: in complex_polar; R: in real) return complex is
variable zL : complex;
begin
zL := POLAR_TO_COMPLEX( L );
return COMPLEX'(zL.Re + R, zL.Im);
end "+";
function "-" ( L: in complex; R: in complex ) return complex is
begin
return COMPLEX'(L.Re - R.Re, L.Im - R.Im);
end "-";
function "-" ( L: in complex_polar; R: in complex_polar) return complex is
variable zL, zR : complex;
begin
zL := POLAR_TO_COMPLEX( L );
zR := POLAR_TO_COMPLEX( R );
return COMPLEX'(zL.Re - zR.Re, zL.Im - zR.Im);
end "-";
function "-" ( L: in complex_polar; R: in complex ) return complex is
variable zL : complex;
begin
zL := POLAR_TO_COMPLEX( L );
return COMPLEX'(zL.Re - R.Re, zL.Im - R.Im);
end "-";
function "-" ( L: in complex; R: in complex_polar) return complex is
variable zR : complex;
begin
zR := POLAR_TO_COMPLEX( R );
return COMPLEX'(L.Re - zR.Re, L.Im - zR.Im);
end "-";
function "-" ( L: in real; R: in complex ) return complex is
begin
return COMPLEX'(L - R.Re, -1.0 * R.Im);
end "-";
function "-" ( L: in complex; R: in real ) return complex is
begin
return COMPLEX'(L.Re - R, L.Im);
end "-";
function "-" ( L: in real; R: in complex_polar) return complex is
variable zR : complex;
begin
zR := POLAR_TO_COMPLEX( R );
return COMPLEX'(L - zR.Re, -1.0*zR.Im);
end "-";
function "-" ( L: in complex_polar; R: in real) return complex is
variable zL : complex;
begin
zL := POLAR_TO_COMPLEX( L );
return COMPLEX'(zL.Re - R, zL.Im);
end "-";
function "*" ( L: in complex; R: in complex ) return complex is
begin
return COMPLEX'(L.Re * R.Re - L.Im * R.Im, L.Re * R.Im + L.Im * R.Re);
end "*";
function "*" ( L: in complex_polar; R: in complex_polar) return complex is
variable zout : complex_polar;
begin
zout.mag := L.mag * R.mag;
zout.arg := L.arg + R.arg;
return POLAR_TO_COMPLEX(zout);
end "*";
function "*" ( L: in complex_polar; R: in complex ) return complex is
variable zL : complex;
begin
zL := POLAR_TO_COMPLEX( L );
return COMPLEX'(zL.Re*R.Re - zL.Im * R.Im, zL.Re * R.Im + zL.Im*R.Re);
end "*";
function "*" ( L: in complex; R: in complex_polar) return complex is
variable zR : complex;
begin
zR := POLAR_TO_COMPLEX( R );
return COMPLEX'(L.Re*zR.Re - L.Im * zR.Im, L.Re * zR.Im + L.Im*zR.Re);
end "*";
function "*" ( L: in real; R: in complex ) return complex is
begin
return COMPLEX'(L * R.Re, L * R.Im);
end "*";
function "*" ( L: in complex; R: in real ) return complex is
begin
return COMPLEX'(L.Re * R, L.Im * R);
end "*";
function "*" ( L: in real; R: in complex_polar) return complex is
variable zR : complex;
begin
zR := POLAR_TO_COMPLEX( R );
return COMPLEX'(L * zR.Re, L * zR.Im);
end "*";
function "*" ( L: in complex_polar; R: in real) return complex is
variable zL : complex;
begin
zL := POLAR_TO_COMPLEX( L );
return COMPLEX'(zL.Re * R, zL.Im * R);
end "*";
function "/" ( L: in complex; R: in complex ) return complex is
variable magrsq : REAL := R.Re ** 2 + R.Im ** 2;
begin
if (magrsq = 0.0) then
assert FALSE report "Attempt to divide by (0,0)"
severity ERROR;
return COMPLEX'(REAL'RIGHT, REAL'RIGHT);
else
return COMPLEX'( (L.Re * R.Re + L.Im * R.Im) / magrsq,
(L.Im * R.Re - L.Re * R.Im) / magrsq);
end if;
end "/";
function "/" ( L: in complex_polar; R: in complex_polar) return complex is
variable zout : complex_polar;
begin
if (R.mag = 0.0) then
assert FALSE report "Attempt to divide by (0,0)"
severity ERROR;
return COMPLEX'(REAL'RIGHT, REAL'RIGHT);
else
zout.mag := L.mag/R.mag;
zout.arg := L.arg - R.arg;
return POLAR_TO_COMPLEX(zout);
end if;
end "/";
function "/" ( L: in complex_polar; R: in complex ) return complex is
variable zL : complex;
variable temp : REAL := R.Re ** 2 + R.Im ** 2;
begin
if (temp = 0.0) then
assert FALSE report "Attempt to divide by (0.0,0.0)"
severity ERROR;
return COMPLEX'(REAL'RIGHT, REAL'RIGHT);
else
zL := POLAR_TO_COMPLEX( L );
return COMPLEX'( (zL.Re * R.Re + zL.Im * R.Im) / temp,
(zL.Im * R.Re - zL.Re * R.Im) / temp);
end if;
end "/";
function "/" ( L: in complex; R: in complex_polar) return complex is
variable zR : complex := POLAR_TO_COMPLEX( R );
variable temp : REAL := zR.Re ** 2 + zR.Im ** 2;
begin
if (R.mag = 0.0) or (temp = 0.0) then
assert FALSE report "Attempt to divide by (0.0,0.0)"
severity ERROR;
return COMPLEX'(REAL'RIGHT, REAL'RIGHT);
else
return COMPLEX'( (L.Re * zR.Re + L.Im * zR.Im) / temp,
(L.Im * zR.Re - L.Re * zR.Im) / temp);
end if;
end "/";
function "/" ( L: in real; R: in complex ) return complex is
variable temp : REAL := R.Re ** 2 + R.Im ** 2;
begin
if (temp = 0.0) then
assert FALSE report "Attempt to divide by (0.0,0.0)"
severity ERROR;
return COMPLEX'(REAL'RIGHT, REAL'RIGHT);
else
temp := L / temp;
return COMPLEX'( temp * R.Re, -temp * R.Im );
end if;
end "/";
function "/" ( L: in complex; R: in real ) return complex is
begin
if (R = 0.0) then
assert FALSE report "Attempt to divide by (0.0,0.0)"
severity ERROR;
return COMPLEX'(REAL'RIGHT, REAL'RIGHT);
else
return COMPLEX'(L.Re / R, L.Im / R);
end if;
end "/";
function "/" ( L: in real; R: in complex_polar) return complex is
variable zR : complex := POLAR_TO_COMPLEX( R );
variable temp : REAL := zR.Re ** 2 + zR.Im ** 2;
begin
if (R.mag = 0.0) or (temp = 0.0) then
assert FALSE report "Attempt to divide by (0.0,0.0)"
severity ERROR;
return COMPLEX'(REAL'RIGHT, REAL'RIGHT);
else
temp := L / temp;
return COMPLEX'( temp * zR.Re, -temp * zR.Im );
end if;
end "/";
function "/" ( L: in complex_polar; R: in real) return complex is
variable zL : complex := POLAR_TO_COMPLEX( L );
begin
if (R = 0.0) then
assert FALSE report "Attempt to divide by (0.0,0.0)"
severity ERROR;
return COMPLEX'(REAL'RIGHT, REAL'RIGHT);
else
return COMPLEX'(zL.Re / R, zL.Im / R);
end if;
end "/";
end MATH_COMPLEX;
<div align="center"><br /><script type="text/javascript"><!--
google_ad_client = "pub-7293844627074885";
//468x60, Created at 07. 11. 25
google_ad_slot = "8619794253";
google_ad_width = 468;
google_ad_height = 60;
//--></script>
<script type="text/javascript" src="http://pagead2.googlesyndication.com/pagead/show_ads.js">
</script><br /> </div>