Collatz sequence generator

collatz.py

#!/usr/bin/env python
# encoding: utf-8
from random import random
from heapq import heappush, nsmallest, nlargest
import time, sys
from sets import Set

def rotl32(a, b):
    return (((a << (b & 0x1f)) & 0xffffffff) |
            ((a >> (32 - (b & 0x1f))) & 0xffffffff))

def fmix(h):
    h = h ^ (h >> 16)
    h = (h * 0x85ebca6b) & 0xffffffff
    h = h ^ (h >> 13)
    h = (h * 0xc2b2ae35) & 0xffffffff
    h = h ^ (h >> 16)
    return h

def mmhash(value):
    import struct
    seed = 0
    length = len(value)
    num_blocks = length / 4
    tail_length = length % 4
    fmt = "<" + ("i" * num_blocks) + ("b" * tail_length)
    vals = struct.unpack(fmt, value)

    h1 = seed
    c1 = 0xcc9e2d51
    c2 = 0x1b873593
    for block in vals[:num_blocks]:
        k1 = block
        k1 = (k1 * c1) & 0xffffffff
        k1 = rotl32(k1, 15)
        k1 = (k1 * c2) & 0xffffffff

        h1 = h1 ^ k1
        h1 = rotl32(h1, 13)
        h1 = (h1 * 5 + 0xe6546b64) & 0xffffffff

    k1 = 0
    if tail_length >= 3:
        k1 = k1 ^ ((vals[num_blocks + 2] << 16) & 0xffffffff)
    if tail_length >= 2:
        k1 = k1 ^ ((vals[num_blocks + 1] << 8) & 0xffffffff)
    if tail_length >= 1:
        k1 = k1 ^ ( vals[num_blocks] & 0xffffffff)
        k1 = (k1 * c1) & 0xffffffff
        k1 = rotl32(k1, 15)
        k1 = (k1 * c2) & 0xffffffff
        h1 = h1 ^ k1

    h1 = h1 ^ (length & 0xffffffff)
    return fmix(h1)



def collatz(value):
    n = value
    yield n
    while n != 1:
        if n % 2 == 0: 
            n = n / 2
            yield n
        else:
            n = 3 * n +1
            yield n

def collatz_single(n):
    if n % 2 == 0: 
        n = n / 2
        return n
    else:
        n = 3 * n +1
        return n

def C(N):
    for n in range(1, N+1):
        for v in collatz(n):
            yield v

def h(x):
    return mmhash(str(x)) & 0xffffffff
    
def calc_with_dict(N):
    values = {}
    iterations = 0
    max_val = 0
    for x in C(N):
        iterations += 1
        if x > max_val: 
            max_val = x
        try:
            values[x] += 1  
        except KeyError:
            values[x] = 1

    carinality = len(values.keys())
    print 'N = {0}: maxCn {1}, |C_N| {2}, len(C_N) {3}'.format(N, max_val, carinality,iterations)

def quadratic_time(N):
    max_val = 1
    iterations = 0
    i = 0
    while i != N:
        i += 1
        iterations += 1
        n = i
        while n != 1:
            n = collatz_single(n)
            iterations += 1
            if n > max_val:
                max_val = n

    carinality = 1
    counter = 0
    while counter <= max_val:
        counter += 1
        i = 0
        inc_counter = False
        while i <= N and not inc_counter:
            i += 1
            n = i
            while n != 1 and not inc_counter:
                if n == counter:
                    carinality += 1
                    inc_counter = True
                n = collatz_single(n)
    print 'N = {0}: maxCn {1}, |C_N| {2}, len(C_N) {3}'.format(N, max_val, carinality,iterations)

def calc_with_limit_memory(N):
    t = 1000
    pheap = []
    size = 0
    iterations = 0
    s = Set()
    for x in C(N):
        v = h(x)
        
        if len(pheap) < t and v not in s: 
            s.add(v)
            heappush(pheap, (h(x), x))
        elif v < nlargest(1, pheap)[0][0] and v not in s:
            s.add(v)
            heappush(pheap, (h(x), x))            
            if size > t:
                s.remove(nlargest(1, pheap)[0][0])                
                pheap = nsmallest(t, pheap)
                
            else:
                size += 1
        iterations += 1

    largest_hash = nlargest(1, pheap)[0][0]
    print nlargest(1, pheap)
    print iterations
    
    R = 0xffffffff   
    print R, len(pheap), largest_hash
    nbr_distinct_elements = R * len(pheap) / largest_hash
    print nbr_distinct_elements

if __name__ == "__main__":
    """
    aa = [10,100,1000,10000,100000,1000000]
    for x in aa:
        calc_with_dict(x)
    """
    """
    old_time = time.time()
    quadratic_time(300)
    print(time.time()-old_time)
    exit()
    """
    calc_with_limit_memory(100000)