Gráfica de la función de Thomae

thomae.py

import matplotlib.pyplot as plt

def mcd(a, b):
    'Algoritmo de Euclides para hallar el máximo común divisor.'
    while b != 0:
        (a, b) = (b, a % b)
    return a

def grafica_thomae(a=0, b=1, denominador_max=1000):
    '''Aproximación de la gráfica de la función de Thomae en el
    intervalo [a, b], con a < b enteros.'''
    for q in range(1, denominador_max + 1):
        for p in range(a*q, b*q + 1):
            if mcd(p, q) == 1:
                yield (p/q, 1/q)    

if __name__ == '__main__':
    X, Y = zip(*grafica_thomae(-1, 1))
    fig, ax = plt.subplots(figsize=(16.18, 10))
    ax.plot(X, Y, '.')
    ax.spines['right'].set_visible(False)
    ax.spines['top'].set_visible(False)
    fig.savefig('thomae.png')