Monsky's theorem: Coloring the plane

monsky.py

#!/usr/bin/python3

""" Illustrates the coloring of the plane used for Monsky's theorem. """


from PIL import Image
import argparse
import os
import sys


def get_cli_args():
    """ Get arguments form command line. """
    parser = argparse.ArgumentParser(description="Create Colormap.")
    parser.add_argument('-p', '--precision', type=int, default=100,
                        help='Number of divisions of the unit square, '
                             'default 100')
    parser.add_argument('-s', '--size', type=int, default=10,
                        help='Pixel size, default=10')
    parser.add_argument('-d', '--divisor', type=int, default=2,
                        help='Divisor, default: 2')
    parser.add_argument('-a', '--auto', action='store_true',
                        help='Automatic preview')
    parser.add_argument("-t", "--test", action="store_true",
                        help="Run doctest.")
    return parser.parse_args()


def divisor_multiplicity(integer, divisor=2):
    """Returns exponent of divisor in prime factorization of an integer.
    Brute force & super slow, but enough for our case.

    >>> divisor_multiplicity(100, 2)
    2
    >>> divisor_multiplicity(1024, 2)
    10
    """
    if not isinstance(integer, int) or integer == 0:
        raise ValueError("{} of type {} not allowed as argument".format(
            integer, type(integer)))
    mult = 0
    while True:
        if integer % divisor == 0:
            mult += 1
            integer /= divisor
        else:
            break
    return mult


def p_adic_value(a, b, divisor=2):
    """ Returns the $divisor-adic value of the fraction a/b, i.e.
    if a/b = divisor^n * c/d with all numbers integers and c/d coprime, then
    this returns divisor^(-n).

    >>> p_adic_value(1,10,2)
    2
    >>> p_adic_value(10,1,2)
    0.5
    >>> p_adic_value(1024,13,2)
    0.0009765625
    >>> p_adic_value(13,1024,2)
    1024
    """
    if b == 0:
        raise ValueError("0 not allowed as argument.")
    elif a == 0:
        return 0
    else:
        return divisor**(divisor_multiplicity(b, divisor) -
                         divisor_multiplicity(a, divisor))


def pixel_color(xcood, ycood, size, divisor=2):
    """ Returns the color of a pixel in our picture.

    >>> pixel_color(10, 10, 1)
    (255, 0, 0)
    >>> pixel_color(10, 11, 1)
    (0, 255, 0)
    >>> pixel_color(1, 11, 1)
    (0, 0, 255)
    """
    red = (255, 0, 0)
    green = (0, 255, 0)
    blue = (0, 0, 255)

    v_x = p_adic_value(xcood, size, divisor)
    v_y = p_adic_value(ycood, size, divisor)

    if v_x >= v_y and v_x >= 1:
        return blue
    elif v_x < v_y and v_y >= 1:
        return green
    elif v_x < 1 and v_y < 1:
        return red
    else:
        raise Exception


def create_picture(pixel_size, precision, divisor):
    """ Create the picture. """
    img = Image.new('RGB', (precision, precision), "black")
    pixels = img.load()

    print("Calculating colors....")

    for i in range(img.size[0]):
        if i > 0 and i % 50 == 0:
            print("{}% completed.".format((i*100)//img.size[0]))
        for j in range(img.size[1]):
            pixels[i, j] = pixel_color(i, img.size[1] - 1 - j,
                                       precision, divisor)

    print("done....")
    print("saving image.....")

    img_sized = img.resize((pixel_size*precision, pixel_size*precision),
                           Image.NEAREST)
    img_sized.save('image.bmp', 'bmp')

    print("finished......")

if __name__ == '__main__':
    args = get_cli_args()

    if args.test:
        import doctest
        doctest.testmod(verbose=True)
        sys.exit(0)

    create_picture(args.size, args.precision, args.divisor)
    if args.auto:
        os.system("eog 'image.bmp' &")

see blog entry http://www.lieret.net/2017/10/24/monsky/