アップデート。
@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 法は、アルゴリズム信号の周期を決定するための理想的なツールのように思えますが、装置が処理できるサンプル数によってのみ制限されます。
、
明らかに、すべての素数に対して降順の数列が表示されます。