diff options
Diffstat (limited to 'libraries/ieee/math_real-body.vhdl')
-rw-r--r-- | libraries/ieee/math_real-body.vhdl | 2157 |
1 files changed, 1837 insertions, 320 deletions
diff --git a/libraries/ieee/math_real-body.vhdl b/libraries/ieee/math_real-body.vhdl index 41e6a084d..c20aab7c9 100644 --- a/libraries/ieee/math_real-body.vhdl +++ b/libraries/ieee/math_real-body.vhdl @@ -1,424 +1,1941 @@ ---------------------------------------------------------------- +------------------------------------------------------------------------ -- --- 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. +-- Copyright 1996 by IEEE. All rights reserved. + +-- This source file is an informative part of IEEE Std 1076.2-1996, IEEE Standard +-- VHDL Mathematical Packages. This source file may not be copied, sold, or +-- included with software that is sold without written permission from the IEEE +-- Standards Department. This source file may be used to implement this standard +-- and may be distributed in compiled form in any manner so long as the +-- compiled form does not allow direct decompilation of the original source file. +-- This source file may be copied for individual use between licensed users. +-- This source file is provided on an AS IS basis. The IEEE disclaims ANY +-- WARRANTY EXPRESS OR IMPLIED INCLUDING ANY WARRANTY OF MERCHANTABILITY +-- AND FITNESS FOR USE FOR A PARTICULAR PURPOSE. The user of the source +-- file shall indemnify and hold IEEE harmless from any damages or liability +-- arising out of the use thereof. + -- --- **************************************************************** --- * * --- * W A R N I N G * --- * * --- * This DRAFT version IS NOT endorsed or approved by IEEE * --- * * --- **************************************************************** +-- Title: Standard VHDL Mathematical Packages (IEEE Std 1076.2-1996, +-- MATH_REAL) -- --- Title: PACKAGE BODY MATH_REAL +-- Library: This package shall be compiled into a library +-- symbolically named IEEE. -- --- Library: This package shall be compiled into a library --- symbolically named IEEE. +-- Developers: IEEE DASC VHDL Mathematical Packages Working Group -- --- Purpose: VHDL declarations for mathematical package MATH_REAL --- which contains common real constants, common real --- functions, and real trascendental functions. +-- Purpose: This package body is a nonnormative implementation of the +-- functionality defined in the MATH_REAL package declaration. -- --- Author: IEEE VHDL Math Package Study Group +-- Limitation: The values generated by the functions in this package may +-- vary from platform to platform, and the precision of results +-- is only guaranteed to be the minimum required by IEEE Std 1076 +-- -1993. -- -- 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). +-- The "package declaration" defines the types, subtypes, and +-- declarations of MATH_REAL. +-- The standard mathematical definition and conventional meaning +-- of the mathematical functions that are part of this standard +-- represent the formal semantics of the implementation of the +-- MATH_REAL package declaration. The purpose of the MATH_REAL +-- package body is to clarify such semantics and provide a +-- guideline for implementations to verify their implementation +-- of MATH_REAL. Tool developers may choose to implement +-- the package body in the most efficient manner available to them. -- --- 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 --- --- GHDL history --- 2005-04-07 Initial version. --- 2005-12-23 I. Curtis : overhaul of log functions to bring in line --- with ieee standard -------------------------------------------------------------- -Library IEEE; - -Package body MATH_REAL is +-- ----------------------------------------------------------------------------- +-- Version : 1.5 +-- Date : 24 July 1996 +-- ----------------------------------------------------------------------------- + +package body MATH_REAL is + + -- + -- Local Constants for Use in the Package Body Only + -- + constant MATH_E_P2 : REAL := 7.38905_60989_30650; -- e**2 + constant MATH_E_P10 : REAL := 22026.46579_48067_17; -- e**10 + constant MATH_EIGHT_PI : REAL := 25.13274_12287_18345_90770_115; --8*pi + constant MAX_ITER: INTEGER := 27; -- Maximum precision factor for cordic + constant MAX_COUNT: INTEGER := 150; -- Maximum count for number of tries + constant BASE_EPS: REAL := 0.00001; -- Factor for convergence criteria + constant KC : REAL := 6.0725293500888142e-01; -- Constant for cordic + + -- + -- Local Type Declarations for Cordic Operations + -- + 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 : in NATURAL_VECTOR; INITIAL_VALUE : in REAL; + NUMBER_OF_VALUES : in NATURAL) return REAL_VECTOR is + -- Description: + -- Returns power of two for a vector of values + -- Notes: + -- None + -- + 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; + exit; + 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 : in REAL; + Y0 : in REAL; + Z0 : in REAL; + N : in NATURAL; -- Precision factor + CORDIC_MODE : in CORDIC_MODE_TYPE -- Rotation (Z -> 0) + -- or vectoring (Y -> 0) + ) return REAL_ARR_3 is + -- Description: + -- Compute cordic values + -- Notes: + -- None + 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 + -- Bodies for Global Mathematical Functions Start Here -- - 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 + function SIGN (X: in REAL ) return REAL is + -- Description: + -- See function declaration in IEEE Std 1076.2-1996 + -- Notes: + -- None begin - assert false severity failure; - end SIGN; + 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 : in REAL ) return REAL is + -- Description: + -- See function declaration in IEEE Std 1076.2-1996 + -- Notes: + -- a) No conversion to an INTEGER type is expected, so truncate + -- cannot overflow for large arguments + -- b) The domain supported by this function is X <= LARGE + -- c) Returns X if ABS(X) >= LARGE + + constant LARGE: REAL := REAL(INTEGER'HIGH); + variable RD: REAL; - function CEIL (X : real ) return real is begin - assert false severity failure; - end CEIL; + if ABS(X) >= LARGE then + return X; + end if; + + RD := REAL ( INTEGER(X)); + if RD = X then + return X; + end if; + + 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 + 1.0; + else + return RD; + end if; + end if; + end CEIL; + + function FLOOR (X : in REAL ) return REAL is + -- Description: + -- See function declaration in IEEE Std 1076.2-1996 + -- Notes: + -- a) No conversion to an INTEGER type is expected, so truncate + -- cannot overflow for large arguments + -- b) The domain supported by this function is ABS(X) <= LARGE + -- c) Returns X if ABS(X) >= LARGE + + constant LARGE: REAL := REAL(INTEGER'HIGH); + variable RD: REAL; - function FLOOR (X : real ) return real is begin - assert false severity failure; + if ABS( X ) >= LARGE then + return X; + end if; + + RD := REAL ( INTEGER(X)); + if RD = X then + return X; + end if; + + 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 - 1.0; + else + return RD; + end if; + end if; end FLOOR; - function ROUND (X : real ) return real is + function ROUND (X : in REAL ) return REAL is + -- Description: + -- See function declaration in IEEE Std 1076.2-1996 + -- Notes: + -- a) Returns 0.0 if X = 0.0 + -- b) Returns FLOOR(X + 0.5) if X > 0 + -- c) Returns CEIL(X - 0.5) if X < 0 + begin - assert false severity failure; + 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 - begin - assert false severity failure; - end FMAX; - function FMIN (X, Y : real ) return real is + function TRUNC (X : in REAL ) return REAL is + -- Description: + -- See function declaration in IEEE Std 1076.2-1996 + -- Notes: + -- a) Returns 0.0 if X = 0.0 + -- b) Returns FLOOR(X) if X > 0 + -- c) Returns CEIL(X) if X < 0 + begin - assert false severity failure; - end FMIN; + if X > 0.0 then + return FLOOR(X); + elsif X < 0.0 then + return CEIL( X); + else + return 0.0; + end if; + end TRUNC; - -- - -- 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 "MOD" (X, Y: in REAL ) return REAL is + -- Description: + -- See function declaration in IEEE Std 1076.2-1996 + -- Notes: + -- a) Returns 0.0 on error - function SRAND (seed: in integer ) return integer is + variable XNEGATIVE : BOOLEAN := X < 0.0; + variable YNEGATIVE : BOOLEAN := Y < 0.0; + variable VALUE : REAL; begin - assert false severity failure; - end SRAND; + -- Check validity of input arguments + if (Y = 0.0) then + assert FALSE + report "MOD(X, 0.0) is undefined" + severity ERROR; + return 0.0; + end if; + + -- Compute value + if ( XNEGATIVE ) then + if ( YNEGATIVE ) then + VALUE := X + (FLOOR(ABS(X)/ABS(Y)))*ABS(Y); + else + VALUE := X + (CEIL(ABS(X)/ABS(Y)))*ABS(Y); + end if; + else + if ( YNEGATIVE ) then + VALUE := X - (CEIL(ABS(X)/ABS(Y)))*ABS(Y); + else + VALUE := X - (FLOOR(ABS(X)/ABS(Y)))*ABS(Y); + end if; + end if; + + return VALUE; + end "MOD"; - function RAND return integer is + + function REALMAX (X, Y : in REAL ) return REAL is + -- Description: + -- See function declaration in IEEE Std 1076.2-1996 + -- Notes: + -- a) REALMAX(X,Y) = X when X = Y + -- begin - assert false severity failure; - end RAND; + if X >= Y then + return X; + else + return Y; + end if; + end REALMAX; - function GET_RAND_MAX return integer is - -- The value this function returns should be the same as - -- RAND_MAX in /usr/include/stdlib.h + function REALMIN (X, Y : in REAL ) return REAL is + -- Description: + -- See function declaration in IEEE Std 1076.2-1996 + -- Notes: + -- a) REALMIN(X,Y) = X when X = Y + -- begin - assert false - report "Be sure to update GET_RAND_MAX in mathpack.vhd" - severity note; - return 2147483647; -- i386 linux - end GET_RAND_MAX; + if X <= Y then + return X; + else + return Y; + end if; + end REALMIN; - -- - -- trascendental and trigonometric functions - -- - function c_sqrt (x : real ) return real; - attribute foreign of c_sqrt : function is "VHPIDIRECT sqrt"; - - function c_sqrt (x : real ) return real is + + procedure UNIFORM(variable SEED1,SEED2:inout POSITIVE;variable X:out REAL) + is + -- Description: + -- See function declaration in IEEE Std 1076.2-1996 + -- Notes: + -- a) Returns 0.0 on error + -- + variable Z, K: INTEGER; + variable TSEED1 : INTEGER := INTEGER'(SEED1); + variable TSEED2 : INTEGER := INTEGER'(SEED2); begin - assert false severity failure; - end c_sqrt; + -- Check validity of arguments + if SEED1 > 2147483562 then + assert FALSE + report "SEED1 > 2147483562 in UNIFORM" + severity ERROR; + X := 0.0; + return; + end if; + + if SEED2 > 2147483398 then + assert FALSE + report "SEED2 > 2147483398 in UNIFORM" + severity ERROR; + X := 0.0; + return; + end if; + + -- Compute new seed values and pseudo-random number + K := TSEED1/53668; + TSEED1 := 40014 * (TSEED1 - K * 53668) - K * 12211; + + if TSEED1 < 0 then + TSEED1 := TSEED1 + 2147483563; + end if; + + K := TSEED2/52774; + TSEED2 := 40692 * (TSEED2 - K * 52774) - K * 3791; + + if TSEED2 < 0 then + TSEED2 := TSEED2 + 2147483399; + end if; + + Z := TSEED1 - TSEED2; + if Z < 1 then + Z := Z + 2147483562; + end if; + + -- Get output values + SEED1 := POSITIVE'(TSEED1); + SEED2 := POSITIVE'(TSEED2); + X := REAL(Z)*4.656613e-10; + end UNIFORM; + + + + function SQRT (X : in REAL ) return REAL is + -- Description: + -- See function declaration in IEEE Std 1076.2-1996 + -- Notes: + -- a) Uses the Newton-Raphson approximation: + -- F(n+1) = 0.5*[F(n) + x/F(n)] + -- b) Returns 0.0 on error + -- + + constant EPS : REAL := BASE_EPS*BASE_EPS; -- Convergence factor + + variable INIVAL: REAL; + variable OLDVAL : REAL ; + variable NEWVAL : REAL ; + variable COUNT : INTEGER := 1; - function SQRT (X : real ) return real is begin - -- check validity of argument + -- Check validity of argument if ( X < 0.0 ) then - assert false report "X < 0 in SQRT(X)" - severity ERROR; - return (0.0); + assert FALSE + report "X < 0.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; + end if; end if; - return c_sqrt(X); + + -- Get the square root for general cases + INIVAL := EXP(LOG(X)*(0.5)); -- Mathematically correct but imprecise + OLDVAL := INIVAL; + NEWVAL := (X/OLDVAL + OLDVAL)*0.5; + + -- Check for relative and absolute error and max count + while ( ( (ABS((NEWVAL -OLDVAL)/NEWVAL) > EPS) OR + (ABS(NEWVAL - OLDVAL) > EPS) ) AND + (COUNT < MAX_COUNT) ) loop + OLDVAL := NEWVAL; + NEWVAL := (X/OLDVAL + OLDVAL)*0.5; + COUNT := COUNT + 1; + end loop; + return NEWVAL; end SQRT; - function CBRT (X : real ) return real is + function CBRT (X : in REAL ) return REAL is + -- Description: + -- See function declaration in IEEE Std 1076.2-1996 + -- Notes: + -- a) Uses the Newton-Raphson approximation: + -- F(n+1) = (1/3)*[2*F(n) + x/F(n)**2]; + -- + constant EPS : REAL := BASE_EPS*BASE_EPS; + + variable INIVAL: REAL; + variable XLOCAL : REAL := X; + variable NEGATIVE : BOOLEAN := X < 0.0; + variable OLDVAL : REAL ; + variable NEWVAL : REAL ; + variable COUNT : INTEGER := 1; + begin - assert false severity failure; + + -- 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; + + INIVAL := EXP(LOG(XLOCAL)/(3.0)); -- Mathematically correct but + -- imprecise + OLDVAL := INIVAL; + NEWVAL := (XLOCAL/(OLDVAL*OLDVAL) + 2.0*OLDVAL)/3.0; + + -- Check for relative and absolute errors and max count + while ( ( (ABS((NEWVAL -OLDVAL)/NEWVAL) > EPS ) OR + (ABS(NEWVAL - OLDVAL) > EPS ) ) AND + ( COUNT < MAX_COUNT ) ) loop + OLDVAL := NEWVAL; + NEWVAL :=(XLOCAL/(OLDVAL*OLDVAL) + 2.0*OLDVAL)/3.0; + COUNT := COUNT + 1; + 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 + function "**" (X : in INTEGER; Y : in REAL) return REAL is + -- Description: + -- See function declaration in IEEE Std 1076.2-1996 + -- Notes: + -- a) Returns 0.0 on error condition + 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); + -- 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); + 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; + + -- Get value for special cases + if ( X = 0 and Y > 0.0 ) then + return 0.0; + end if; + + if ( X = 1 ) then + return 1.0; + end if; + + if ( Y = 0.0 and X /= 0 ) then + return 1.0; end if; - -- compute the result + if ( Y = 1.0) then + return (REAL(X)); + end if; + + -- Get value for general case 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 + function "**" (X : in REAL; Y : in REAL) return REAL is + -- Description: + -- See function declaration in IEEE Std 1076.2-1996 + -- Notes: + -- a) Returns 0.0 on error condition + 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); + -- 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); + 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; + + -- Get value for special cases + if ( X = 0.0 and Y > 0.0 ) then + return 0.0; + end if; + + if ( X = 1.0 ) then + return 1.0; end if; - -- compute the result + if ( Y = 0.0 and X /= 0.0 ) then + return 1.0; + end if; + + if ( Y = 1.0) then + return (X); + end if; + + -- Get value for general case return EXP (Y * LOG (X)); end "**"; - function EXP (X : real ) return real is + function EXP (X : in REAL ) return REAL is + -- Description: + -- See function declaration in IEEE Std 1076.2-1996 + -- Notes: + -- a) This function computes the exponential using the following + -- series: + -- exp(x) = 1 + x + x**2/2! + x**3/3! + ... ; |x| < 1.0 + -- and reduces argument X to take advantage of exp(x+y) = + -- exp(x)*exp(y) + -- + -- b) This implementation limits X to be less than LOG(REAL'HIGH) + -- to avoid overflow. Returns REAL'HIGH when X reaches that + -- limit + -- + constant EPS : REAL := BASE_EPS*BASE_EPS*BASE_EPS;-- Precision criteria + + variable RECIPROCAL: BOOLEAN := X < 0.0;-- Check sign of argument + variable XLOCAL : REAL := ABS(X); -- Use positive value + variable OLDVAL: REAL ; + variable COUNT: INTEGER ; + variable NEWVAL: REAL ; + variable LAST_TERM: REAL ; + variable FACTOR : REAL := 1.0; + + begin + -- Compute value for special cases + if X = 0.0 then + return 1.0; + end if; + + if XLOCAL = 1.0 then + if RECIPROCAL then + return MATH_1_OVER_E; + else + return MATH_E; + end if; + end if; + + if XLOCAL = 2.0 then + if RECIPROCAL then + return 1.0/MATH_E_P2; + else + return MATH_E_P2; + end if; + end if; + + if XLOCAL = 10.0 then + if RECIPROCAL then + return 1.0/MATH_E_P10; + else + return MATH_E_P10; + end if; + end if; + + if XLOCAL > LOG(REAL'HIGH) then + if RECIPROCAL then + return 0.0; + else + assert FALSE + report "X > LOG(REAL'HIGH) in EXP(X)" + severity NOTE; + return REAL'HIGH; + end if; + end if; + + -- Reduce argument to ABS(X) < 1.0 + while XLOCAL > 10.0 loop + XLOCAL := XLOCAL - 10.0; + FACTOR := FACTOR*MATH_E_P10; + end loop; + + while XLOCAL > 1.0 loop + XLOCAL := XLOCAL - 1.0; + FACTOR := FACTOR*MATH_E; + end loop; + + -- Compute value for case 0 < XLOCAL < 1 + OLDVAL := 1.0; + LAST_TERM := XLOCAL; + NEWVAL:= OLDVAL + LAST_TERM; + COUNT := 2; + + -- Check for relative and absolute errors and max count + while ( ( (ABS((NEWVAL - OLDVAL)/NEWVAL) > EPS) OR + (ABS(NEWVAL - OLDVAL) > EPS) ) AND + (COUNT < MAX_COUNT ) ) loop + OLDVAL := NEWVAL; + LAST_TERM := LAST_TERM*(XLOCAL / (REAL(COUNT))); + NEWVAL := OLDVAL + LAST_TERM; + COUNT := COUNT + 1; + end loop; + + -- Compute final value using exp(x+y) = exp(x)*exp(y) + NEWVAL := NEWVAL*FACTOR; + + if RECIPROCAL then + NEWVAL := 1.0/NEWVAL; + end if; + + return NEWVAL; + end EXP; + + + -- + -- Auxiliary Functions to Compute LOG + -- + function ILOGB(X: in REAL) return INTEGER IS + -- Description: + -- Returns n such that -1 <= ABS(X)/2^n < 2 + -- Notes: + -- None + + variable N: INTEGER := 0; + variable Y: REAL := ABS(X); + begin - assert false severity failure; - end EXP; + if(Y = 1.0 or Y = 0.0) then + return 0; + end if; - function c_log (x : real ) return real; - attribute foreign of c_log : function is "VHPIDIRECT log"; + if( Y > 1.0) then + while Y >= 2.0 loop + Y := Y/2.0; + N := N+1; + end loop; + return N; + end if; + + -- O < Y < 1 + while Y < 1.0 loop + Y := Y*2.0; + N := N -1; + end loop; + return N; + end ILOGB; - function c_log (x : real ) return real is + function LDEXP(X: in REAL; N: in INTEGER) RETURN REAL IS + -- Description: + -- Returns X*2^n + -- Notes: + -- None begin - assert false severity failure; - end c_log; + return X*(2.0 ** N); + end LDEXP; - function LOG (X : real ) return real is - -- returns natural logarithm of X; X > 0 + function LOG (X : in REAL ) return REAL IS + -- Description: + -- See function declaration in IEEE Std 1076.2-1996 -- - -- 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 - -- - 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; - return c_log(x); - end LOG; + -- Notes: + -- a) Returns REAL'LOW on error + -- + -- Copyright (c) 1992 Regents of the University of California. + -- All rights reserved. + -- + -- Redistribution and use in source and binary forms, with or without + -- modification, are permitted provided that the following conditions + -- are met: + -- 1. Redistributions of source code must retain the above copyright + -- notice, this list of conditions and the following disclaimer. + -- 2. Redistributions in binary form must reproduce the above copyright + -- notice, this list of conditions and the following disclaimer in the + -- documentation and/or other materials provided with the distribution. + -- 3. All advertising materials mentioning features or use of this + -- software must display the following acknowledgement: + -- This product includes software developed by the University of + -- California, Berkeley and its contributors. + -- 4. Neither the name of the University nor the names of its + -- contributors may be used to endorse or promote products derived + -- from this software without specific prior written permission. + -- + -- THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' + -- AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, + -- THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A + -- PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR + -- CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, + -- EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, + -- PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR + -- PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY + -- OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT + -- (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE + -- USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH + -- DAMAGE. + -- + -- NOTE: This VHDL version was generated using the C version of the + -- original function by the IEEE VHDL Mathematical Package + -- Working Group (CS/JT) - function LOG (X : in real; BASE: in real) return real is - -- returns logarithm base BASE of X; X > 0 - begin - -- check validity of argument - if ( BASE <= 0.0 ) or ( x <= 0.0 ) then - assert false report "BASE <= 0.0 or X <= 0.0 in LOG(BASE, X)" - severity ERROR; - return(REAL'LOW); - end if; - -- compute the value - return (LOG(X)/LOG(BASE)); + constant N: INTEGER := 128; + + -- Table of log(Fj) = logF_head[j] + logF_tail[j], for Fj = 1+j/128. + -- Used for generation of extend precision logarithms. + -- The constant 35184372088832 is 2^45, so the divide is exact. + -- It ensures correct reading of logF_head, even for inaccurate + -- decimal-to-binary conversion routines. (Everybody gets the + -- right answer for INTEGERs less than 2^53.) + -- Values for LOG(F) were generated using error < 10^-57 absolute + -- with the bc -l package. + + type REAL_VECTOR is array (NATURAL range <>) of REAL; + + constant A1:REAL := 0.08333333333333178827; + constant A2:REAL := 0.01250000000377174923; + constant A3:REAL := 0.002232139987919447809; + constant A4:REAL := 0.0004348877777076145742; + + constant LOGF_HEAD: REAL_VECTOR(0 TO N) := ( + 0.0, + 0.007782140442060381246, + 0.015504186535963526694, + 0.023167059281547608406, + 0.030771658666765233647, + 0.038318864302141264488, + 0.045809536031242714670, + 0.053244514518837604555, + 0.060624621816486978786, + 0.067950661908525944454, + 0.075223421237524235039, + 0.082443669210988446138, + 0.089612158689760690322, + 0.096729626458454731618, + 0.103796793681567578460, + 0.110814366340264314203, + 0.117783035656430001836, + 0.124703478501032805070, + 0.131576357788617315236, + 0.138402322859292326029, + 0.145182009844575077295, + 0.151916042025732167530, + 0.158605030176659056451, + 0.165249572895390883786, + 0.171850256926518341060, + 0.178407657472689606947, + 0.184922338493834104156, + 0.191394852999565046047, + 0.197825743329758552135, + 0.204215541428766300668, + 0.210564769107350002741, + 0.216873938300523150246, + 0.223143551314024080056, + 0.229374101064877322642, + 0.235566071312860003672, + 0.241719936886966024758, + 0.247836163904594286577, + 0.253915209980732470285, + 0.259957524436686071567, + 0.265963548496984003577, + 0.271933715484010463114, + 0.277868451003087102435, + 0.283768173130738432519, + 0.289633292582948342896, + 0.295464212893421063199, + 0.301261330578199704177, + 0.307025035294827830512, + 0.312755710004239517729, + 0.318453731118097493890, + 0.324119468654316733591, + 0.329753286372579168528, + 0.335355541920762334484, + 0.340926586970454081892, + 0.346466767346100823488, + 0.351976423156884266063, + 0.357455888922231679316, + 0.362905493689140712376, + 0.368325561158599157352, + 0.373716409793814818840, + 0.379078352934811846353, + 0.384411698910298582632, + 0.389716751140440464951, + 0.394993808240542421117, + 0.400243164127459749579, + 0.405465108107819105498, + 0.410659924985338875558, + 0.415827895143593195825, + 0.420969294644237379543, + 0.426084395310681429691, + 0.431173464818130014464, + 0.436236766774527495726, + 0.441274560805140936281, + 0.446287102628048160113, + 0.451274644139630254358, + 0.456237433481874177232, + 0.461175715122408291790, + 0.466089729924533457960, + 0.470979715219073113985, + 0.475845904869856894947, + 0.480688529345570714212, + 0.485507815781602403149, + 0.490303988045525329653, + 0.495077266798034543171, + 0.499827869556611403822, + 0.504556010751912253908, + 0.509261901790523552335, + 0.513945751101346104405, + 0.518607764208354637958, + 0.523248143765158602036, + 0.527867089620485785417, + 0.532464798869114019908, + 0.537041465897345915436, + 0.541597282432121573947, + 0.546132437597407260909, + 0.550647117952394182793, + 0.555141507540611200965, + 0.559615787935399566777, + 0.564070138285387656651, + 0.568504735352689749561, + 0.572919753562018740922, + 0.577315365035246941260, + 0.581691739635061821900, + 0.586049045003164792433, + 0.590387446602107957005, + 0.594707107746216934174, + 0.599008189645246602594, + 0.603290851438941899687, + 0.607555250224322662688, + 0.611801541106615331955, + 0.616029877215623855590, + 0.620240409751204424537, + 0.624433288012369303032, + 0.628608659422752680256, + 0.632766669570628437213, + 0.636907462236194987781, + 0.641031179420679109171, + 0.645137961373620782978, + 0.649227946625615004450, + 0.653301272011958644725, + 0.657358072709030238911, + 0.661398482245203922502, + 0.665422632544505177065, + 0.669430653942981734871, + 0.673422675212350441142, + 0.677398823590920073911, + 0.681359224807238206267, + 0.685304003098281100392, + 0.689233281238557538017, + 0.693147180560117703862); + + constant LOGF_TAIL: REAL_VECTOR(0 TO N) := ( + 0.0, + -0.00000000000000543229938420049, + 0.00000000000000172745674997061, + -0.00000000000001323017818229233, + -0.00000000000001154527628289872, + -0.00000000000000466529469958300, + 0.00000000000005148849572685810, + -0.00000000000002532168943117445, + -0.00000000000005213620639136504, + -0.00000000000001819506003016881, + 0.00000000000006329065958724544, + 0.00000000000008614512936087814, + -0.00000000000007355770219435028, + 0.00000000000009638067658552277, + 0.00000000000007598636597194141, + 0.00000000000002579999128306990, + -0.00000000000004654729747598444, + -0.00000000000007556920687451336, + 0.00000000000010195735223708472, + -0.00000000000017319034406422306, + -0.00000000000007718001336828098, + 0.00000000000010980754099855238, + -0.00000000000002047235780046195, + -0.00000000000008372091099235912, + 0.00000000000014088127937111135, + 0.00000000000012869017157588257, + 0.00000000000017788850778198106, + 0.00000000000006440856150696891, + 0.00000000000016132822667240822, + -0.00000000000007540916511956188, + -0.00000000000000036507188831790, + 0.00000000000009120937249914984, + 0.00000000000018567570959796010, + -0.00000000000003149265065191483, + -0.00000000000009309459495196889, + 0.00000000000017914338601329117, + -0.00000000000001302979717330866, + 0.00000000000023097385217586939, + 0.00000000000023999540484211737, + 0.00000000000015393776174455408, + -0.00000000000036870428315837678, + 0.00000000000036920375082080089, + -0.00000000000009383417223663699, + 0.00000000000009433398189512690, + 0.00000000000041481318704258568, + -0.00000000000003792316480209314, + 0.00000000000008403156304792424, + -0.00000000000034262934348285429, + 0.00000000000043712191957429145, + -0.00000000000010475750058776541, + -0.00000000000011118671389559323, + 0.00000000000037549577257259853, + 0.00000000000013912841212197565, + 0.00000000000010775743037572640, + 0.00000000000029391859187648000, + -0.00000000000042790509060060774, + 0.00000000000022774076114039555, + 0.00000000000010849569622967912, + -0.00000000000023073801945705758, + 0.00000000000015761203773969435, + 0.00000000000003345710269544082, + -0.00000000000041525158063436123, + 0.00000000000032655698896907146, + -0.00000000000044704265010452446, + 0.00000000000034527647952039772, + -0.00000000000007048962392109746, + 0.00000000000011776978751369214, + -0.00000000000010774341461609578, + 0.00000000000021863343293215910, + 0.00000000000024132639491333131, + 0.00000000000039057462209830700, + -0.00000000000026570679203560751, + 0.00000000000037135141919592021, + -0.00000000000017166921336082431, + -0.00000000000028658285157914353, + -0.00000000000023812542263446809, + 0.00000000000006576659768580062, + -0.00000000000028210143846181267, + 0.00000000000010701931762114254, + 0.00000000000018119346366441110, + 0.00000000000009840465278232627, + -0.00000000000033149150282752542, + -0.00000000000018302857356041668, + -0.00000000000016207400156744949, + 0.00000000000048303314949553201, + -0.00000000000071560553172382115, + 0.00000000000088821239518571855, + -0.00000000000030900580513238244, + -0.00000000000061076551972851496, + 0.00000000000035659969663347830, + 0.00000000000035782396591276383, + -0.00000000000046226087001544578, + 0.00000000000062279762917225156, + 0.00000000000072838947272065741, + 0.00000000000026809646615211673, + -0.00000000000010960825046059278, + 0.00000000000002311949383800537, + -0.00000000000058469058005299247, + -0.00000000000002103748251144494, + -0.00000000000023323182945587408, + -0.00000000000042333694288141916, + -0.00000000000043933937969737844, + 0.00000000000041341647073835565, + 0.00000000000006841763641591466, + 0.00000000000047585534004430641, + 0.00000000000083679678674757695, + -0.00000000000085763734646658640, + 0.00000000000021913281229340092, + -0.00000000000062242842536431148, + -0.00000000000010983594325438430, + 0.00000000000065310431377633651, + -0.00000000000047580199021710769, + -0.00000000000037854251265457040, + 0.00000000000040939233218678664, + 0.00000000000087424383914858291, + 0.00000000000025218188456842882, + -0.00000000000003608131360422557, + -0.00000000000050518555924280902, + 0.00000000000078699403323355317, + -0.00000000000067020876961949060, + 0.00000000000016108575753932458, + 0.00000000000058527188436251509, + -0.00000000000035246757297904791, + -0.00000000000018372084495629058, + 0.00000000000088606689813494916, + 0.00000000000066486268071468700, + 0.00000000000063831615170646519, + 0.00000000000025144230728376072, + -0.00000000000017239444525614834); + + variable M, J:INTEGER; + variable F1, F2, G, Q, U, U2, V: REAL; + variable ZERO: REAL := 0.0;--Made variable so no constant folding occurs + variable ONE: REAL := 1.0; --Made variable so no constant folding occurs + + -- double logb(), ldexp(); + + variable U1:REAL; + + begin + + -- Check validity of argument + if ( X <= 0.0 ) then + assert FALSE + report "X <= 0.0 in LOG(X)" + severity ERROR; + return(REAL'LOW); + end if; + + -- Compute value for special cases + if ( X = 1.0 ) then + return 0.0; + end if; + + if ( X = MATH_E ) then + return 1.0; + end if; + + -- Argument reduction: 1 <= g < 2; x/2^m = g; + -- y = F*(1 + f/F) for |f| <= 2^-8 + + M := ILOGB(X); + G := LDEXP(X, -M); + J := INTEGER(REAL(N)*(G-1.0)); -- C code adds 0.5 for rounding + F1 := (1.0/REAL(N)) * REAL(J) + 1.0; --F1*128 is an INTEGER in [128,512] + F2 := G - F1; + + -- Approximate expansion for log(1+f2/F1) ~= u + q + G := 1.0/(2.0*F1+F2); + U := 2.0*F2*G; + V := U*U; + Q := U*V*(A1 + V*(A2 + V*(A3 + V*A4))); + + -- Case 1: u1 = u rounded to 2^-43 absolute. Since u < 2^-8, + -- u1 has at most 35 bits, and F1*u1 is exact, as F1 has < 8 bits. + -- It also adds exactly to |m*log2_hi + log_F_head[j] | < 750. + -- + if ( J /= 0 or M /= 0) then + U1 := U + 513.0; + U1 := U1 - 513.0; + + -- Case 2: |1-x| < 1/256. The m- and j- dependent terms are zero + -- u1 = u to 24 bits. + -- + else + U1 := U; + --TRUNC(U1); --In c this is u1 = (double) (float) (u1) + end if; + + U2 := (2.0*(F2 - F1*U1) - U1*F2) * G; + -- u1 + u2 = 2f/(2F+f) to extra precision. + + -- log(x) = log(2^m*F1*(1+f2/F1)) = + -- (m*log2_hi+LOGF_HEAD(j)+u1) + (m*log2_lo+LOGF_TAIL(j)+q); + -- (exact) + (tiny) + + U1 := U1 + REAL(M)*LOGF_HEAD(N) + LOGF_HEAD(J); -- Exact + U2 := (U2 + LOGF_TAIL(J)) + Q; -- Tiny + U2 := U2 + LOGF_TAIL(N)*REAL(M); + return (U1 + U2); end LOG; - function LOG2 (X : in real) return real is - -- returns logarithm BASE 2 of X; X > 0 + + function LOG2 (X: in REAL) return REAL is + -- Description: + -- See function declaration in IEEE Std 1076.2-1996 + -- Notes: + -- a) Returns REAL'LOW on error begin - return LOG(X) / MATH_LOG_OF_2; + -- Check validity of arguments + if ( X <= 0.0 ) then + assert FALSE + report "X <= 0.0 in LOG2(X)" + severity ERROR; + return(REAL'LOW); + end if; + + -- Compute value for special cases + if ( X = 1.0 ) then + return 0.0; + end if; + + if ( X = 2.0 ) then + return 1.0; + end if; + + -- Compute value for general case + return ( MATH_LOG2_OF_E*LOG(X) ); end LOG2; - function LOG10 (X : in real) return real is - -- returns logarithm BASE 10 of X; X > 0 + + function LOG10 (X: in REAL) return REAL is + -- Description: + -- See function declaration in IEEE Std 1076.2-1996 + -- Notes: + -- a) Returns REAL'LOW on error begin - return LOG(X) / MATH_LOG_OF_10; + -- Check validity of arguments + if ( X <= 0.0 ) then + assert FALSE + report "X <= 0.0 in LOG10(X)" + severity ERROR; + return(REAL'LOW); + end if; + + -- Compute value for special cases + if ( X = 1.0 ) then + return 0.0; + end if; + + if ( X = 10.0 ) then + return 1.0; + end if; + + -- Compute value for general case + return ( MATH_LOG10_OF_E*LOG(X) ); end LOG10; - - function SIN (X : real ) return real is - begin - assert false severity failure; - end SIN; - - function COS (x : REAL) return REAL is - begin - assert false severity failure; - end COS; - - function TAN (x : REAL) return REAL is + + function LOG (X: in REAL; BASE: in REAL) return REAL is + -- Description: + -- See function declaration in IEEE Std 1076.2-1996 + -- Notes: + -- a) Returns REAL'LOW on error begin - assert false severity failure; - end TAN; - - function c_asin (x : real ) return real; - attribute foreign of c_asin : function is "VHPIDIRECT asin"; + -- Check validity of arguments + if ( X <= 0.0 ) then + assert FALSE + report "X <= 0.0 in LOG(X, BASE)" + severity ERROR; + return(REAL'LOW); + end if; + + if ( BASE <= 0.0 or BASE = 1.0 ) then + assert FALSE + report "BASE <= 0.0 or BASE = 1.0 in LOG(X, BASE)" + severity ERROR; + return(REAL'LOW); + end if; + + -- Compute value for special cases + if ( X = 1.0 ) then + return 0.0; + end if; + + if ( X = BASE ) then + return 1.0; + end if; + + -- Compute value for general case + return ( LOG(X)/LOG(BASE)); + end LOG; + + + function SIN (X : in REAL ) return REAL is + -- Description: + -- See function declaration in IEEE Std 1076.2-1996 + -- Notes: + -- a) SIN(-X) = -SIN(X) + -- b) SIN(X) = X if ABS(X) < EPS + -- c) SIN(X) = X - X**3/3! if EPS < ABS(X) < BASE_EPS + -- d) SIN(MATH_PI_OVER_2 - X) = COS(X) + -- e) COS(X) = 1.0 - 0.5*X**2 if ABS(X) < EPS + -- f) COS(X) = 1.0 - 0.5*X**2 + (X**4)/4! if + -- EPS< ABS(X) <BASE_EPS + + constant EPS : REAL := BASE_EPS*BASE_EPS; -- Convergence criteria + + variable N : INTEGER; + variable NEGATIVE : BOOLEAN := X < 0.0; + variable XLOCAL : REAL := ABS(X) ; + variable VALUE: REAL; + variable TEMP : REAL; - function c_asin (x : real ) return real is begin - assert false severity failure; - end c_asin; - - 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; + -- Make XLOCAL < MATH_2_PI + if XLOCAL > MATH_2_PI then + TEMP := FLOOR(XLOCAL/MATH_2_PI); + XLOCAL := XLOCAL - TEMP*MATH_2_PI; + end if; + + if XLOCAL < 0.0 then + assert FALSE + report "XLOCAL <= 0.0 after reduction in SIN(X)" + severity ERROR; + XLOCAL := -XLOCAL; + end if; + + -- Compute value for special cases + if XLOCAL = 0.0 or XLOCAL = MATH_2_PI or XLOCAL = MATH_PI then + return 0.0; + end if; + + if XLOCAL = MATH_PI_OVER_2 then + if NEGATIVE then + return -1.0; + else + return 1.0; + end if; + end if; + + if XLOCAL = MATH_3_PI_OVER_2 then + if NEGATIVE then + return 1.0; + else + return -1.0; + end if; + end if; + + if XLOCAL < EPS then + if NEGATIVE then + return -XLOCAL; + else + return XLOCAL; + end if; + else + if XLOCAL < BASE_EPS then + TEMP := XLOCAL - (XLOCAL*XLOCAL*XLOCAL)/6.0; + if NEGATIVE then + return -TEMP; + else + return TEMP; + end if; + end if; + end if; + + TEMP := MATH_PI - XLOCAL; + if ABS(TEMP) < EPS then + if NEGATIVE then + return -TEMP; + else + return TEMP; + end if; + else + if ABS(TEMP) < BASE_EPS then + TEMP := TEMP - (TEMP*TEMP*TEMP)/6.0; + if NEGATIVE then + return -TEMP; + else + return TEMP; + end if; + end if; + end if; + + TEMP := MATH_2_PI - XLOCAL; + if ABS(TEMP) < EPS then + if NEGATIVE then + return TEMP; + else + return -TEMP; + end if; + else + if ABS(TEMP) < BASE_EPS then + TEMP := TEMP - (TEMP*TEMP*TEMP)/6.0; + if NEGATIVE then + return TEMP; + else + return -TEMP; + end if; + end if; + end if; + + TEMP := ABS(MATH_PI_OVER_2 - XLOCAL); + if TEMP < EPS then + TEMP := 1.0 - TEMP*TEMP*0.5; + if NEGATIVE then + return -TEMP; + else + return TEMP; + end if; + else + if TEMP < BASE_EPS then + TEMP := 1.0 -TEMP*TEMP*0.5 + TEMP*TEMP*TEMP*TEMP/24.0; + if NEGATIVE then + return -TEMP; + else + return TEMP; + end if; + end if; + end if; + + TEMP := ABS(MATH_3_PI_OVER_2 - XLOCAL); + if TEMP < EPS then + TEMP := 1.0 - TEMP*TEMP*0.5; + if NEGATIVE then + return TEMP; + else + return -TEMP; + end if; + else + if TEMP < BASE_EPS then + TEMP := 1.0 -TEMP*TEMP*0.5 + TEMP*TEMP*TEMP*TEMP/24.0; + if NEGATIVE then + return TEMP; + else + return -TEMP; + end if; + end if; + end if; + + -- Compute value for general cases + if ((XLOCAL < MATH_PI_OVER_2 ) and (XLOCAL > 0.0)) then + VALUE:= CORDIC( KC, 0.0, x, 27, ROTATION)(1); + end if; + + N := INTEGER ( FLOOR(XLOCAL/MATH_PI_OVER_2)); + case QUADRANT( N mod 4) is + when 0 => + VALUE := CORDIC( KC, 0.0, XLOCAL, 27, ROTATION)(1); + when 1 => + VALUE := CORDIC( KC, 0.0, XLOCAL - MATH_PI_OVER_2, 27, + ROTATION)(0); + when 2 => + VALUE := -CORDIC( KC, 0.0, XLOCAL - MATH_PI, 27, ROTATION)(1); + when 3 => + VALUE := -CORDIC( KC, 0.0, XLOCAL - MATH_3_PI_OVER_2, 27, + ROTATION)(0); + end case; + + if NEGATIVE then + return -VALUE; else - return c_asin(x); - end if; - end ASIN; - - function c_acos (x : real ) return real; - attribute foreign of c_acos : function is "VHPIDIRECT acos"; - - function c_acos (x : real ) return real is + return VALUE; + end if; + end SIN; + + + function COS (X : in REAL) return REAL is + -- Description: + -- See function declaration in IEEE Std 1076.2-1996 + -- Notes: + -- a) COS(-X) = COS(X) + -- b) COS(X) = SIN(MATH_PI_OVER_2 - X) + -- c) COS(MATH_PI + X) = -COS(X) + -- d) COS(X) = 1.0 - X*X/2.0 if ABS(X) < EPS + -- e) COS(X) = 1.0 - 0.5*X**2 + (X**4)/4! if + -- EPS< ABS(X) <BASE_EPS + -- + constant EPS : REAL := BASE_EPS*BASE_EPS; + + variable XLOCAL : REAL := ABS(X); + variable VALUE: REAL; + variable TEMP : REAL; + begin - assert false severity failure; - end c_acos; - - 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; + -- Make XLOCAL < MATH_2_PI + if XLOCAL > MATH_2_PI then + TEMP := FLOOR(XLOCAL/MATH_2_PI); + XLOCAL := XLOCAL - TEMP*MATH_2_PI; + end if; + + if XLOCAL < 0.0 then + assert FALSE + report "XLOCAL <= 0.0 after reduction in COS(X)" + severity ERROR; + XLOCAL := -XLOCAL; + end if; + + -- Compute value for special cases + if XLOCAL = 0.0 or XLOCAL = MATH_2_PI then + return 1.0; + end if; + + if XLOCAL = MATH_PI then + return -1.0; + end if; + + if XLOCAL = MATH_PI_OVER_2 or XLOCAL = MATH_3_PI_OVER_2 then + return 0.0; + end if; + + TEMP := ABS(XLOCAL); + if ( TEMP < EPS) then + return (1.0 - 0.5*TEMP*TEMP); + else + if (TEMP < BASE_EPS) then + return (1.0 -0.5*TEMP*TEMP + TEMP*TEMP*TEMP*TEMP/24.0); + end if; + end if; + + TEMP := ABS(XLOCAL -MATH_2_PI); + if ( TEMP < EPS) then + return (1.0 - 0.5*TEMP*TEMP); + else + if (TEMP < BASE_EPS) then + return (1.0 -0.5*TEMP*TEMP + TEMP*TEMP*TEMP*TEMP/24.0); + end if; + end if; + + TEMP := ABS (XLOCAL - MATH_PI); + if TEMP < EPS then + return (-1.0 + 0.5*TEMP*TEMP); + else + if (TEMP < BASE_EPS) then + return (-1.0 +0.5*TEMP*TEMP - TEMP*TEMP*TEMP*TEMP/24.0); + end if; + end if; + + -- Compute value for general cases + return SIN(MATH_PI_OVER_2 - XLOCAL); + end COS; + + function TAN (X : in REAL) return REAL is + -- Description: + -- See function declaration in IEEE Std 1076.2-1996 + -- Notes: + -- a) TAN(0.0) = 0.0 + -- b) TAN(-X) = -TAN(X) + -- c) Returns REAL'LOW on error if X < 0.0 + -- d) Returns REAL'HIGH on error if X > 0.0 + + variable NEGATIVE : BOOLEAN := X < 0.0; + variable XLOCAL : REAL := ABS(X) ; + variable VALUE: REAL; + variable TEMP : REAL; + + begin + -- Make 0.0 <= XLOCAL <= MATH_2_PI + if XLOCAL > MATH_2_PI then + TEMP := FLOOR(XLOCAL/MATH_2_PI); + XLOCAL := XLOCAL - TEMP*MATH_2_PI; + end if; + + if XLOCAL < 0.0 then + assert FALSE + report "XLOCAL <= 0.0 after reduction in TAN(X)" + severity ERROR; + XLOCAL := -XLOCAL; + end if; + + -- Check validity of argument + if XLOCAL = MATH_PI_OVER_2 then + assert FALSE + report "X is a multiple of MATH_PI_OVER_2 in TAN(X)" + severity ERROR; + if NEGATIVE then + return(REAL'LOW); + else + return(REAL'HIGH); + end if; + end if; + + if XLOCAL = MATH_3_PI_OVER_2 then + assert FALSE + report "X is a multiple of MATH_3_PI_OVER_2 in TAN(X)" + severity ERROR; + if NEGATIVE then + return(REAL'HIGH); + else + return(REAL'LOW); + end if; + end if; + + -- Compute value for special cases + if XLOCAL = 0.0 or XLOCAL = MATH_PI then + return 0.0; + end if; + + -- Compute value for general cases + VALUE := SIN(XLOCAL)/COS(XLOCAL); + if NEGATIVE then + return -VALUE; + else + return VALUE; + end if; + end TAN; + + function ARCSIN (X : in REAL ) return REAL is + -- Description: + -- See function declaration in IEEE Std 1076.2-1996 + -- Notes: + -- a) ARCSIN(-X) = -ARCSIN(X) + -- b) Returns X on error + + variable NEGATIVE : BOOLEAN := X < 0.0; + variable XLOCAL : REAL := ABS(X); + variable VALUE : REAL; + + begin + -- Check validity of arguments + if XLOCAL > 1.0 then + assert FALSE + report "ABS(X) > 1.0 in ARCSIN(X)" + severity ERROR; + return X; + end if; + + -- Compute value for special cases + if XLOCAL = 0.0 then + return 0.0; + elsif XLOCAL = 1.0 then + if NEGATIVE then + return -MATH_PI_OVER_2; + else + return MATH_PI_OVER_2; + end if; + end if; + + -- Compute value for general cases + if XLOCAL < 0.9 then + VALUE := ARCTAN(XLOCAL/(SQRT(1.0 - XLOCAL*XLOCAL))); else - return c_acos(x); + VALUE := MATH_PI_OVER_2 - ARCTAN(SQRT(1.0 - XLOCAL*XLOCAL)/XLOCAL); end if; - end ACOS; - - function ATAN (x : REAL) return REAL is - -- returns -PI/2 < atan X < PI/2 + + if NEGATIVE then + VALUE := -VALUE; + end if; + + return VALUE; + end ARCSIN; + + function ARCCOS (X : in REAL) return REAL is + -- Description: + -- See function declaration in IEEE Std 1076.2-1996 + -- Notes: + -- a) ARCCOS(-X) = MATH_PI - ARCCOS(X) + -- b) Returns X on error + + variable NEGATIVE : BOOLEAN := X < 0.0; + variable XLOCAL : REAL := ABS(X); + variable VALUE : REAL; + begin - assert false severity failure; - end ATAN; + -- Check validity of argument + if XLOCAL > 1.0 then + assert FALSE + report "ABS(X) > 1.0 in ARCCOS(X)" + severity ERROR; + return X; + end if; + + -- Compute value for special cases + if X = 1.0 then + return 0.0; + elsif X = 0.0 then + return MATH_PI_OVER_2; + elsif X = -1.0 then + return MATH_PI; + end if; + + -- Compute value for general cases + if XLOCAL > 0.9 then + VALUE := ARCTAN(SQRT(1.0 - XLOCAL*XLOCAL)/XLOCAL); + else + VALUE := MATH_PI_OVER_2 - ARCTAN(XLOCAL/SQRT(1.0 - XLOCAL*XLOCAL)); + end if; - function c_atan2 (x : real; y : real) return real; - attribute foreign of c_atan2 : function is "VHPIDIRECT atan2"; - function c_atan2 (x : real; y: real) return real is - begin - assert false severity failure; - end c_atan2; - - 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 and x = 0.0 then - assert false - report "atan2(0.0, 0.0) is undetermined, returned 0,0" - severity NOTE; - return 0.0; + if NEGATIVE then + VALUE := MATH_PI - VALUE; + end if; + + return VALUE; + end ARCCOS; + + + function ARCTAN (Y : in REAL) return REAL is + -- Description: + -- See function declaration in IEEE Std 1076.2-1996 + -- Notes: + -- a) ARCTAN(-Y) = -ARCTAN(Y) + -- b) ARCTAN(Y) = -ARCTAN(1.0/Y) + MATH_PI_OVER_2 for |Y| > 1.0 + -- c) ARCTAN(Y) = Y for |Y| < EPS + + constant EPS : REAL := BASE_EPS*BASE_EPS*BASE_EPS; + + variable NEGATIVE : BOOLEAN := Y < 0.0; + variable RECIPROCAL : BOOLEAN; + variable YLOCAL : REAL := ABS(Y); + variable VALUE : REAL; + + begin + -- Make argument |Y| <=1.0 + if YLOCAL > 1.0 then + YLOCAL := 1.0/YLOCAL; + RECIPROCAL := TRUE; + else + RECIPROCAL := FALSE; + end if; + + -- Compute value for special cases + if YLOCAL = 0.0 then + if RECIPROCAL then + if NEGATIVE then + return (-MATH_PI_OVER_2); + else + return (MATH_PI_OVER_2); + end if; + else + return 0.0; + end if; + end if; + + if YLOCAL < EPS then + if NEGATIVE then + if RECIPROCAL then + return (-MATH_PI_OVER_2 + YLOCAL); + else + return -YLOCAL; + end if; + else + if RECIPROCAL then + return (MATH_PI_OVER_2 - YLOCAL); + else + return YLOCAL; + end if; + end if; + end if; + + -- Compute value for general cases + VALUE := CORDIC( 1.0, YLOCAL, 0.0, 27, VECTORING )(2); + + if RECIPROCAL then + VALUE := MATH_PI_OVER_2 - VALUE; + end if; + + if NEGATIVE then + VALUE := -VALUE; + end if; + + return VALUE; + end ARCTAN; + + + function ARCTAN (Y : in REAL; X : in REAL) return REAL is + -- Description: + -- See function declaration in IEEE Std 1076.2-1996 + -- Notes: + -- a) Returns 0.0 on error + + variable YLOCAL : REAL; + variable VALUE : REAL; + begin + + -- Check validity of arguments + if (Y = 0.0 and X = 0.0 ) then + assert FALSE report + "ARCTAN(0.0, 0.0) is undetermined" + severity ERROR; + return 0.0; + end if; + + -- Compute value for special cases + if Y = 0.0 then + if X > 0.0 then + return 0.0; + else + return MATH_PI; + end if; + end if; + + if X = 0.0 then + if Y > 0.0 then + return MATH_PI_OVER_2; else - return c_atan2(x,y); - end if; - end ATAN2; + return -MATH_PI_OVER_2; + end if; + end if; + + + -- Compute value for general cases + YLOCAL := ABS(Y/X); + + VALUE := ARCTAN(YLOCAL); + + if X < 0.0 then + VALUE := MATH_PI - VALUE; + end if; + + if Y < 0.0 then + VALUE := -VALUE; + end if; + return VALUE; + end ARCTAN; + + + function SINH (X : in REAL) return REAL is + -- Description: + -- See function declaration in IEEE Std 1076.2-1996 + -- Notes: + -- a) Returns (EXP(X) - EXP(-X))/2.0 + -- b) SINH(-X) = SINH(X) + + variable NEGATIVE : BOOLEAN := X < 0.0; + variable XLOCAL : REAL := ABS(X); + variable TEMP : REAL; + variable VALUE : REAL; - function SINH (X : real) return real is - -- hyperbolic sine; returns (e**X - e**(-X))/2 begin - assert false severity failure; + -- Compute value for special cases + if XLOCAL = 0.0 then + return 0.0; + end if; + + -- Compute value for general cases + TEMP := EXP(XLOCAL); + VALUE := (TEMP - 1.0/TEMP)*0.5; + + if NEGATIVE then + VALUE := -VALUE; + end if; + + return VALUE; end SINH; - function COSH (X : real) return real is - -- hyperbolic cosine; returns (e**X + e**(-X))/2 + function COSH (X : in REAL) return REAL is + -- Description: + -- See function declaration in IEEE Std 1076.2-1996 + -- Notes: + -- a) Returns (EXP(X) + EXP(-X))/2.0 + -- b) COSH(-X) = COSH(X) + + variable XLOCAL : REAL := ABS(X); + variable TEMP : REAL; + variable VALUE : REAL; begin - assert false severity failure; + -- Compute value for special cases + if XLOCAL = 0.0 then + return 1.0; + end if; + + + -- Compute value for general cases + TEMP := EXP(XLOCAL); + VALUE := (TEMP + 1.0/TEMP)*0.5; + + return VALUE; end COSH; - function TANH (X : real) return real is - -- hyperbolic tangent; -- returns (e**X - e**(-X))/(e**X + e**(-X)) + function TANH (X : in REAL) return REAL is + -- Description: + -- See function declaration in IEEE Std 1076.2-1996 + -- Notes: + -- a) Returns (EXP(X) - EXP(-X))/(EXP(X) + EXP(-X)) + -- b) TANH(-X) = -TANH(X) + + variable NEGATIVE : BOOLEAN := X < 0.0; + variable XLOCAL : REAL := ABS(X); + variable TEMP : REAL; + variable VALUE : REAL; + begin - assert false severity failure; + -- Compute value for special cases + if XLOCAL = 0.0 then + return 0.0; + end if; + + -- Compute value for general cases + TEMP := EXP(XLOCAL); + VALUE := (TEMP - 1.0/TEMP)/(TEMP + 1.0/TEMP); + + if NEGATIVE then + return -VALUE; + else + return VALUE; + end if; end TANH; - - function ASINH (X : real) return real is - -- returns ln( X + sqrt( X**2 + 1)) - begin - assert false severity failure; - end ASINH; - function c_acosh (x : real ) return real; - attribute foreign of c_acosh : function is "VHPIDIRECT acosh"; + function ARCSINH (X : in REAL) return REAL is + -- Description: + -- See function declaration in IEEE Std 1076.2-1996 + -- Notes: + -- a) Returns LOG( X + SQRT( X*X + 1.0)) - function c_acosh (x : real ) return real is begin - assert false severity failure; - end c_acosh; + -- Compute value for special cases + if X = 0.0 then + return 0.0; + end if; + + -- Compute value for general cases + return ( LOG( X + SQRT( X*X + 1.0)) ); + end ARCSINH; + + + + function ARCCOSH (X : in REAL) return REAL is + -- Description: + -- See function declaration in IEEE Std 1076.2-1996 + -- Notes: + -- a) Returns LOG( X + SQRT( X*X - 1.0)); X >= 1.0 + -- b) Returns X on error - 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 c_acosh(x); - end ACOSH; - - function c_atanh (x : real ) return real; - attribute foreign of c_atanh : function is "VHPIDIRECT atanh"; - - function c_atanh (x : real ) return real is begin - assert false severity failure; - end c_atanh; + -- Check validity of arguments + if X < 1.0 then + assert FALSE + report "X < 1.0 in ARCCOSH(X)" + severity ERROR; + return X; + end if; + + -- Compute value for special cases + if X = 1.0 then + return 0.0; + end if; - function ATANH (X : real) return real is - -- returns (ln( (1 + X)/(1 - X)))/2 ; | X | < 1 + -- Compute value for general cases + return ( LOG( X + SQRT( X*X - 1.0))); + end ARCCOSH; + + function ARCTANH (X : in REAL) return REAL is + -- Description: + -- See function declaration in IEEE Std 1076.2-1996 + -- Notes: + -- a) Returns (LOG( (1.0 + X)/(1.0 - X)))/2.0 ; | X | < 1.0 + -- b) Returns X on error begin - if abs x < 1.0 then - assert false report "Out of range parameter passed to ATANH" - severity ERROR; - return x; - end if; - return c_atanh(x); - end ATANH; + -- Check validity of arguments + if ABS(X) >= 1.0 then + assert FALSE + report "ABS(X) >= 1.0 in ARCTANH(X)" + severity ERROR; + return X; + end if; + + -- Compute value for special cases + if X = 0.0 then + return 0.0; + end if; + + -- Compute value for general cases + return( 0.5*LOG( (1.0+X)/(1.0-X) ) ); + end ARCTANH; end MATH_REAL; |