Plot complex rational numbers grouped by the absolute value of the denominator in reduced form.

complex_rationals.py

from __future__ import division
from argparse import ArgumentParser
from itertools import cycle

from pylab import *
tau = 2 * pi


if __name__ == '__main__':
    parser = ArgumentParser(description='Plot complex rational numbers')
    parser.add_argument('N', default=4, nargs='?', type=int, metavar='Q', help='Maximum denominator size')
    parser.add_argument('-e', '--eisenstein', dest='eisenstein', action='store_true', help='Use Eisenstein integers')
    parser.add_argument('-s', '--size', dest='size', default=0.1, type=float, help='Size of the circles')
    parser.add_argument('-d', '--decay', dest='decay', default=1.1, type=float, help='Circle size decay constant')
    args = parser.parse_args()
    N = args.N
    size = args.size
    decay = args.decay
    ps = []
    qs = []
    w = exp(1j * tau / 3) if args.eisenstein else 1j
    for i in arange(-N, N+1):
        for j in arange(-N, N+1):
            ps.append(i + w * j)
            # Keep denominators in the first quadrant.
            if i > 0 and j > 0:
                qs.append(i + w * j)

    ps.sort(key=abs)
    qs.sort(key=abs)

    xmin = -0.5
    xmax = 0.5
    ymin = -0.5
    ymax = 0.5
    rs = []
    abss = []
    epsilon = 1e-8
    for q in qs:
        # Numerator determines quadrant.
        for p in ps:
            r = p / q
            if xmin <= r.real <= xmax and ymin <= r.imag <= ymax:
                # As the denominators are ordered by absolute value this captures only reduced forms.
                if all([abs(r - o) > epsilon for o in rs]):
                    rs.append(r)
                    abss.append(abs(q))

    print "Plotting %d rational complex numbers." % len(rs)

    fig = gcf()
    fig.gca().set_aspect(1)
    fig.gca().set_xlim(-0.5, 0.5)
    fig.gca().set_ylim(-0.5, 0.5)
    colors = cycle(['r', 'g', 'b', 'y', 'c', 'm'])
    old_a = None
    color = None
    for r, a in zip(rs, abss):
        if a != old_a:
            old_a = a
            color = next(colors)
        circle = Circle((r.real, r.imag), size * a ** -decay, color=color)
        fig.gca().add_artist(circle)
    show()