"Three-bonacchi" 序列
"Three-bonacchi" sequence
我们有一个包含前 3 个元素的序列：
t_1 = t_2 = t_3 = 1
序列的其余部分由规则定义：
t_n = t_(n-1) + t_(n-2) + t_(n-3)（类似于斐波那契数列，但有 3 个数字）。
t = {1; 1; 1; 3; 5; 9; 17; 31; ...}
任务是找到不是序列中任何元素的除数的第 N 个奇数。
Input: N (1 <= N <= 10^4 )
Output: N-th number which satisfies the condition.
示例：
Input: 125
Output: 2025
我的直接解决方案太慢了。我应该如何 improve/change 在给定限制条件下 (N <= 10^4) 在 1 秒内完成工作的算法？
t = [1, 1, 1]
N = int(input()) # find N-th odd number which isn't a divider of any number in the sequence
counter = 0 # how many appropriate numbers we've already found
curr_number = 1 # number to check
for i in range(100000):
t.append(t[-1] + t[-2] + t[-3])
while counter < N:
curr_number += 2
for i in range(len(t)):
if t[i] % curr_number == 0:
break
else:
counter += 1
print(curr_number)
如 Project Euler 中所述 description 这个问题：
It can be shown that 27 does not divide any terms of this sequence.
为了证明这一点，您显然不会将 tribonacci 序列计算为无穷大来检查 27 不会整除任何数字。必须有一个数学捷径来证明这一点，如果我们能找到这个捷径，我们就可以用它来检查其他数字是否能划分 tribonacci 数列。
检查一个数是否被 27 除与检查一个数是否以 27 为模等于 0 是一样的。
如果我们对 tribonacci 序列取模 27，我们得到：
1 % 27 = 1
1 % 27 = 1
1 % 27 = 1
3 % 27 = 3
5 % 27 = 5
9 % 27 = 9
17 % 27 = 17
31 % 27 = 4
57 % 27 = 3
105 % 27 = 24
193 % 27 = 4
...
你会注意到，为了求出 193 % 27 = 4，我们不需要使用数字 193（因为它等于 31 + 57 + 105），我们可以只使用前三个数字：
(4 + 3 + 24) % 27 = 4
这意味着我们不需要实际的 tribonacci 序列来检查 27 是否整除它。我们只需要查看模的序列来检查我们是否找到零：
1 % 27 = 1
1 % 27 = 1
1 % 27 = 1
( 1 + 1 + 1) % 27 = 3
( 1 + 1 + 3) % 27 = 5
( 1 + 3 + 5) % 27 = 9
( 3 + 5 + 9) % 27 = 17
( 5 + 9 + 17) % 27 = 4
( 9 + 17 + 4) % 27 = 3
(17 + 4 + 3) % 27 = 24
( 4 + 3 + 24) % 27 = 4
...
由于这个序列只包含27以下的数字，任何三个连续数字的可能性都是有限的，并且在某个时候会出现三个连续的数字，这些数字已经出现在序列的前面。
任何三个特定数字总是会产生相同的第四个数字，这意味着如果重复三个连续数字的组合，则将重复此后的整个序列。所以如果没有找到零，并且序列开始重复，你知道永远不会有零，并且数字不会整除 tribonacci 序列。
同样需要注意的是，序列中任意三个连续的数字只能是特定前一个数字的结果；例如序列 [3, 24, 4] 只能以数字 4 开头。这意味着重复将从序列的开头开始。因此，要在找到零之前检查序列是否重复，我们只需要找到前三个数字的重复：[1, 1, 1]。
这也意味着我们在计算时不必存储整个序列，我们可以继续用 [b, c, (a + b + c) % 替换 [a, b, c] n]，并将它们与 [1, 1, 1] 进行比较。
以27为例，序列在117个数字后重复：
1, 1, 1, 3, 5, 9, 17, 4, 3, 24, 4, 4, 5, 13, 22, 13, 21, 2, 9, 5, 16, 3, 24, 16, 16, 2, 7, 25, 7, 12, 17, 9, 11, 10, 3, 24, 10, 10, 17, 10, 10, 10, 3, 23, 9, 8, 13, 3, 24, 13, 13, 23, 22, 4, 22, 21, 20, 9, 23, 25, 3, 24, 25, 25, 20, 16, 7, 16, 12, 8, 9, 2, 19, 3, 24, 19, 19, 8, 19, 19, 19, 3, 14, 9, 26, 22, 3, 24, 22, 22, 14, 4, 13, 4, 21, 11, 9, 14, 7, 3, 24, 7, 7, 11, 25, 16, 25, 12, 26, 9, 20, 1, 3, 24, 1, 1, 26, 1, 1, 1 ...
所以检查数字 n 是否整除 tribonacci 序列的算法是：
- Start with the numbers a = 1, b = 1, c = 1
- Calucate d = (a + b + c) % n
- If d = 0 return true (n divides a number from the tribonacci sequence)
- Set a = b, b = c, c = d
- If a = 1 and b = 1 and c = 1 return false (beginning of repetition found)
- Repeat with new values for a, b and c
此代码示例适用于任何 N 值，但它显然不够快，无法在不到一秒的时间内找到第 10000 个奇数 non-dividing 数字（即 134241）。
var N = 125, n = 1, count = 0;
while (count < N) {
var a = 1, b = 1, c = 1, d;
n += 2;
while (d = (a + b + c) % n) {
a = b; b = c; c = d;
if (a == 1 && b == 1 && c == 1) {
++count;
break;
}
}
}
document.write(N + ": " + n);
我发现第一个零总是在第一个相同的三元组 [a=b=c] 之前，而不是在 [1,1,1] 之前，所以你可以将测试更改为 a == b && b == c 以让它 运行 快三倍。
var N = 125, n = 1, count = 0;
while (count < N) {
var a = 1, b = 1, c = 1, d;
n += 2;
while (d = (a + b + c) % n) {
a = b; b = c; c = d;
if (a == b && b == c) {
++count;
break;
}
}
}
document.write(N + ": " + n);
但即使使用 C 而不是 JavaScript，用这种方法找到第 10000 个奇数 non-dividing 也需要几分钟而不是几秒钟。可以使用@rici 的答案中的筛子思想进行进一步的改进，但如果真的有可能将其降低到一秒以下，那么一定有一个额外的数学捷径仍然在躲避我们。
下面的代码示例使用增量筛，因此它不需要有预定义的大小，并且可以用于任何 N 值。当测试一个值 n 并发现它是一个 non-divider，它的第一个奇数倍数 3×n 被设置为值 n，或者如果它已经被标记为另一个 non-divider 的倍数，则 5×n 或 7×n 或 ... 被设置为 n .当一个值 n 被认为在筛子中标记为 non-divider 的倍数时，标记将移动到 non-divider 的下一个奇数倍数。
function Sieve() { // incremental sieve
this.array = []; // associative array
}
Sieve.prototype.add = function(n) {
var base = n;
while (this.array[n += (2 * base)]); // find first unset odd multiple of n
this.array[n] = base; // set to base value
}
Sieve.prototype.check = function(n) {
var base = this.array[n]; // get base value
if (! base) return false; // if not set, return
delete this.array[n]; // delete current multiple
while (this.array[n += (2 * base)]); // find next unset odd multiple
this.array[n] = base; // set to base value
return true;
}
function dividesTribonacci(n) {
var a = 1, b = 1, c = 1, d;
while (d = (a + b + c) % n) {
a = b; b = c; c = d;
if (a == b && b == c) return false; // identical triple found
}
return true; // zero found, n divides tribonacci
}
function NthOddNonDivider(N) {
var n = 1, count = 0, sieve = new Sieve();
while (count < N) {
while (sieve.check(n += 2)) { // skip multiples of non-dividers
if (++count == N) return n;
}
if (! dividesTribonacci(n)) {
++count;
sieve.add(n);
}
}
return n;
}
document.write(NthOddNonDivider(125)); // 2025
我不知道如何在不到一秒的时间内解决 N=104 的问题，至少没有并行处理。但是我用至少在附近的商用笔记本电脑在 10 秒内完成了。
天真的 brute-force 解决方案简单地 运行 通过系列对每个奇数取模依次寻找 0 或重复的起始配置对于较小的 N 值当然足够了，但它的时间复杂度（显然）是 O(N3);在我的笔记本电脑上计算前 10,000 non-divisors 大约需要 22 分钟。问题是许多模数的周期长度为 O(N2)，因此一些试验涉及数亿步。即使内部循环是微不足道的，当你执行数十万次计算时，每次计算数亿步，也需要一段时间。
但是有一种非常有效的方法来修剪搜索，从简单的观察开始，如果 k 不除任何 tribonacci 数，那么 3k 或 5k 也不会或者，实际上，k 的任何其他倍数。 [注1]
所以一个简单的筛子就可以避免大量的试验。每次我们找到一个 non-divisor，我们将它的所有倍数（达到某个限制）标记为 non-divisors，然后我们就不必研究这些模数的周期了。这个简单的解决方案使计算时间减少了两位小数，从 1320 秒减少到 10 秒。
此外，使用筛子，观察到的时间复杂度从 O(N3) 变为 O(N2)。这有几个原因，我还没有做足够的分析来提供明确的证据，所以我不知道随着 N 变大，复杂度是否会保持二次方。 （事实上 ，我什至没有证据证明朴素算法的复杂性是三次的；如果没有某种数值分析，就不可能说某些模数 i 的 Tribonacci 周期小于 i<sup>3</sup>，这将使测试 n 模数为 O(n4)，但在测试之后几十万个模数，我也没找到周期长度接近模数平方的。
（我打算用 π(k) 表示 tribonacci 序列的周期长度模 k；这个函数的通常名称是 Pisano period，π(k) 是一个常见的符号。）
Chinese remainder theorem is sufficient to show that if j and k are coprime (i.e. they have no prime factors in common), then π(jk) = LCM(π(j), π(k)) (LCM is the Least Common Multiple)。这意味着我们可以很容易地计算 π(k) 给定 k 的质因数分解，如果我们对所有素数和素数的幂有 π(p) 的值。但是，对于此算法，重要的一点是复合数往往具有较大的周期，因为与最小公倍数组合的周期长度往往具有很少的公因子，因此 LCM 近似于乘积。此外，如果 p 是素数，则 π(p<sup>e</sup>) 通常是 p<sup> e−1</sup>×π(p)。 （它始终是 p<sup>e−1</sup>×π(p) 的约数，但在绝大多数情况下等式成立。）
筛子在过滤掉检查化合物数方面相当好；一旦找到 non-divisor 的幂，它还会过滤掉素数的幂。最终结果是实际调查的大部分模数的周期较短。特别是，由于我还不明白的原因，所有大多数non-divisor素数的周期（据我能够测试，唯一的例外是1,574,029）更少而不是素数本身，因此可以快速找到它。并且由于在任何包含零的序列中通常有很多零，因此也可以相当快地找到除数（尽管有时需要检查序列模 k 的前 k 个值以外的值, 它永远不会接近 k².)
所有这些都意味着在筛选 non-divisors 的已知倍数和当模数显示为除数时停止扫描之间，很少进行大周期扫描；在计算前10000个non-divisors时，单个模数的最长扫描只有几百万
备注
- 对于除数则不能这样说：
k 确实能整除一些 tribonacci 数这一事实并没有说明 k 的倍数是否是除数。
我们有一个包含前 3 个元素的序列：
t_1 = t_2 = t_3 = 1
序列的其余部分由规则定义：
t_n = t_(n-1) + t_(n-2) + t_(n-3)（类似于斐波那契数列，但有 3 个数字）。
t = {1; 1; 1; 3; 5; 9; 17; 31; ...}
任务是找到不是序列中任何元素的除数的第 N 个奇数。
Input: N (1 <= N <= 10^4 )
Output: N-th number which satisfies the condition.
示例：
Input: 125
Output: 2025
我的直接解决方案太慢了。我应该如何 improve/change 在给定限制条件下 (N <= 10^4) 在 1 秒内完成工作的算法？
t = [1, 1, 1]
N = int(input()) # find N-th odd number which isn't a divider of any number in the sequence
counter = 0 # how many appropriate numbers we've already found
curr_number = 1 # number to check
for i in range(100000):
t.append(t[-1] + t[-2] + t[-3])
while counter < N:
curr_number += 2
for i in range(len(t)):
if t[i] % curr_number == 0:
break
else:
counter += 1
print(curr_number)
如 Project Euler 中所述 description 这个问题：
It can be shown that 27 does not divide any terms of this sequence.
为了证明这一点，您显然不会将 tribonacci 序列计算为无穷大来检查 27 不会整除任何数字。必须有一个数学捷径来证明这一点，如果我们能找到这个捷径，我们就可以用它来检查其他数字是否能划分 tribonacci 数列。
检查一个数是否被 27 除与检查一个数是否以 27 为模等于 0 是一样的。
如果我们对 tribonacci 序列取模 27，我们得到：
1 % 27 = 1
1 % 27 = 1
1 % 27 = 1
3 % 27 = 3
5 % 27 = 5
9 % 27 = 9
17 % 27 = 17
31 % 27 = 4
57 % 27 = 3
105 % 27 = 24
193 % 27 = 4
...
你会注意到，为了求出 193 % 27 = 4，我们不需要使用数字 193（因为它等于 31 + 57 + 105），我们可以只使用前三个数字：
(4 + 3 + 24) % 27 = 4
这意味着我们不需要实际的 tribonacci 序列来检查 27 是否整除它。我们只需要查看模的序列来检查我们是否找到零：
1 % 27 = 1
1 % 27 = 1
1 % 27 = 1
( 1 + 1 + 1) % 27 = 3
( 1 + 1 + 3) % 27 = 5
( 1 + 3 + 5) % 27 = 9
( 3 + 5 + 9) % 27 = 17
( 5 + 9 + 17) % 27 = 4
( 9 + 17 + 4) % 27 = 3
(17 + 4 + 3) % 27 = 24
( 4 + 3 + 24) % 27 = 4
...
由于这个序列只包含27以下的数字，任何三个连续数字的可能性都是有限的，并且在某个时候会出现三个连续的数字，这些数字已经出现在序列的前面。
任何三个特定数字总是会产生相同的第四个数字，这意味着如果重复三个连续数字的组合，则将重复此后的整个序列。所以如果没有找到零，并且序列开始重复，你知道永远不会有零，并且数字不会整除 tribonacci 序列。
同样需要注意的是，序列中任意三个连续的数字只能是特定前一个数字的结果；例如序列 [3, 24, 4] 只能以数字 4 开头。这意味着重复将从序列的开头开始。因此，要在找到零之前检查序列是否重复，我们只需要找到前三个数字的重复：[1, 1, 1]。
这也意味着我们在计算时不必存储整个序列，我们可以继续用 [b, c, (a + b + c) % 替换 [a, b, c] n]，并将它们与 [1, 1, 1] 进行比较。
以27为例，序列在117个数字后重复：
1, 1, 1, 3, 5, 9, 17, 4, 3, 24, 4, 4, 5, 13, 22, 13, 21, 2, 9, 5, 16, 3, 24, 16, 16, 2, 7, 25, 7, 12, 17, 9, 11, 10, 3, 24, 10, 10, 17, 10, 10, 10, 3, 23, 9, 8, 13, 3, 24, 13, 13, 23, 22, 4, 22, 21, 20, 9, 23, 25, 3, 24, 25, 25, 20, 16, 7, 16, 12, 8, 9, 2, 19, 3, 24, 19, 19, 8, 19, 19, 19, 3, 14, 9, 26, 22, 3, 24, 22, 22, 14, 4, 13, 4, 21, 11, 9, 14, 7, 3, 24, 7, 7, 11, 25, 16, 25, 12, 26, 9, 20, 1, 3, 24, 1, 1, 26, 1, 1, 1 ...
所以检查数字 n 是否整除 tribonacci 序列的算法是：
- Start with the numbers a = 1, b = 1, c = 1
- Calucate d = (a + b + c) % n
- If d = 0 return true (n divides a number from the tribonacci sequence)
- Set a = b, b = c, c = d
- If a = 1 and b = 1 and c = 1 return false (beginning of repetition found)
- Repeat with new values for a, b and c
此代码示例适用于任何 N 值，但它显然不够快，无法在不到一秒的时间内找到第 10000 个奇数 non-dividing 数字（即 134241）。
var N = 125, n = 1, count = 0;
while (count < N) {
var a = 1, b = 1, c = 1, d;
n += 2;
while (d = (a + b + c) % n) {
a = b; b = c; c = d;
if (a == 1 && b == 1 && c == 1) {
++count;
break;
}
}
}
document.write(N + ": " + n);
我发现第一个零总是在第一个相同的三元组 [a=b=c] 之前，而不是在 [1,1,1] 之前，所以你可以将测试更改为 a == b && b == c 以让它 运行 快三倍。
var N = 125, n = 1, count = 0;
while (count < N) {
var a = 1, b = 1, c = 1, d;
n += 2;
while (d = (a + b + c) % n) {
a = b; b = c; c = d;
if (a == b && b == c) {
++count;
break;
}
}
}
document.write(N + ": " + n);
但即使使用 C 而不是 JavaScript，用这种方法找到第 10000 个奇数 non-dividing 也需要几分钟而不是几秒钟。可以使用@rici 的答案中的筛子思想进行进一步的改进，但如果真的有可能将其降低到一秒以下，那么一定有一个额外的数学捷径仍然在躲避我们。
下面的代码示例使用增量筛，因此它不需要有预定义的大小，并且可以用于任何 N 值。当测试一个值 n 并发现它是一个 non-divider，它的第一个奇数倍数 3×n 被设置为值 n，或者如果它已经被标记为另一个 non-divider 的倍数，则 5×n 或 7×n 或 ... 被设置为 n .当一个值 n 被认为在筛子中标记为 non-divider 的倍数时，标记将移动到 non-divider 的下一个奇数倍数。
function Sieve() { // incremental sieve
this.array = []; // associative array
}
Sieve.prototype.add = function(n) {
var base = n;
while (this.array[n += (2 * base)]); // find first unset odd multiple of n
this.array[n] = base; // set to base value
}
Sieve.prototype.check = function(n) {
var base = this.array[n]; // get base value
if (! base) return false; // if not set, return
delete this.array[n]; // delete current multiple
while (this.array[n += (2 * base)]); // find next unset odd multiple
this.array[n] = base; // set to base value
return true;
}
function dividesTribonacci(n) {
var a = 1, b = 1, c = 1, d;
while (d = (a + b + c) % n) {
a = b; b = c; c = d;
if (a == b && b == c) return false; // identical triple found
}
return true; // zero found, n divides tribonacci
}
function NthOddNonDivider(N) {
var n = 1, count = 0, sieve = new Sieve();
while (count < N) {
while (sieve.check(n += 2)) { // skip multiples of non-dividers
if (++count == N) return n;
}
if (! dividesTribonacci(n)) {
++count;
sieve.add(n);
}
}
return n;
}
document.write(NthOddNonDivider(125)); // 2025
我不知道如何在不到一秒的时间内解决 N=104 的问题，至少没有并行处理。但是我用至少在附近的商用笔记本电脑在 10 秒内完成了。
天真的 brute-force 解决方案简单地 运行 通过系列对每个奇数取模依次寻找 0 或重复的起始配置对于较小的 N 值当然足够了，但它的时间复杂度（显然）是 O(N3);在我的笔记本电脑上计算前 10,000 non-divisors 大约需要 22 分钟。问题是许多模数的周期长度为 O(N2)，因此一些试验涉及数亿步。即使内部循环是微不足道的，当你执行数十万次计算时，每次计算数亿步，也需要一段时间。
但是有一种非常有效的方法来修剪搜索，从简单的观察开始，如果 k 不除任何 tribonacci 数，那么 3k 或 5k 也不会或者，实际上，k 的任何其他倍数。 [注1]
所以一个简单的筛子就可以避免大量的试验。每次我们找到一个 non-divisor，我们将它的所有倍数（达到某个限制）标记为 non-divisors，然后我们就不必研究这些模数的周期了。这个简单的解决方案使计算时间减少了两位小数，从 1320 秒减少到 10 秒。
此外，使用筛子，观察到的时间复杂度从 O(N3) 变为 O(N2)。这有几个原因，我还没有做足够的分析来提供明确的证据，所以我不知道随着 N 变大，复杂度是否会保持二次方。 （事实上 ，我什至没有证据证明朴素算法的复杂性是三次的；如果没有某种数值分析，就不可能说某些模数 i 的 Tribonacci 周期小于 i<sup>3</sup>，这将使测试 n 模数为 O(n4)，但在测试之后几十万个模数，我也没找到周期长度接近模数平方的。
（我打算用 π(k) 表示 tribonacci 序列的周期长度模 k；这个函数的通常名称是 Pisano period，π(k) 是一个常见的符号。）
Chinese remainder theorem is sufficient to show that if j and k are coprime (i.e. they have no prime factors in common), then π(jk) = LCM(π(j), π(k)) (LCM is the Least Common Multiple)。这意味着我们可以很容易地计算 π(k) 给定 k 的质因数分解，如果我们对所有素数和素数的幂有 π(p) 的值。但是，对于此算法，重要的一点是复合数往往具有较大的周期，因为与最小公倍数组合的周期长度往往具有很少的公因子，因此 LCM 近似于乘积。此外，如果 p 是素数，则 π(p<sup>e</sup>) 通常是 p<sup> e−1</sup>×π(p)。 （它始终是 p<sup>e−1</sup>×π(p) 的约数，但在绝大多数情况下等式成立。）
筛子在过滤掉检查化合物数方面相当好；一旦找到 non-divisor 的幂，它还会过滤掉素数的幂。最终结果是实际调查的大部分模数的周期较短。特别是，由于我还不明白的原因，所有大多数non-divisor素数的周期（据我能够测试，唯一的例外是1,574,029）更少而不是素数本身，因此可以快速找到它。并且由于在任何包含零的序列中通常有很多零，因此也可以相当快地找到除数（尽管有时需要检查序列模 k 的前 k 个值以外的值, 它永远不会接近 k².)
所有这些都意味着在筛选 non-divisors 的已知倍数和当模数显示为除数时停止扫描之间，很少进行大周期扫描；在计算前10000个non-divisors时，单个模数的最长扫描只有几百万
备注
- 对于除数则不能这样说：
k确实能整除一些 tribonacci 数这一事实并没有说明k的倍数是否是除数。