#435 Thomae’s function using python

import numpy as np
import matplotlib.pyplot as plt
xlist = []
ylist = []
for q in range(2,1000):
    for p in range(1,q):
        if np.gcd(p,q) == 1:
            xlist.append(p/q)
            ylist.append(1/q)

plt.plot(xlist, ylist, '.r', ms=3)
plt.xlabel('x')
plt.ylabel('y')
plt.xlim(0,1)
plt.grid('on')
plt.title("Thomae's function")
plt.savefig("ex435.png", dpi=100)
plt.show()
lim = 0.001
num = sum(y > lim for y in ylist)
print(f'Found {num} points above y={lim}')
lim = 0.01
num = sum(y > lim for y in ylist)
print(f'Found {num} points above y={lim}')
lim = 0.1
num = sum(y > lim for y in ylist)
print(f'Found {num} points above y={lim}')

Output:

Found 303791 points above y=0.001
Found 3003 points above y=0.01
Found 27 points above y=0.1