package FloatingPoint;





// ************************************************************************** //
//                                                                            //
//    eses                   eses                                             //
//   eses                     eses                                            //
//  eses    eseses  esesese    eses   Embedded Systems Group                  //
//  ese    ese  ese ese         ese                                           //
//  ese    eseseses eseseses    ese   Department of Computer Science          //
//  eses   eses          ese   eses                                           //
//   eses   eseses  eseseses  eses    University of Kaiserslautern            //
//    eses                   eses                                             //
//                                                                            //
// ************************************************************************** //
// This module implements multiplication of IEEE-like floating point numbers. //
// The algorithm below deals with the full IEEE standard including denormal   //
// and normal numbers, infinity and NaN values.                               //
// ************************************************************************** //

// parameters of the floating-point format
macro MantWidth = 4;
macro ExpWidth  = 3;
macro ExpMax    = exp(2,ExpWidth)-1;
macro ExpBias   = ExpMax / 2;
macro ExpMin    = +1 - ExpBias;
macro Width     = 1+ExpWidth+MantWidth;

macro MantMinNorm = exp(2,MantWidth-1);  // smallest normalized mantissa
macro MantMax     = exp(2,MantWidth)-1;  // largest mantissa value

// macros for extracting parts of the floating point number
macro Sign(x) = x{-1};
macro Exp(x)  = x{-2:MantWidth};
macro Mant(x) = x{MantWidth-1:0};
macro RmHBit(x) = Mant(x);            // remove hidden bit

// macros for interpretation of the floating point numbers
macro ValSign(s)    = (s?-1.0:1.0);
macro ValNMant(m)   = bv2nat(true@m)  * exp(2.0,-1*MantWidth);
macro ValDMant(m)   = bv2nat(false@m) * exp(2.0,-1*MantWidth);
macro ValMant(e,m)  = (e=={false::ExpWidth} ? ValDMant(m) : ValNMant(m));
macro ValNExp(e)    = +bv2nat(e)-ExpBias;
macro ValDExp(e)    = +1-ExpBias;
macro ValExp(e)     = (e=={false::ExpWidth} ? ValDExp(e) : ValNExp(e));
macro Float2Real(x) = ValSign(Sign(x)) * ValMant(Exp(x),Mant(x)) * exp(2.0,ValExp(Exp(x)));

// special values
macro posZero = false@{false::ExpWidth}@{false::MantWidth};
macro negZero =  true@{false::ExpWidth}@{false::MantWidth};
macro posInf  = false@{ true::ExpWidth}@{false::MantWidth};
macro negInf  =  true@{ true::ExpWidth}@{false::MantWidth};
macro posNaN  = false@{ true::ExpWidth}@{ true::MantWidth};
macro negNaN  =  true@{ true::ExpWidth}@{ true::MantWidth};


// Definition of some operations: declaration as macro
macro isZero(x)      = ( ( Exp(x) == {false::ExpWidth} ) & ( Mant(x) == {false::MantWidth} ) );
macro isNormal(x)    = ( 1 <= bv2nat(Exp(x)) ) & ( bv2nat(Exp(x)) <= (exp(2,ExpWidth)-2) );
macro isSubnormal(x) = ( Exp(x) == {false::ExpWidth} ) & ( Mant(x) != {false::MantWidth} );
macro isFinite(x)    = ( isZero(x) | isSubnormal(x) | isNormal(x) );
macro isInfinite(x)  = ( Exp(x) == {true::ExpWidth} ) & ( Mant(x) == {false::MantWidth} );
macro isNaN(x)       = ( Exp(x) == {true::ExpWidth} ) & ( Mant(x) != {false::MantWidth} );
macro isSignMinus(x) = ( Sign(x) != false );

// Macros that return a float
macro fAbs(x)         = (false@Exp(x)@Mant(x));
macro fCopySign(x,y)  = (Sign(y)@Exp(x)@Mant(x));
macro fNegate(x)      = ((Sign(x) ? false : true)@Exp(x)@Mant(x));


// here comes the main module: it first deals with special cases like NaN, Zero and Infinities
// according to the IEEE standard and then (in the largest part of case construct) handles
// normal and denormal numbers

module FloatMult(bv{Width} ?x,?y,!z,event !rdy) {
    bv{2*MantWidth+2} p;
    int e;
    bool round,sticky;
    bv{ExpWidth}  z_exp;
    bv{MantWidth} z_mnt;

    case // define special cases according to IEEE standard
        (isNaN(x) | isNaN(y) | isZero(x) & isInfinite(y) | isZero(y) & isInfinite(x)) do {
            z = (Sign(x) xor Sign(y) ? negNaN : posNaN);
            emit(rdy);
            }
        (isInfinite(x) | isInfinite(y)) do {
            z = (Sign(x) xor Sign(y) ? negInf : posInf);
            emit(rdy);
            }
        (isZero(x) | isZero(y)) do {
            z = (Sign(x) xor Sign(y) ? negZero : posZero);
            emit(rdy);
            }
    default {
        // x and y are (de)normal numbers different to zero

        // compute exact product by adding hidden bits; note that the radix
        // point is at now between positions 2*MantWidth and 2*MantWidth-1
        let(nx = isNormal(x) )
        let(ny = isNormal(y) )
        let(mx = bv2nat(nx@Mant(x)))
        let(my = bv2nat(ny@Mant(y)))
        p = nat2bv(mx*my, 2*MantWidth+2);

        // determine suitable exponent: we add 2 to shift the radix point
        // to the left of p{2*MantWidth+1}
        e = 2 + ExpBias + ValExp(Exp(x)) + ValExp(Exp(y));

        // if exponent should be less than 1, we shift the product to the right
        while(e<1) {
            next(p) = false@p{-1:1};
            next(e) = e+1;
            next(round)  = p{0};
            next(sticky) = sticky | round;
            pause;
        }
        
        // now we have e>=1, but p might have lost all 1s
        if(p=={false::2*MantWidth+2}) {
            z_mnt = {false::MantWidth};
            z_exp = nat2bv(abs(e),ExpWidth);
        } else {
            // If p still has 1s, we try to shift p to the left until it starts
            // with a 1. However, we do not want to destroy the property e>=1.
            while(p{-1}!=true & e>1){
                next(p) = p{-2:0}@false;
                next(e) = e-1;
                pause;
            }
    
            // now we have p{-1}&e>=1 or e==1
            if(e==1) {
                // product is a denormal number
                z_mnt  = p{2*MantWidth+1:MantWidth+2};
                round  = p{MantWidth+1};
                sticky = bv2nat(p{MantWidth:0})!=0;
                z_exp  = {false::ExpWidth};
            } else if(e>bv2nat({true::ExpWidth})) {
                // product is larger than largest representable number
                z_mnt  = {false::MantWidth};
                round  = false;
                sticky = false;
                z_exp  = {true::ExpWidth};
            } else {
                // product is a normal number: shift p to the left to
                // skip the hidden bit; thus decrease e once more
                z_mnt  = p{2*MantWidth:MantWidth+1};
                round  = p{MantWidth};
                sticky = bv2nat(p{MantWidth-1:0})!=0;
                z_exp  = nat2bv(abs(e-1),ExpWidth);
            }
        }

        // possible rounding of mantissa
        let(zf = (Sign(x) xor Sign(y))@z_exp@z_mnt)
        if(round & (sticky | z_mnt{0})) {
            // round towards infinity
            z = nat2bv(bv2nat(zf)+1,Width);
            // TODO: check overflow
        } else {
            // round towards zero
            z = zf;
        }
        emit(rdy);
    }
}
drivenby {
    // two normal numbers yield denormal result with round & sticky
    x = (false@0b001@0b1010);
    y = (false@0b010@0b0011);
    immediate abort { await(false); } when (rdy);
    assert(z==(false@0b000@0b1111));
}
drivenby {
    // two normal numbers yield denormal result with !round & sticky
    x = (false@0b001@0b1010);
    y = (false@0b010@0b0010);
    immediate abort { await(false); } when (rdy);
    assert(z==(false@0b000@0b1111));
}
drivenby {
    // two normal numbers yield denormal result with round & !sticky
    x = (false@0b001@0b0000);
    y = (false@0b010@0b0000);
    immediate abort { await(false); } when (rdy);
    assert(z==(false@0b000@0b1000));
}