Skip to the content.
import numpy as np
import copy

class rational(object):
    @staticmethod
    def gcd(p,q):
        if(q>p):
            return rational.gcd(q,p)
        if(q<0):
            return rational.gcd(p,-q)
        if(q==0):
            return p
        return rational.gcd(q,p%q)
    
    @staticmethod
    def lcm(p,q):
        return p*(q//rational.gcd(p,q))

    @staticmethod
    def continuedFraction(num,eps):
        p = 1
        q = 0
        a = [int(num)]
        num -= a[0]
        if(num < 0):
            num += 1
            a = [a[0] - 1]
        

        r = [a[0],1]

        while(abs(num) > eps):
            num = 1/num
            a = int(num) 
            rNext = [r[0]*a+p,r[1]*a+q]
            
            if(rNext[1]*eps > 1):
                break

            p = r[0]
            q = r[1]
            r = copy.copy(rNext)
            num -= a

        return r

    def __init__(self,p,q=1,eps=2**-8):
        if(isinstance(p,rational) and q==1):
            self.p = p.p
            self.q = p.q
            self.eps = p.eps
            return
        self.eps = eps
        if(isinstance(p,int)):
            g = rational.gcd(p,q)
            rat = [p//g, q//g]
        else:
            rat = rational.continuedFraction(p,eps=eps)
        self.p,self.q = rat
        
    def __float__(self):    
        return self.p/self.q

    def _cmp(self,r):
        r = rational(r,eps=self.eps)
        diff = self - r
        return diff.p

    def __gt__(self,r):
        return self._cmp(r) > 0  

    def __ge__(self,r):
        return self._cmp(r) >= 0  

    def __le__(self,r):
        return self._cmp(r) <= 0  

    def __lt__(self,r):
        return self._cmp(r) < 0 

    def __eq__(self,r):
        return self._cmp(r) == 0

    def __add__(self,r):
        r = rational(r,eps=self.eps)
        l = rational.lcm(self.q,r.q)
        p1 = self.p * (l//self.q)
        p2 = r.p * (l//r.q)
        p = p1+p2
        q = l
        return rational(p,q)

    def __mul__(self,r):
        r = rational(r,eps=self.eps)
        g1 = rational.gcd(self.p,r.q)
        g2 = rational.gcd(self.q,r.p)
        p1 = self.p // g1
        p2 = r.p // g2
        q1 = self.q // g2
        q2 = r.q // g1
        p = p1*p2
        q = q1*q2
        return rational(p,q)

    def __truediv__(self,r):
        r = rational(r,eps=self.eps)
        rat = rational(r.q,r.p)
        return self*rat

    def __neg__(self):
        return rational(-self.p,self.q)

    def __sub__(self,r):
        r = rational(r,eps=self.eps)
        return self + (-r)

    def __isub__(self,rat):
        return self - rat

    def __imul__(self,rat):
        return self * rat

    def __iadd__(self,rat):
        return self + rat

    def __idiv__(self,rat):
        return self / rat
   
    def __repr__(self):
        return ("%d/%d(%f)" %(self.p,self.q,float(self)))
    
    def __str__(self):
        return ("%d/%d(%f)" %(self.p,self.q,float(self)))

    @property
    def inv(self):
        return rational(self.q,self.p)

    
class rationalMatrix(object): 
    @staticmethod
    def toRational(array,eps):
        shape = array.shape
        if(len(shape)==0):
            return rational(array,eps=eps)
        rat = []
        for a in array:
            rat.append(rationalMatrix.toRational(a,eps))
        return np.array(rat)
        
    def __init__(self,matrix,eps=2**-8):
        if(isinstance(matrix,rationalMatrix)):
            self.matrix = matrix
            self.eps = eps
            return
        self.matrix = rationalMatrix.toRational(matrix,eps=eps)
        self.eps = eps

    def __repr__(self):
        return self.matrix.__repr__()
   
    def __str__(self):
        return self.matrix.__str__()

    def __add__(self,other):
        return self.matrix + other.matrix

    def __sub__(self,other):
        return self.matrix - other.matrix

    def __mul__(self,other):
        return rationalMatrix(np.dot(self.matrix,other.matrix))

    @staticmethod
    def identity(n):
       return rationalMatrix(np.eye(n))