posits.py

from decimal import Decimal 

class Posit():
    
    @staticmethod
    def _twosComp(x,n):
        return ((x^((1<<n)-1))+1)&((1<<n)-1)

    @staticmethod
    def _getuptoRegime(k, n):  
        if k < 0:
            if abs(k) >= n - 1: 
                return 0, 0
            else: 
                rem_bits = n+k-2 
                return rem_bits, 1<< rem_bits 
        if k >= 0:
            if abs(k) >= n - 2: 
                return 0, ((1 << (n-1))-1)
            rem_bits = n-3-k 
            return rem_bits, ((1 << (k+2)) - 2) << (rem_bits)
        
    @staticmethod
    def _extractPosit(p: 'Posit'):
        signbit = ((p.value >> (p.n-1)) & 1)

        value = p.value
        if(signbit == 1):
            value=Posit._twosComp(value, p.n)
        # print(f"value:{value:08b}")

        r1 = ((value >> p.n-2) & 1) #first regime bit
        k = 0 if r1 else -1
        increment = 1 if r1 else -1

        m = 3
        while((m <= p.n) and (r1 == ((value >> p.n-m) & 1))):
            k += increment
            m += 1

        exp_bits = min(p.es, p.n-m)
        if(exp_bits < 0):
            return signbit, k, 0, 0, 0
        fraclength = max(0, p.n-m-exp_bits)
        
        fracbits = value & ((1<<fraclength)-1)
        exp = (value>>fraclength)&((1<<exp_bits)-1)
        return signbit, k, exp, fracbits, fraclength

    
    @staticmethod
    def fromBits(binary, n, es):
        temp = Posit(0,n,es)
        temp.value = binary
        return temp
        
    def __init__(self, number: Decimal, n: int, es: int):
        self.n = n
        self.es = es
        self.inf = 1<<(n-1)

        if number == 0: 
            self.value = 0
            return
        signbit = 1 if number < 0 else 0
        number=abs(number)
        exponent = 0
        if(number < 1):
            while(number < 1):
                exponent-=1
                number=number*2
        else:
            while(number >= 2):
                exponent+=1
                number=number/2
        number-=1
    
        e = exponent & ((1 << es) - 1)
        k = (exponent >> es)
        rem_bits, regime = Posit._getuptoRegime(k, n)
        # print(f"regime:{regime:08b}")
        if rem_bits <= es: 
            self.value = regime + (e>>(es-rem_bits))
        
        else: 
            frac_bits = rem_bits - es
            frac = int(round(max(number*2**frac_bits,0), 0)) #is max necessary?
            self.value = regime + (e << (frac_bits)) + frac

        if(signbit):
            # print(f"{self.value:08b}")
            self.value = Posit._twosComp(self.value, self.n)
            # print(f"{self.value:08b}")
            # self.value &= 1<<(self.n-1)

    def __neg__(self):
        # print(f"{self.value:0{self.n}b}")
        # value = (self.value)^(1<<(self.n-1))
        value = Posit._twosComp(self.value, self.n)
        # print(f"{(self.value^(1<<self.n)-1)+1:0{self.n}b}")
        return Posit.fromBits(value, self.n, self.es)
    

    def __add__(self, p2: 'Posit'):
        if(not isinstance(p2, Posit)):
            p2 = Posit(p2, self.n, self.es)
        if(self.value == self.inf or p2.value == self.inf):
            return Posit.fromBits(self.inf, self.n, self.es)
        if(self.value == 0):
            return p2
        if(p2.value == 0):
            return self
        signbit1, k1, exponent1, frac1, fraclength1 = Posit._extractPosit(self)
        signbit2, k2, exponent2, frac2, fraclength2 = Posit._extractPosit(p2)
        # print(signbit1, k1, exponent1, f"frac1:{frac1:0{fraclength1}b}")
        # print(signbit2, k2, exponent2, f"frac2:{frac2:0{fraclength2}b}")

        if(signbit1^signbit2):
            return self-(-p2)

        scale1 = (k1<<self.es)+exponent1
        scale2 = (k2<<self.es)+exponent2
        # print("scales:", scale1, scale2)

        # get fractions with max precision
        frac_bits = max(fraclength1, fraclength2)
        if(fraclength1 > fraclength2):
            frac2 <<= (fraclength1-fraclength2)
        else:
            frac1 <<= (fraclength2-fraclength1)

        # calculate offset before adding
        offset = (scale1 - scale2)
        if offset == 0:
            resultFrac = (frac1 + (1<<(frac_bits+1)) + frac2)
            resultScale = scale1
        elif offset > 0:
            resultFrac = (frac1 + (1<<frac_bits)) + ((frac2 + (1<<frac_bits)) >> offset)
            resultScale = scale1
        else: 
            resultFrac = (frac2 + (1<<frac_bits)) + ((frac1 + (1<<frac_bits)) >> -offset)
            resultScale = scale2
                
        # print(f"resultFrac: {resultFrac:0b}")
        if resultFrac.bit_length() > frac_bits+1:
            resultScale += 1

        resultExponent = resultScale%(1<<self.es)
        resultk= resultScale >> self.es

        # print(f"frac_bits: {frac_bits}")
        # print(f"resultFrac: {resultFrac:0b}")
        # print(f"resultscale: {resultScale}")

        rem_bits, regime = Posit._getuptoRegime(resultk, self.n)
        if rem_bits <= self.es: 
            return Posit.fromBits(regime + (resultExponent<<(self.es-rem_bits)),self.n, self.es)
        
        else: 
            frac_bits = rem_bits - self.es
            # while(resultFrac >= (1<<(frac_bits+1))):
            #     resultFrac//=2
            roundup = 0
            if(resultFrac.bit_length() > frac_bits+1):
                roundup = (resultFrac >> (resultFrac.bit_length()-frac_bits-2)) & 1
                resultFrac = resultFrac >> (resultFrac.bit_length()-frac_bits-1)
            
            resultFrac-=1<<frac_bits
            value = regime + (resultExponent << (frac_bits)) + resultFrac
            if(roundup and value != (1<<self.n)-1):
                value += 1
            if(signbit1):
                value = Posit._twosComp(value, self.n)
            return Posit.fromBits(value, self.n, self.es)
        
    def __radd__(self, other):
        return self+other
        
    def __sub__(self,p2):
        if(self.value == self.inf or p2.value == self.inf):
            return Posit.fromBits(self.inf, self.n, self.es)
        if(self.value == 0):
            return -p2
        if(p2.value == 0):
            return self 
        signbit1, k1, exponent1, frac1, fraclength1 = Posit._extractPosit(self)
        signbit2, k2, exponent2, frac2, fraclength2 = Posit._extractPosit(p2)
        # print(signbit1, k1, exponent1, f"frac1:{frac1:0{fraclength1}b}")
        # print(signbit2, k2, exponent2, f"frac2:{frac2:0{fraclength2}b}")

        if(signbit1^signbit2):
            return self+(-p2)
        
        if(k1==k2 and exponent1==exponent2 and frac1==frac2):
            return Posit(0,self.n, self.es)

        scale1 = (k1<<self.es)+exponent1
        scale2 = (k2<<self.es)+exponent2
        # print("scales:", scale1, scale2)

        # get fractions with max precision
        frac_bits = max(fraclength1, fraclength2)
        if(fraclength1 > fraclength2):
            frac2 <<= (fraclength1-fraclength2)
        else:
            frac1 <<= (fraclength2-fraclength1)
        # print(f"frac1:{frac1:0{frac_bits}b} frac2:{frac2:0{frac_bits}b}")

        # calculate offset before adding
        offset = (scale1 - scale2)
        if offset == 0:
            resultFrac = (frac1 - frac2)
            resultScale = scale1
        elif offset > 0:
            resultFrac = ((frac1 + (1<<frac_bits)) << offset) - (frac2 + (1<<frac_bits))
            resultScale = scale1
            # if((frac2>>(offset-1))&1):
            #     resultFrac-=1
        else: 
            signbit1 = signbit1^1
            resultFrac = ((frac2 + (1<<frac_bits)) << (-offset)) - (frac1 + (1<<frac_bits))
            resultScale = scale2
            # if((frac1>>(1-offset))&1):
            #     resultFrac-=1

        if(resultFrac < 0):
            signbit1 = signbit1^1
            resultFrac = abs(resultFrac)

        # bitlength = 1
        # while(resultFrac>>bitlength):
        #     bitlength+=1
        bitlength = resultFrac.bit_length() 

        # print(bitlength, frac_bits)

        # while (bitlength <= frac_bits):
        #     resultFrac = resultFrac<<1
        #     resultScale-=1
        #     bitlength+=1

        if(frac_bits > bitlength):
            resultFrac = resultFrac << (frac_bits-bitlength)
            resultScale-=(frac_bits-bitlength)
            bitlength+=(frac_bits-bitlength)

        resultScale += bitlength - 1 - abs(offset) - frac_bits

        # print(f"resultFrac: {resultFrac:0b}")
        # print(f"resultscale: {resultScale}")

        resultExponent = resultScale & resultScale%(1<<self.es)
        resultk= resultScale >> self.es

        rem_bits, regime = Posit._getuptoRegime(resultk, self.n)
        if rem_bits <= self.es: 
            return Posit.fromBits(regime + (resultExponent>>(self.es-rem_bits)),self.n,self.es)
        
        else: 
            frac_bits = rem_bits - self.es   
            roundup = 0         
            if(bitlength > frac_bits+1):
                # print(bin(resultFrac), frac_bits)
                roundup = (resultFrac >> (resultFrac.bit_length()-frac_bits-2)) & 1
                resultFrac = resultFrac>>(bitlength-frac_bits-1)
            if(bitlength < frac_bits+1):
                resultFrac = resultFrac<<(frac_bits+1-bitlength)
            resultFrac-=1<<frac_bits
            value = regime + (resultExponent << (frac_bits)) + resultFrac
            if(roundup and value != (1<<self.n)-1):
                value += 1
            if(signbit1):
                value = Posit._twosComp(value, self.n)
            return Posit.fromBits(value, self.n, self.es)
    
    def __mul__(self, p2):
        if(self.value == self.inf or p2.value == p2.inf):
            return Posit.fromBits(self.inf, self.n, self.es)
        if(self.value == 0 or p2.value == 0):
            return Posit.fromBits(0, self.n, self.es)
        sign1, k1, exponent1, frac1, fraclength1 = self._extractPosit(self) 
        sign2, k2, exponent2, frac2, fraclength2 = self._extractPosit(p2)
        
        # getting the result sign 
        signbit1 = sign1^sign2

        # calculating 
        scale1 = (k1<<self.es)+exponent1
        scale2 = (k2<<self.es)+exponent2
        resultScale = (scale1 + scale2)

        resultFrac = ((1<<fraclength1) + frac1) * ((1<<fraclength2) + frac2)    
        roundup = 0 # this is to "round up" the result 
        
        # Normalization
        while resultFrac >= 1<<(fraclength2+fraclength1 + 1):
            # roundup = resultFrac&1
            resultFrac>>= 1
            resultScale += 1
        
        # removing the extra 1
        resultFrac = ((resultFrac ) & ((1<<(fraclength2 + fraclength1)) - 1)) 
        
        # extracting k and e 
        resultk= resultScale>>self.es
        resultExponent = resultScale%(1<<self.es)

        # get sign + regime 
        rem_bits, sign_w_regime = self._getuptoRegime(resultk, self.n)

        # rebuild the posit 
        if rem_bits <= self.es: 
            value = sign_w_regime + (resultExponent>>(self.es-rem_bits))
        
        else: 
            frac_bits = rem_bits - self.es
            roundup = 0
            if(fraclength1+fraclength2 > frac_bits):
                roundup = (resultFrac >> (fraclength2 + fraclength1 - frac_bits - 1)) & 1
                frac = resultFrac >> (fraclength1 + fraclength2 - frac_bits)
            else:
                frac = resultFrac << (frac_bits - fraclength1 - fraclength2)
            # print(bin(frac), resultk, resultExponent)
            value = sign_w_regime + (resultExponent << (frac_bits)) + frac
            if(roundup and value != (1<<self.n)-1):
                value += 1
        
        if signbit1:
                value = Posit._twosComp(value, self.n)
        return Posit.fromBits(value, self.n, self.es)


    def __truediv__(self, p2):
        if(self.value == self.inf or p2.value == p2.inf):
            return Posit.fromBits(self.inf, self.n, self.es)
        if(self.value == 0):
            return Posit.fromBits(0, self.n, self.es)
        if(p2.value == 0):
            raise ZeroDivisionError
        sign1, k1, exponent1, frac1, fraclength1 = Posit._extractPosit(self)
        sign2, k2, exponent2, frac2, fraclength2 = Posit._extractPosit(p2)
        signbit = sign1 ^ sign2

        scale1 = (k1<<self.es)+exponent1
        scale2 = (k2<<self.es)+exponent2

        resultScale = scale1 - scale2

        if fraclength2 > fraclength1:
            resultFracLength = fraclength2
            frac1 <<= (fraclength2 - fraclength1)
        else:
            resultFracLength = fraclength1
            frac2 <<= (fraclength1 - fraclength2)

        numerator = frac1 + (1 << resultFracLength)
        denominator = frac2 + (1 << resultFracLength)
        resultFrac = (numerator << resultFracLength) // denominator

        if resultFrac < (1 << resultFracLength):
            resultScale -= 1
            resultFrac <<= 1

        resultFrac -= (1 << resultFracLength)
        resultExponent = resultScale%(1<<self.es)
        resultk = resultScale >> self.es
        rem_bits, sign_w_regime = Posit._getuptoRegime(resultk, self.n)

        if rem_bits <= self.es:
            value = sign_w_regime + (resultExponent >> (self.es - rem_bits))
        else:
            frac_bits = rem_bits - self.es
            if resultFracLength > frac_bits:
                roundup = (resultFrac >> (resultFracLength - frac_bits - 1)) & 1
                frac = resultFrac >> (resultFracLength - frac_bits)
                frac = frac + roundup
            else:
                frac = resultFrac << (frac_bits - resultFracLength)
            value = sign_w_regime + (resultExponent << frac_bits) + frac

        if signbit:
            value = Posit._twosComp(value, self.n)
        return Posit.fromBits(value, self.n, self.es)

    def __repr__(self):
        if self.value == 0:
            return "0"
        if self.value == self.inf: 
            return "inf"
        sign, k, exp, fracbits, fraclength = Posit._extractPosit(self)

        sign = -1 if sign else 1
        useed = 2 ** (2 ** self.es)
        value = sign * Decimal(useed ** k) * Decimal(2**exp) * (1+ fracbits/Decimal(2**fraclength))
        return str(value)
    
    def __float__(self):
        if self.value == 0:
            return 0.0
        if self.value == self.inf: 
            return float("inf")
        sign, k, exp, fracbits, fraclength = Posit._extractPosit(self)

        sign = -1 if sign else 1
        useed = 2 ** (2 ** self.es)
        value = sign * (useed ** k) * (2**exp) * (1+ fracbits/(2**fraclength))
        return value
        
    def toDecimal(self):
        if self.value == 0:
            return Decimal("0")
        if self.value == self.inf: 
            return Decimal('Infinity')
        sign, k, exp, fracbits, fraclength = Posit._extractPosit(self)

        sign = -1 if sign else 1
        useed = 2 ** (2 ** self.es)
        value = sign * Decimal(useed ** k) * 2**exp * (1+ Decimal(fracbits)/(2**Decimal(fraclength)))
        return value

    def toBytes(self):
        return self.value.to_bytes((self.n+7)//8)
    
    def is_nar(self):
        return self.value == (1 << (self.n - 1))  # same condition as in __repr__ for "inf"

    
    def sigmoid(self):
        return Posit.fromBits((((self.value ^ (1 << (self.n - 1))) >> 2)), self.n, self.es)