Summary
=======
This answer discusses four cases:
- limit $f/f_s \to 0$,
- synchronous sampling of a cosine with rational $f/f_s$,
- synchronous sampling of a sinusoid of arbitrary phase with rational $f/f_s$, and
- asynchronous sampling of a sinusoid with rational $f/f_s$,
where $f$ is the sinusoid frequency and $f_s$ is the sampling frequency.
---
Limit $f/f_s \to 0$
===================
Let the zero-to-peak amplitude of the continuous-time *prototypical sinusoid*, which is being sampled and quantized, be an integer $A$. In the limit $f/f_s \to 0,$ where $f$ is the sinusoid frequency and $f_s$ is the sampling frequency, the amplitudes of the harmonics are found in the [Fourier series][1] representation of the continuous-time piece-wise constant quantized waveform. The terms of the Fourier series are directly the harmonics of the fundamental frequency. We can construct the waveform as a sum of components like the red curve here, illustrated for sinusoid amplitude $A = 7$:
[![enter image description here][2]][2]<br>*Figure 1. prototypical sinusoid, of amplitude $A=7$, and its quantization in blue, with a component highlighted in red.*
We are working with cosine instead of a sine because the math is nicer this way. Because of the symmetries, each component has only odd harmonics and only cosine terms in its Fourier series. The non-zero coefficients of the Fourier series are given by:
$$a_n = \frac{4}{\pi}\int_{x_0}^{x_1}\cos(nx)\,k\,dx = \frac{4k}{\pi n}\big(\sin(x_1n) - \sin(x_0n)\big),\quad n\text{ odd},\tag{1}$$
with integer amplitude $k \in 1\ldots A$ in range $x_0 < x < x_1$ in the first quarter-period of the cosine. I did not bother to explicitly write the contributions from the symmetrically arranged pieces in the other quarters, because they will contribute identically to the odd harmonic cosine terms in the Fourier series. Instead I simply included the implicit factor 4 in the equation.
The Fourier series of the full piece-wise constant waveform is a sum of the Fourier series of the components. The boundaries of the pieces that we need to include in the sum are:
$$\begin{gather}0 < x < \operatorname{acos}\left(\frac{A - 0.5}{A}\right),\quad \text{if }k = A,\\
\operatorname{acos}\left(\frac{k + 0.5}{A}\right) < x < \operatorname{acos}\left(\frac{k - 0.5}{A}\right),\quad \text{if }k \in 1\ldots A-1,\end{gather}\tag{2}$$
where $\operatorname{acos}\left(\frac{s}{A}\right)$ comes from solving $A\cos(x) = s$ in the first quadrant. The non-zero coefficients of the Fourier series of the full piece-wise constant waveform are then given by:
$$\begin{gather}a_n = \frac{4A}{\pi n}\sin\bigg(\operatorname{acos}\left(\tfrac{A - 0.5}{A}\right)n\bigg) +\\\frac{4}{\pi n}\sum_{k=1}^{A-1} \Bigg(k\sin\bigg(\operatorname{acos}\left(\tfrac{k - 0.5}{A}\right)n\bigg) - k\sin\bigg(\operatorname{acos}\left(\tfrac{k + 0.5}{A}\right)n\bigg)\Bigg),\quad n\text{ odd}.\end{gather}\tag{3}$$
16-bit quantization
-------------------
It is not out of reach to compute Eq. 3 for $A=2^{23} - 1$ and $n \in {5, 7},$ as parameterized in the question, using something like Python's `mpmath`. Let's try that but with a 16-bit sinusoid first:
import mpmath as mp
def a_n(A, n):
if int(n) & 1:
return 4*A/(mp.pi*n)*mp.sin(mp.acos((A-0.5)/A)*n) + 4/(mp.pi*n)*mp.nsum(lambda k: k*( mp.sin(mp.acos((k-0.5)/A)*n) - mp.sin(mp.acos((k+0.5)/A)*n) ), [1, A-1])
else:
return 0
A = 2**15-1 # amplitude (integer)
n = 5 # number of the harmonic (odd integer)
mp.mp.prec = 53 # default precision
20*mp.log10(mp.fabs(a_n(A, n))/A) # Fourier series term amplitude in dB
mp.mp.prec = 106 # Compute again with double the default precision, see if we get the same result
20*mp.log10(mp.fabs(a_n(A, n))/A) # Fourier series term amplitude in dB
mp.mp.prec = 53 # Restore default precision
I'm only reporting as many digits as are agreed about by repeated computation using two different `mpmath` precision settings, including also the digits that got changed by $\pm1$ due to larger changes in later digits. For a 16-bit full-scale $A = 2^{15}-1$ sinusoid and $f/f_s \to 0$, the result is **-154.66449 dBFS** for the 5th harmonic and, after setting `n = 7` in the script, **-154.66506 dBFS** for the 7th harmonic, in about 15 seconds of computation on my PC for a single result with the higher precision setting. As a sanity check, I also computed the amplitude of the fundamental frequency (1st harmonic), which is interestingly 0.0000001605530 dBFS (or 1.605530E-7 dBFS in scientific notation) compared to the peak value $2^{15}-1$ at 0 dBFS. I think the rounding at the top of the sinusoid is "pulling it up", giving it a higher amplitude of the fundamental.
Let's have a closer look at the decay of the amplitudes of the odd harmonics:
for m in range(1, 20):
n = mp.mpf(10)**m+1 # number of the harmonic
mp.mp.prec = 70 # precision
a_n_low_prec = a_n(A, n)
mp.mp.prec = 90 # precision
a_n_high_prec = a_n(A, n)
print(str(n)+","+str(20*mp.log10(mp.fabs(a_n_low_prec)/A))+","+str(20*mp.log10(mp.fabs(a_n_high_prec)/A)))
[![enter image description here][3]][3]<br>*Figure 2. Amplitudes of a selection of harmonics of the quantized prototypical sinusoid, with $A=2^{15}-1$.*
The graph in Fig. 2 is not very smooth, and it would probably look more jagged if we were to plot all points. I think this is because the waveform is similar to a sawtooth riding on a sinusoid, and near the zero crossing the sawtooth frequency stays the same over a large portion of the waveform. Actually, the sawtooth frequency is $2\pi A$ times the fundamental frequency, which approximately equals that of the 2.05881E5th harmonic. Indeed, the amplitude of that harmonic turns out as -135.49405259 dBFS, seemingly much higher than that of a "random" early harmonic. It would also seem that the amplitudes of the harmonics start to roll off after this frequency, with the factor $1/n$ in Eq. 3 making the asymptotic decay -20 dB/decade (~-6 dB/octave).
24-bit quantization
-------------------
The following continuation of the earlier Python script computes in high precision your desired numbers for $A=2^{23} - 1$, in the limiting case $f/f_s \to 0$:
A = 2**23-1 # amplitude (integer)
n = 5 # number of the harmonic (odd integer)
mp.mp.prec = 70 # Precision 1
20*mp.log10(mp.fabs(a_n(A, n))/A) # in dBFS
mp.mp.prec = 140 # Precision 2
20*mp.log10(mp.fabs(a_n(A, n))/A) # in dBFS
n = 7 # number of the harmonic (odd integer)
mp.mp.prec = 70 # Precision 1
20*mp.log10(mp.fabs(a_n(A, n))/A) # in dB
mp.mp.prec = 140 # Precision 2
20*mp.log10(mp.fabs(a_n(A, n))/A) # in dB
mp.mp.prec = 53 # Restore default precision
After a total of a few hours of computation on my PC, the result is that with $A=2^{23} - 1$ and $f/f_s \to 0,$ the amplitude of the 5th harmonic is **-226.91150085 dBFS** and the amplitude of the 7th harmonic is **-226.9115030 dBFS**. I also computed the amplitude of the fundamental frequency, 3.9195785E-11 dBFS, and in a full-day computation the amplitudes of a number of higher harmonics, in Python:
A = 2**23-1 # amplitude (integer)
for m in range(1, 20):
n = mp.mpf(10)**m+1 # number of the harmonic
mp.mp.prec = 70 # precision
a_n_low_prec = a_n(A, n)
mp.mp.prec = 90 # precision
a_n_high_prec = a_n(A, n)
print(str(n)+","+str(20*mp.log10(mp.fabs(a_n_low_prec)/A))+","+str(20*mp.log10(mp.fabs(a_n_high_prec)/A)))
[![!\[enter image description here][4]][4]<br>*Figure 3. Amplitudes of a selection of harmonics of the quantized prototypical sinusoid, with $A=2^{23}-1$.*
The corner frequency before the roll-off is again approximately the sawtooth frequency at the zero crossing of the quantized prototypical sinusoid, $2\pi A = 2\pi(2^{23}-1) \approx$ 5.270717225E7 times the fundamental frequency $f$.
Synchronous sampling of a cosine with rational $f/f_s$
======================================================
Let the prototypical sinusoid, which is being sampled and quantized, be a cosine without a phase shift. That is to say, one of the samples is at the peak of the sinusoid. Let $f/f_s$ be equal to a rational number $c/d$, with $c$ and $d$ integer. The fundamental period $L$ of the [periodic sequence][5] of quantized samples is:
$$L = \frac{d}{\operatorname{gcd}(c, d)},\tag{4}$$
where gcd denotes the greatest common divisor. A real discrete Fourier transform (real-DFT) of length $L$ has bin frequencies that are the harmonic frequencies of the period-$L$ frequency. The period-$L$ frequency is equivalent to the frequency of the prototypical sinusoid only if the numerator of the [irreducible fraction][6] representation of $f/f_s$ equals 1. Otherwise the real-DFT bin frequencies contain also other frequencies than those equivalent to in-band harmonics of the sinusoid.
What we have learned from the results so far is that the amplitudes of the harmonics of the quantized prototypical sinusoid do not really start to decay until after the approximately $2\pi A$th harmonic. This means that with rational $f/f_s$, the amplitude of an in-band harmonic frequency of the quantized sequence is typically almost fully determined by aliased frequencies. As an example of such aliasing, with $f_s =$ 48 kHz and $f =$ 1000 Hz, the 5th harmonic gets summed with the aliases of the 43th, 53th, 91th, 101th, 139th, 149th, etc. harmonic:
$$\begin{eqnarray}
&& \ldots\\
&=& 48000\text{ Hz}\times -3 + 1000\text{ Hz}\times 149\\
&=& 48000\text{ Hz}\times -2 + 1000\text{ Hz}\times 101\\
&=& 48000\text{ Hz}\times -1 + 1000\text{ Hz}\times 53\\
&=& 1000\text{ Hz}\times 5\\
&=& 48000\text{ Hz}\times 1 - 1000\text{ Hz}\times 43\\
&=& 48000\text{ Hz}\times 2 - 1000\text{ Hz}\times 91\\
&=& 48000\text{ Hz}\times 3 - 1000\text{ Hz}\times 139\\
&=& \ldots\end{eqnarray}\tag{5}$$
The effect of summation of the harmonics aliased to a frequency depends also on the phase of the sinusoid being quantized, because the phases of the aliased harmonics control whether there is constructive or destructive interference. Aliasing does not change the phase of a harmonic. A phase shift of the quantized prototypical sinusoid is equivalent to shifting it in time, which accounts to a phase shift of each unaliased harmonic proportional to its number. In this answer we will only be sampling and quantizing a prototypical cosine without a phase shift, so to compute a single bin of the real-DFT of the discrete-time periodic sequence of samples, it would suffices to sum those $a$ coefficients of the Fourier series of the quantized prototypical sinusoid that alias to the equivalent frequency of the bin. If the numerator of $f/f_s$ expressed as an [irreducible fraction][6] equals 1, we have:
$$\operatorname{bin}_n = a_n + \sum_{k=1}^{\infty}\left(a_{kL-n} + a_{kL+n}\right),\tag{6}$$
where $a_n$ is given by Eq. 3 for odd $n$ and $a_n$ is zero for even $n$. With $f = $ 1000 Hz, $f_s = $ 48 kHz, on my PC it would take perhaps thousands of years to compute Eq. 6 to reasonable accuracy for $A=2^{23}-1$, and perhaps weeks for $A=2^{15}-1$. I also have difficulties evaluating the series even with $A=2^3-1$, with the sequence of partial sums oscillating. I think this difficulty arises from that with $A = 2^3 - 1$ and $f/f_s = 1/48$, the sequence before quantization contains values $-3.5$ and $3.5$ exactly at a discontinuity in the quantized prototypical sinusoid. In any case, with rational $f/f_s$, the real-DFT-based approach to determine the amplitudes of the harmonics is much more practical. It would also work with non-integer amplitude $A$, which Eq. 2 cannot handle in its current form.
24-bit quantization
-------------------
The following Python script computes using the real-DFT-based method the amplitudes of the harmonics of a synchronously sampled cosine for $A = 2^{23}-1$, $f =$ 1000 Hz, and $f_s =$ 48 kHz:
import numpy as np
from sympy import Rational
def bins(A, c, d):
L = Rational(d, np.gcd(c, d)) # Fundamental period (samples)
waveform = np.around(A*np.cos(np.dot(range(L), 2*np.pi*c/d)))
return np.fft.rfft(waveform)/(A*L/2)
A = 2**23-1 # amplitude (integer)
c = 1000 # numerator of f/f_s (positive integer)
d = 48000 # denominator of f/f_s (positive integer)
[20*np.log10(float(abs(x))) for x in bins(A, c, d)] # real-DFT magnitudes (dBFS)
While not required by the script, in this case $c/d = 1/48,$ so each real-DFT bin corresponds to a harmonic of the prototypical sinusoid. The fundamental frequency's amplitude measures as -1.75E-7 dBFS, the 5th harmonic as **-160.90 dBFS** and the 7th as **-160.75 dBFS**. Unfortunately I could not find a suitable multiple precision Fast Fourier Transform (FFT) library, but I think those three numbers came out in sufficient precision using NumPy's `rfft`. The precision is not sufficient to calculate the amplitude of say the 3rd harmonic seemingly at about -333 dBFS (Fig. 4).
[![enter image description here][7]][7]<br>*Figure 4. Magnitudes of the real-DFT bins of the quantized sequence of a synchronously sampled cosine with $A=2^{23}-1$, $f =$ 1000 Hz, and $f_s =$ 48 kHz. These results are only applicable when one of the samples is at the peak of the sinusoid, so not applicable to arbitrary phase shifts or asynchronous sampling.*
With this $f/f_s$ the fundamental period is $L = 48$. The computed amplitudes of the harmonics are also correct for phase shifts of the prototypical cosine that are multiples of $2\pi/L$, because the cosine is sampled also at those phases, or such a multiple plus $\pi$, because of the symmetry properties of cosine. In the presented case the compatible phase shifts of the cosine are multiples of $7.5°$.
We have encountered here a troublesome combination of $L = 48$ and $A = 2^{23}-1$. We have samples $A\cos(8\times2\pi/L) =$ $A\cos(40\times2\pi/L) =$ $A/2$ and $A\cos(16\times2\pi/L) =$ $A\cos(32\times2\pi/L) =$ $-A/2$, and as $A$ is odd, those samples will always have a fractional part of $1/2$. The above reported amplitudes of the harmonics may thus be dependent on the rounding mode, or, the computation of the phase and/or evaluation of the cosine may have given a small numerical error which has determined the direction of rounding. This may even break the symmetry properties of the sequence, affecting which amplitudes of the harmonics of the period-$L$ frequency are zero.
Method validation
=================
We can validate the two approaches against each other by choosing a pair of $A$ and $f/f_s$ that does not result in sampling of the quantized prototypical cosine at any of its discontinuities. Firstly, this removes having to choose arbitrarily which way a number exactly half-way between integers should be rounded, and secondly, should help in convergence of the infinite series in Eq. 6. Let's try $A = 8$ and $f/f_s = 1/48$. First, the real-DFT-based method, in Python:
import numpy as np
from sympy import Rational
def bins(A, c, d):
L = Rational(d, np.gcd(c, d)) # Fundamental period (samples)
waveform = np.around(A*np.cos(np.dot(range(L), 2*np.pi*c/d)))
return np.fft.rfft(waveform)/(A*L/2)
A = 8 # amplitude (integer)
c = 1000 # numerator of f/f_s (positive integer)
d = 48000 # denominator of f/f_s (positive integer)
[20*np.log10(float(abs(x))) for x in bins(A, c, d)] # real-DFT magnitudes (dBFS)
This results in the following amplitudes of the harmonics, printed in precision beyond their accuracy, in dBFS:
-inf, 0.005622747208892056, -inf, -331.19944825653926, -inf, -50.41376796795221, -inf, -35.41599672115829, -inf, -327.2878768807577, -inf, -38.14783193548751, -inf, -45.42857299685606, -inf, -324.8188873214554, -inf, -40.41625528549065, -inf, -36.12873038111221, -inf, -337.55682735433885, -inf, -33.04816436790989, -inf
Then, aliasing the harmonics of the quantized prototypical cosine by evaluating Eq. 6, with the series truncated to the first approximately $10^4\pi A$ terms, in Python:
import mpmath as mp
def a_n(A, n):
if int(n) & 1:
return 4*A/(mp.pi*n)*mp.sin(mp.acos((A-0.5)/A)*n) + 4/(mp.pi*n)*mp.nsum(lambda k: k*( mp.sin(mp.acos((k-0.5)/A)*n) - mp.sin(mp.acos((k+0.5)/A)*n) ), [1, A-1])
else:
return 0
mp.mp.prec = 100
L = Rational(d, np.gcd(c, d)) # Fundamental period (samples)
for n in range(L/2):
bin = a_n(A, n) + mp.nsum(lambda k: a_n(A, k*L - n) + a_n(A, k*L + n), [1, int(10000*mp.pi*A)])
float(20*mp.log10(mp.fabs(bin)/A)) # harmonic amplitude (dBFS)
mp.mp.prec = 53
which results in the following amplitudes, printed in precision beyond their accuracy, in dBFS:
-inf, 0.005622452922818016, -inf, -165.35811090601203, -inf, -50.41370472648242, -inf, -35.41599394709451, -inf, -155.34167605811237, -inf, -38.14782298030069, -inf, -45.42859376902997, -inf, -164.6586753088176, -inf, -40.41623886311205, -inf, -36.12872396082294, -inf, -149.78722963456818, -inf, -33.04815743325275
The results seem to agree within their numerical accuracy (Fig 6.).
[![enter image description here][8]][8]<br>*Figure 6. Comparison of amplitudes, in dBFS, of harmonics of a quantized cosine sampled synchronously, with $A = 8$ and $f/f_s = 1/48$, calculated by the real-DFT-based method (blue squares) and by summing coefficients $a_n$ of harmonics that alias to the same frequency (orange crosses).*
The 3rd harmonic amplitude calculated by truncating the sum in Eq. 6 converges very slowly as function of the number of terms (Fig. 6). This is probably the case also with the other low-amplitude harmonics. It does not seem feasible to reach the accuracy of the real-DFT-based method using the harmonic aliasing method.
[![enter image description here][9]][9]<br>*Figure 7. $\operatorname{bin}_3$ as function of the number of terms in truncation of the sum in Eq. 6.*
Synchronous sampling of a sinusoid of arbitrary phase with rational $f/f_s$
========================================================================
24-bit quantization
-------------------
The real-DFT-based method is simple to modify for arbitrary phase shifts of the prototypical cosine. For $A = 2^{23}-1$, $f =$ 1000 Hz, $f_s =$ 48 kHz, and a phase shift $\phi = 0.123$ in radians, the amplitudes of the harmonics of the fundamental period can be computed by the following Python script:
import numpy as np
from sympy import Rational
def bins(A, c, d, phi):
L = Rational(d, np.gcd(c, d)) # Fundamental period (samples)
waveform = np.around(A*np.cos(np.dot(range(L), 2*np.pi*c/d) + phi))
return np.fft.rfft(waveform)/(A*L/2)
A = 2**23-1 # amplitude (integer)
c = 1000 # numerator of f/f_s (positive integer)
d = 48000 # denominator of f/f_s (positive integer)
phi = 0.123 # phase shift of the cosine (radians)
[20*np.log10(float(abs(x))) for x in bins(A, c, d, phi)] # real-DFT magnitudes (dBFS)
[![enter image description here][10]][10]<br>*Figure 8. Magnitudes of the real-DFT bins of the quantized sequence of a synchronously sampled cosine with phase shift $\phi = 0.123$, $A=2^{23}−1$, $f=$ 1000 Hz, and $f_s$= 48 kHz.*
Asynchronous sampling of a sinusoid with rational $f/f_s$
=========================================================
Asynchronous sampling with rational $f/f_s$ can be modeled by assuming that $f/f_s$ is very close to a rational number. Then, locally the sequence of samples appears periodic with a computable harmonic spectrum. After sufficient time has elapsed, the phase of the prototypical sinusoid will have drifted to a different value and with a different local harmonic spectrum. For a random such sufficient time, the drift is a random variable between $0$ and $2\pi$. A reasonable probabilistic analog of the amplitude of a real-DFT bin of the periodic sequence is the expected value of the square of its absolute value, $\operatorname{E}[|\operatorname{bin}_n|^2]$, which can be expressed in the same dBFS scale. We can estimate its value using uniform statistical sampling over phase shifts $\phi$ in a non-redundant range $0 < \phi < 2\pi/L$, where $L$ is the fundamental period of the locally periodic sequence. The end points of the range shall be deliberately left out so that we don't have the previously encountered problem with rounding modes and such.
24-bit quantization
-------------------
$\operatorname{E}[|\operatorname{bin}_n|^2]$ is estimated for $A = 2^{23}-1$, $f =$ 1000 Hz, and $f_s =$ 48 kHz, by the following Python script, which also estimates $\max|\operatorname{bin}_n|$ over random phase shifts, representing the worst-case non-fundamental harmonic amplitudes in the local spectra:
import numpy as np
from sympy import Rational
def E_sq_abs_bins_and_E_max_sq_abs_bins(A, c, d, N):
L = Rational(d, np.gcd(c, d)) # Fundamental period (samples)
num_bins = int((L + 2)/2)
sum_sq_bins = np.zeros(num_bins)
max_sq_bins = np.zeros(num_bins)
for n in range(N):
phi = float((0.5 + n/N)*2*np.pi/L)
waveform = np.around(A*np.cos(np.dot(range(L), 2*np.pi*c/d) + phi))
sq_bins = np.square(np.absolute(np.fft.rfft(waveform)/(A*L/2)))
sum_sq_bins = sum_sq_bins + sq_bins
max_sq_bins = np.maximum(max_sq_bins, sq_bins)
return [sum_sq_bins/N, max_sq_bins]
A = 2**23-1 # amplitude (integer)
c = 1000 # numerator of f/f_s (positive integer)
d = 48000 # denominator of f/f_s (positive integer)
N = 2**12 # size of statistical sample (positive integer)
results = E_sq_abs_bins_and_E_max_sq_abs_bins(A, c, d, N)
[10*np.log10(float(x)) for x in results[0]] # est. expected value of real-DFT squared magnitudes (dBFS)
[10*np.log10(float(x)) for x in results[1]] # est. max real-DFT squared magnitudes (dBFS)
Unofortunately, the results convergence quite slowly as the size of the statistical sample `N` is increased, and the accuracy requirements stated in the question are not fulfilled in a reasonable time on my PC for $\operatorname{E}[|\operatorname{bin}_n|^2]$ or for $\max|\operatorname{bin}_n|$. The results (Fig. 9) calculated in about half and hour are that for these parameter settings, $\operatorname{E}[|\operatorname{bin}_n|^2]$ is about -157 dB and $\max|\operatorname{bin}_n|$ about -149 to -146 dBFS for non-fundamental harmonics. For the fundamental frequency, $\operatorname{E}[|\operatorname{bin}_1|^2]$ is in the order of 5E-10 or less, and $\max|\operatorname{bin}_1|$ is in the order of 2.5E-7 or higher.
[![enter image description here][11]][11]<br>*Figure 9. Estimated $\operatorname{E}[|\operatorname{bin}_n|^2]$ (blue squares) and $\max|\operatorname{bin}_n|$ (red lines) for local $\operatorname{bin}_n$ in slow-drifting asynchronous sampling, for $A = 2^{23}-1$, $f =$ 1000 Hz, and $f_s =$ 48 kHz, estimated by a uniform statistical sample of phase shifts with a sample size of $2^{12}$.*
The estimates show that asynchronous sampling gives a smoothing effect on the statistics of the local amplitudes of the harmonics, compared to synchronous sampling with a fixed phase shift of the prototypical cosine with respect to the sampling grid. However, with asynchronous sampling, sometimes the worst case local harmonic peaks occur, whereas with synchronous sampling the harmonic amplitudes are well controlled and the worst cases are avoided by a suitable choice of phase shift.
[1]: https://en.wikipedia.org/wiki/Fourier_series
[2]: https://i.sstatic.net/1a4gD.png
[3]: https://i.sstatic.net/opb6G.png
[4]: https://i.sstatic.net/Gh3Jt.png
[5]: https://en.wikipedia.org/wiki/Periodic_sequence
[6]: https://en.wikipedia.org/wiki/Irreducible_fraction
[7]: https://i.sstatic.net/npZyu.png
[8]: https://i.sstatic.net/M4Rrg.png
[9]: https://i.sstatic.net/eqGx8.png
[10]: https://i.sstatic.net/LgsBj.png
[11]: https://i.sstatic.net/pTNKN.png