Comment to MathOverflow question

mo301110.py

# for python2/python3 portability
from __future__ import (print_function, division)

# parameters

limit = 10000000 # Check up to limit
trace = False    # Show trace of positives

# code

_checked = {0: (False, ()), 1: (True, (1,))}

def Check(n):
    if not isinstance(n,int) or n < 0:
        raise ValueError
    if not n in _checked:
        r2, r3 = n % 2, n % 3
        if n == 0:
            rv = (False, ())
        elif n == 1:
            rv = (True, ())
        elif r2 == 0 and r3 == 0:
            rv = Check(n // 2)
        elif r2 == 0 and r3 == 1:
            rv = Check(n // 2)
        elif r2 == 0 and r3 == 2:
            rv = Check(n // 3)
            if not rv[0]:
                 rv = Check(n // 2)
        elif r2 == 1 and r3 == 0:
            rv = Check(n // 6)
        elif r2 == 1 and r3 == 1:
            rv = (False, ())
        elif r2 == 1 and r3 == 2:
            rv = Check((n-2) // 3)
        else:
            raise ValueError
        _checked[n] = (rv[0], (n,) + rv[1])
    return _checked[n]


if __name__ == '__main__':
    
    count = 0
    for n in range(1,limit):
        value = Check(n)
        count += int(value[0])
        if trace and value[0]:
            print(n, ':', *value[1:])

    print("density = ", count / limit, '=', count, '/', limit)