from math import gcd
import numpy as np
import matplotlib.pyplot as plt

## Thomae's function
xthom = []
ythom = []
for q in range(1, 101):
    for p in range(1, 3*q):
        xthom.append(p * 1.0 / q * 1.0)
        ythom.append(1 / (q / gcd(p, q)))

fig, axs = plt.subplots(1, 1, figsize=(5, 5))
axs.scatter(xthom, ythom, color="k", marker='.')
#axs.set_title("Negative Thomae's popcorn function")
fig.tight_layout()

# plt.show()
plt.savefig("popcorn.png")