Bitwise operations for fractions in Python.

fractions_bitwise.py

"""This module defines bitwise operations for fractions.

In the same way that rational numbers have repeating decimal expansions ie. 1/3 = 0.333...
they also have repeating binary expansions: 1/3 = 0.010101...
Bitwise operators can be extended to such expansions and then converted back to rational numbers.
Just like we usually don't write 1 as 0.999... this module follows the convention that 
1 = 1.000... even when expanded in binary.
The inversion operator '~' is not defined as ~1 = ...1110.111... would violate this convention.
For the same reason bitwise operations between negative numbers result in undefined behaviour.
The option to use one's complement for negative numbers ie. -1 = ...1110.111... = ~1
is not pursued here because fractions don't have a signed zero.
The inverse of zero ~0 = ...111.111... could not be converted back to a fraction.
"""
from fractions import *

__author__ = "Pyry Pakkanen"
__copyright__ = "Copyright 2011"
__credits__ = ["Pyry Pakkanen"]
__license__ = "MIT"
__version__ = "0.1"
__maintainer__ = "Pyry Pakkanen"
__email__ = "frostburn@suomi24.fi"
__status__ = "initial release"

def fraction_decompose(x):
    """Decompose a rational number to an "integer" part and a repeating binary expansion.

    x = 2^n ( a + b/(2^m - 1) )
    The returned coefficients n, a, b, m are integers that are as small 
    as possible in absolute value with b < 2^m - 1.
    a is the "integer" part of x
    b is the repeating fraction
    m is the period of the repetition
    n indicates the scale by moving the fraction point.
    """
    x = Fraction(x)
    p = x.numerator
    q = x.denominator
    n = 0
    #Get rid of factors of 2.
    while p%2 == 0:
        p //= 2
        n += 1
    while q%2 == 0:
        q //= 2
        n -= 1
    
    a = p//q
    p = p%q
    l = []
    b = []
    #Find the period and store the bits of the fractional part.
    while p not in l:
        l.append(p)
        b.append(2*p > q)
        p = (2*p)%q
        
    return n, a, sum(bit*2**i for i, bit in enumerate(b[::-1])), len(l)

def expand(x, y):
    """Expands x, y so that they have common decompositions.

    x = 2^n ( a_x + b_x/(2^m 1) ),
    y = 2^n ( a_y + b_y/(2^m-1) ), 
    with common n and m.

    returns n, a_x, b_x, a_y, b_y, m
    """
    nx, ax, bx, mx = fraction_decompose(x)
    ny, ay, by, my = fraction_decompose(y)

    nr = min(nx, ny)
    ax = ax*2**(nx-nr)
    ay = ay*2**(ny-nr)
    bx = bx*2**(nx-nr)
    by = by*2**(ny-nr)

    ax += bx//(2**mx-1)
    bx = bx%(2**mx-1)
    ay += by//(2**my-1)
    by = by%(2**my-1)

    mr = mx*my//gcd(mx, my)
    bx = sum(bx*2**(mx*i) for i in range(mr//mx))
    by = sum(by*2**(my*i) for i in range(mr//my))
    
    return nr, ax, bx, ay, by, mr

#Unfortunately the not operator cannot be defined unambiguously.
#def fraction_not(x):
#    nx, ax, bx, mx = fraction_decompose(x)
#    return Fraction(2)**nx*((~ax) + Fraction((~bx)&(2**mx-1),2**mx-1))

def fraction_and(x,y):
    nr, ax, bx, ay, by, mr = expand(x, y)
    return Fraction(2)**nr*((ax&ay) + Fraction(bx&by,2**mr-1))

def fraction_or(x,y):
    nr, ax, bx, ay, by, mr = expand(x, y)
    return Fraction(2)**nr*((ax|ay) + Fraction(bx|by,2**mr-1))

def fraction_xor(x,y):
    nr, ax, bx, ay, by, mr = expand(x, y)
    return Fraction(2)**nr*((ax^ay) + Fraction(bx^by,2**mr-1))

class BitFraction(Fraction):
    """Utility class such that &, |, ^ work as expected.

     The operators << and >> work reasonably with integer shifts but may
     produce unexpected results otherwise due to conversion to float and back.
     """
    def __and__(self,other):
        return BitFraction(fraction_and(self,other))
    def __rand__(self,other):
        return self&other
    def __or__(self,other):
        return BitFraction(fraction_or(self,other))
    def __ror__(self,other):
        return self|other
    def __xor__(self,other):
        return BitFraction(fraction_xor(self,other))
    def __rxor__(self,other):
        return self^other
    def __lshift__(self,other):
        return BitFraction(self*BitFraction(2)**other)
    def __rlshift__(self,other):
        return BitFraction(other*BitFraction(2)**self)
    def __rshift__(self,other):
        return BitFraction(self*BitFraction(1,2)**other)
    def __rrshift__(self,other):
        return BitFraction(other*BitFraction(1,2)**self)
    def __repr__(self):
        return ('BitFraction(%s, %s)' % (self._numerator, self._denominator))