Hello! I wrote this post. I know your comment was a lighthearted one, but I found it interesting, so I just wanted to share that it is quite straightforward to extend the solution to numbers divisible by 7. The indicator functions-based solution simply becomes:

for n in range(1, 1001):
    s = [n, 'Fizz', 'Buzz', 'FizzBuzz', 'Jazz', 'FizzJazz', 'BuzzJazz', 'FizzBuzzJazz']
    i = (n % 3 == 0) + 2 * (n % 5 == 0) + 4 * (n % 7 == 0)
    print(s[i])

For the cosine-based solution, we only need to translate the additional term 4 * (n % 7 == 0) to 4 * (1/7) * sum_{k=0}^6 e^(2 * pi * i * k * n / 7) = 4 * (1/7 + (2/7) * cos(2 * pi * n / 7) + (2/7) * cos(4 * pi * n / 7) + (2/7) * cos(6 * pi * n / 7)). When we add it to the existing solution, we get this:

from math import cos, pi
for n in range(1, 1001):
    s = [n, 'Fizz', 'Buzz', 'FizzBuzz', 'Jazz', 'FizzJazz', 'BuzzJazz', 'FizzBuzzJazz']
    i = round((137 / 105) + (2 / 3) * cos(2 * pi * n / 3)
                          + (4 / 5) * cos(2 * pi * n / 5)
                          + (4 / 5) * cos(4 * pi * n / 5)
                          + (8 / 7) * cos(2 * pi * n / 7)
                          + (8 / 7) * cos(4 * pi * n / 7)
                          + (8 / 7) * cos(6 * pi * n / 7))
    print(s[i])

It is worth noting that i can be viewed as a flag in which each bit corresponds to whether Fizz, Buzz or Jazz appears in the nth term of the sequence. All other details follow from this idea.