------------------------------------------------------------------------ -- -- 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;