更新。
@tom10 对自相关进行了彻底的调查，并解释了为什么自相关中的第一个凸起可以给出周期性信号的周期。
我尝试了两种方法，FFT 和自相关。他们的结果同意，尽管我更喜欢 FFT 而不是自相关，因为它更直接地为您提供了周期。
使用自相关时，我们只需确定第一个峰值的坐标。手动检查自相关图将显示您是否有“正确”峰值，因为您可以注意到周期（尽管对于 7 以上的素数，这变得不太清楚）。我相信你也可以制定一个简单的算法来计算“正确”的峰值。也许有人可以详细说明一些可以完成这项工作的简单算法？
例如，请参见以下序列自相关旁边的图。
代码：
1 # Plotting sequences satisfying, x_{i+1} = p-1 - (p*i-1 mod x_i)
2 # with p prime and x_0 = 1, next to their autocorrelation.
3
4 from __future__ import print_function
5 import numpy as np
6 from matplotlib import pyplot as plt
7
8 # The length of the sequences.
9 seq_length = 10000
10
11 upperbound_primes = 12
12
13 # Computing a list of prime numbers up to n
14 def primes(n):
15 sieve = [True] * n
16 for i in xrange(3,int(n**0.5)+1,2):
17 if sieve[i]:
18 sieve[i*i::2*i]=[False]*((n-i*i-1)/(2*i)+1)
19 return [2] + [i for i in xrange(3,n,2) if sieve[i]]
20
21 # The list of prime numbers up to upperbound_primes
22 p = primes(upperbound_primes)
23
24 # The amount of primes numbers
25 no_primes = len(p)
26
27 # Generate the sequence for the prime number p
28 def sequence(p):
29 x = np.empty(seq_length)
30 x[0] = 1
31 for i in range(1,seq_length):
32 x[i] = p - 1 - (p * (i-1) - 1) % x[i-1]
33 return x
34
35 # List with the sequences.
36 seq = [sequence(i) for i in p]
37
38 # Autocorrelation function.
39 def autocor(x):
40 result = np.correlate(x,x,mode='full')
41 return result[result.size/2:]
42
43 fig = plt.figure("The sequences next to their autocorrelation")
44 plt.suptitle("The sequences next to their autocorrelation")
45
46 # Proper spacing between subplots.
47 fig.subplots_adjust(hspace=1.2)
48
49 # Set up pyplot to use TeX.
50 plt.rc('text',usetex=True)
51 plt.rc('font',family='serif')
52
53 # Maximize plot window by command.
54 mng = plt.get_current_fig_manager()
55 mng.resize(*mng.window.maxsize())
56
57 k = 0
58 for s in seq:
59 k = k + 1
60 fig.add_subplot(no_primes,2,2*(k-1)+1)
61 plt.title("Sequence of the prime %d" % p[k-1])
62 plt.plot(s)
63 plt.xlabel(r"Index $i$")
64 plt.ylabel(r"Sequence number $x_i$")
65 plt.xlim(0,100)
66
67 # Constrain the number of ticks on the y-axis, for clarity.
68 plt.locator_params(axis='y',nbins=4)
69
70 fig.add_subplot(no_primes,2,2*k)
71 plt.title(r"Autocorrelation of the sequence $^{%d}x$" % p[k-1])
72 plt.plot(autocor(s))
73 plt.xlabel(r"Index $i$")
74 plt.xticks
75 plt.ylabel("Autocorrelation")
76
77 # Proper scaling of the y-axis.
78 ymin = autocor(s)[1]-int(autocor(s)[1]/10)
79 ymax = autocor(s)[1]+int(autocor(s)[1]/10)
80 plt.ylim(ymin,ymax)
81 plt.xlim(0,500)
82
83 plt.locator_params(axis='y',nbins=4)
84
85 # Use scientific notation when 0< t < 1 or t > 10
86 plt.ticklabel_format(style='sci',axis='y',scilimits=(0,1))
87
88 plt.show()
使用 FFT 时，我们对序列进行傅里叶变换并寻找第一个峰值。第一个峰值的坐标给出了代表我们信号的最粗略的频率。这将给出我们的周期，因为最粗略的频率是我们的序列（理想情况下）振荡的频率。
请参阅下面的傅立叶变换序列图。
代码：
1 # Plotting sequences satisfying, x_{i+1} = p-1 - (p*i-1 mod x_i)
2 # with p prime and x_0 = 1, next to their Fourier transforms.
3
4 from __future__ import print_function
5 import numpy as np
6 from matplotlib import pyplot as plt
7
8 # The length of the sequences.
9 seq_length = 10000
10
11 upperbound_primes = 12
12
13 # Computing a list of prime numbers up to n
14 def primes(n):
15 sieve = [True] * n
16 for i in xrange(3,int(n**0.5)+1,2):
17 if sieve[i]:
18 sieve[i*i::2*i]=[False]*((n-i*i-1)/(2*i)+1)
19 return [2] + [i for i in xrange(3,n,2) if sieve[i]]
20
21 # The list of prime numbers up to upperbound_primes
22 p = primes(upperbound_primes)
23
24 # The amount of primes numbers
25 no_primes = len(p)
26
27 # Generate the sequence for the prime number p
28 def sequence(p):
29 x = np.empty(seq_length)
30 x[0] = 1
31 for i in range(1,seq_length):
32 x[i] = p - 1 - (p * (i-1) - 1) % x[i-1]
33 return x
34
35 # List with the sequences.
36 seq = [sequence(i) for i in p]
37
38 fig = plt.figure("The sequences next to their FFT")
39 plt.suptitle("The sequences next to their FFT")
40
41 # Proper spacing between subplots.
42 fig.subplots_adjust(hspace=1.2)
43
44 # Set up pyplot to use TeX.
45 plt.rc('text',usetex=True)
46 plt.rc('font',family='serif')
47
48
49 # Maximize plot window by command.
50 mng = plt.get_current_fig_manager()
51 mng.resize(*mng.window.maxsize())
52
53 k = 0
54 for s in seq:
55 f = np.fft.rfft(s)
56 f[0] = 0
57 freq = np.fft.rfftfreq(seq_length)
58 k = k + 1
59 fig.add_subplot(no_primes,2,2*(k-1)+1)
60 plt.title("Sequence of the prime %d" % p[k-1])
61 plt.plot(s)
62 plt.xlabel(r"Index $i$")
63 plt.ylabel(r"Sequence number $x_i$")
64 plt.xlim(0,100)
65
66 # Constrain the number of ticks on the y-axis, for clarity.
67 plt.locator_params(nbins=4)
68
69 fig.add_subplot(no_primes,2,2*k)
70 plt.title(r"FFT of the sequence $^{%d}x$" % p[k-1])
71 plt.plot(freq,abs(f))
72 plt.xlabel("Frequency")
73 plt.ylabel("Amplitude")
74 plt.locator_params(nbins=4)
75
76 # Use scientific notation when 0 < t < 0 or t > 10
77 plt.ticklabel_format(style='sci',axis='y',scilimits=(0,1))
78
79 plt.show()
要了解为什么 FFT 方法比自相关方法更方便，请注意我们有一个明确的确定周期的算法：找到傅立叶变换的第一个峰值。对于足够数量的样本，这总是有效的。
见下表，由 FFT 方法得到，与自相关方法一致。
prime frequency period
2 0.00 1000.00
3 0.50 2.00
5 0.08 12.00
7 0.02 59.88
11 0.00 1000.00
下面的代码实现了该算法，打印一个表格，指定每个素数序列的频率和周期。
1 # Print a table of periods, determined by the FFT method,
2 # of sequences satisfying,
3 # x_{i+1} = p-1 - (p*i-1 mod x_i) with p prime and x_0 = 1.
4
5 from __future__ import print_function
6 import numpy as np
7 from matplotlib import pyplot as plt
8
9 # The length of the sequences.
10 seq_length = 10000
11
12 upperbound_primes = 12
13
14 # Computing a list of prime numbers up to n
15 def primes(n):
16 sieve = [True] * n
17 for i in xrange(3,int(n**0.5)+1,2):
18 if sieve[i]:
19 sieve[i*i::2*i]=[False]*((n-i*i-1)/(2*i)+1)
20 return [2] + [i for i in xrange(3,n,2) if sieve[i]]
21
22 # The list of prime numbers up to upperbound_primes
23 p = primes(upperbound_primes)
24
25 # The amount of primes numbers
26 no_primes = len(p)
27
28 # Generate the sequence for the prime number p
29 def sequence(p):
30 x = np.empty(seq_length)
31 x[0] = 1
32 for i in range(1,seq_length):
33 x[i] = p - 1 - (p * (i-1) - 1) % x[i-1]
34 return x
35
36 # List with the sequences.
37 seq = [sequence(i) for i in p]
38
39 # Function that finds the first peak.
40 # Assumption: seq_length >> 10 so the Fourier transformed
41 # signal is sufficiently smooth.
42 def firstpeak(x):
43 for i in range(10,len(x)-1):
44 if x[i+1] < x[i]:
45 return i
46 return len(x)-1
47
48 k = 0
49 for s in seq:
50 f = np.fft.rfft(s)
51 freq = np.fft.rfftfreq(seq_length)
52 k = k + 1
53 if k == 1:
54 print("prime \t frequency \t period")
55 print(p[k-1],'\t %.2f' % float(freq[firstpeak(abs(f))]), \
56 '\t\t %.2f' % float(1/freq[firstpeak(abs(f))]))
我在上述所有代码中使用了 10000 个样本（seq_length）。随着我们增加样本数量，可以看到周期收敛到某个整数值（使用 FFT 方法）。
在我看来，FFT 方法是确定算法信号周期的理想工具，它仅受设备可以处理的样本数量的限制。
,
很明显，我们看到所有素数的数字降序排列。
