是否有更 Pythonic 的方式来编码这种递归关系 re: OEIS A077947
Is there a more Pythonic way of coding this recurrance relation re: OEIS A077947
我正在写一篇关于 Jacobsthal 序列 (A001045) 以及如何将它们视为由一些不同的子序列组成的论文。我在 A077947 上发表了评论，指出了这一点，并包含了一个 python 程序。不幸的是，编写的程序还有很多不足之处，所以我当然想求助于 Stack，看看这里是否有人知道如何改进代码！
代码如下：
a=1
b=1
c=2
d=5
e=9
f=18
for x in range(0, 100):
print(a, b, c, d, e, f)
a = a+(36*64**x)
b = b+(72*64**x)
c = c+(144*64**x)
d = d+(288*64**x)
e = e+(576*64**x)
f = f+(1152*64**x)
我解释一下这背后的原因：
The sequence A077947 is generated by 6 digital root preserving sequences
stitched together; per the Python code these sequences initiate at the
seed values a-f. The number of iterations required to calculate a given
A077947 a(n) is ~n/6. The code when executed returns all the values for
A077947 up to range(x), or ~x*6 terms of A077947. I find the repeated
digital roots interesting as I look for periodic digital root preservation
within sequences as a method to identify patterns within data. For
example, digital root preserving sequences enable time series analysis of
large datasets when estimating true-or-false status for alarms in large IT
ecosystems that undergo maintenance (mod7 environments); such analysis is
also related to predicting consumer demand / patterns of behavior.
Appropriating those methods of analysis, carving A077947 into 6 digital
root preserving sequences was meant to reduce complexity; the Python code
reproduces A077947 across 6 "channels" with seed values a-f. This long
paragraph boils down to statement, "The digital roots of the terms of the
sequence repeat in the pattern (1, 1, 2, 5, 9, 9)." The bigger statement
is that all sequences whose digital roots repeat with a pattern can be
partitioned/separated into an equal number of distinct sequences and those
sequences can be calculated independently. There was a bounty related to
this sequence.
这段代码很难看，但如果不这样编码我似乎无法得到正确的答案；
我还没有想出如何将它写成一个函数，因为我似乎无法将重复值正确地存储在一个函数中。
当然，如果这产生了好的结果，我们希望 link 对 OEIS 参考文献进行讨论。
这是一个link序列：
https://oeis.org/A077947
这与您的代码的行为相同，而且可以说更漂亮。您可能会看到使魔法常量不那么随意的方法。
factors = [ 1, 1, 2, 5, 9, 18 ]
cofactors = [ 36*(2**n) for n in range(6) ]
for x in range(10):
print(*factors)
for i in range(6):
factors[i] = factors[i] + cofactors[i] * 64**x
要仅计算其中一个子序列，在迭代时保持 i 固定就足够了。
这是一种无需第二个 for 循环的替代方法：
sequences = [ 1, 1, 2, 5, 9, 18 ]
multipliers = [ 36, 72, 144, 288, 576, 1152 ]
for x in range(100):
print(*sequences)
sequences = [ s + m*64**x for s,m in zip(sequences,multipliers) ]
[编辑] 查看这些值我注意到这个特定的序列也可以通过以下方式获得：
N[i+1] = 2 * N[i] + (-1,0,1 循环)
或
N[i+1] = 2 * N[i] + i mod 3 - 1 （假设索引从零开始）
i N[i] [-1,0,1] N[i+1]
0 1 -1 --> 2*1 - 1 --> 1
1 1 0 --> 2*1 + 0 --> 2
2 2 1 --> 2*2 + 1 --> 5
3 5 -1 --> 2*5 - 1 --> 9
4 9 0 --> 2*9 + 0 --> 18
...
因此，生成序列的更简单的循环可以是：
n = 1
for i in range(100):
print(n)
n = 2*n + i % 3 - 1
使用 functools 的 reduce 函数可以使这个更加简洁：
from functools import reduce
sequence = reduce(lambda s,i: s + [s[-1]*2 + i%3 - 1],range(20),[1])
print(sequence)
>>> [1, 1, 2, 5, 9, 18, 37, 73, 146, 293, 585, 1170, 2341, 4681, 9362, 18725, 37449, 74898, 149797, 299593, 599186]
使用您的多渠道方法和我建议的公式，这将得出：
sequences = [ 1, 1, 2, 5, 9, 18 ]
multipliers = [ 36, 72, 144, 288, 576, 1152 ]
allSequences = reduce(lambda ss,x: ss + [[ s + m*64**x for s,m in zip(ss[-1],multipliers) ]],range(100),[sequences])
for seq in allSequences: print(*seq) # print 6 by 6
[EDIT2] 如果您的所有序列都将具有相似的模式（即起始通道、乘数和计算公式），您可以将此类序列的打印概括为一个函数，因此每个序列只需要一行：
def printSeq(calcNext,sequence,multipliers,count):
for x in range(count):
print(*sequence)
sequence = [ calcNext(x,s,m) for s,m in zip(sequence,multipliers) ]
printSeq(lambda x,s,m:s*2+m*64**x,[1,1,2,5,9,18],multipliers=[36,72,144,288,576,1152],count=100)
[EDIT3] 改进 printSeq 函数。
我相信您不会总是需要一组乘数来计算每个通道中的下一个值。对该函数的改进是为 lambda 函数提供通道索引而不是乘数。这将允许您在需要时使用乘数数组，但也会让您使用更通用的计算。
def printSeq(name,count,calcNext,sequence):
p = len(sequence)
for x in range(count):
print(name, x,":","\t".join(str(s) for s in sequence))
sequence = [ calcNext(x,s,c,p) for c,s in enumerate(sequence) ]
lambda 函数有 4 个参数，预计 return 指定通道的下一个序列值：
s : current sequence value for the channel
x : iteration number
c : channel index (zero based)
p : number of channels
因此，在公式中使用数组可以这样表示：
printSeq("A077947",100,lambda x,s,c,p: s + [36,72,144,288,576,1152][c] * 64**x, [1,1,2,5,9,18])
但您也可以使用基于频道索引（和频道数）的更通用的公式：
printSeq("A077947",100,lambda x,s,c,p: s + 9 * 2**(p*x+c+2), [1,1,2,5,9,18])
或（基于 2*S + i%3 - 1 的 6 个通道）：
printSeq("A077947",100,lambda x,s,c,p: 64*s + 9*(c%3*2 - (c+2)%3 - 1) ,[1,1,2,5,9,18])
printSeq("A077947",100,lambda x,s,c,p: 64*s + 9*[-3,1,2][c%3],[1,1,2,5,9,18])
我的理由是，如果你有一个函数可以根据序列中的当前索引和值计算下一个值，你应该能够定义一个跨步函数来计算 N 个索引的值更远。
给定 F(i,S[i]) --> i+1,S[i+1]
F2(i,S[i]) --> i+2,S[i+2] = F(F(i,S[i]))
F3(i,S[i]) --> i+3,S[i+3] = F(F(F(i,S[i])))
...
F6(i,S[i]) --> i+6,S[i+6] = F(F(F(F(F(F(i,S[i]))))))
...
Fn(i,S[i]) --> i+n,S[i+n] = ...
这将始终有效，不需要乘法器数组。大多数情况下，应该可以仅使用代数来简化 Fn。
例如 A001045 : F(i,S) = i+1, 2*S + (-1)**i
printSeq("A001045",20,lambda x,s,c,p: 64*s + 21*(-1)**(x*p+c),[0,1,1,3,5,11])
请注意，从第 3 个值开始，可以在不知道索引的情况下计算该序列中的下一个值：
A001045: F(S) = 2*S + 1 - 2*0**((S+1)%4)
这是一个使用生成器的替代版本：
def Jacobsthal():
roots = [1, 1, 2, 5, 9, 18]
x = 0
while True:
yield roots
for i in range(6):
roots[i] += 36 * 2**i * 64**x
x += 1
这里有一个安全的使用方法：
j = Jacobsthal()
for _ in range(10):
print(j.send(None))
从OEIS注解可以看出，只需要3个初始值和1次递归3个前面的序列元素a(n-1)、a(n-2)和a(n-3)。然后您就可以轻松获得该系列。
# a(n-1)+a(n-2)+2*a(n-3)
# start values: 1, 1, 2
a1, a2, a3 = 1,1,2 #initial vales
m3, m2, m1 = 1,1,2 #multipliers for a(n-1):m3 a(n-2):m2 and a(n-3):m1
print a1, a2, a3,
for i in range(16):
a1,a2,a3=a2,a3,a3*m3+a2*m2+a1*m1
print a3,
这给你 1 1 2 5 9 18 37 73 146 293 585 1170 2341 4681 9362 18725 ...
我正在写一篇关于 Jacobsthal 序列 (A001045) 以及如何将它们视为由一些不同的子序列组成的论文。我在 A077947 上发表了评论，指出了这一点，并包含了一个 python 程序。不幸的是，编写的程序还有很多不足之处，所以我当然想求助于 Stack，看看这里是否有人知道如何改进代码！
代码如下：
a=1
b=1
c=2
d=5
e=9
f=18
for x in range(0, 100):
print(a, b, c, d, e, f)
a = a+(36*64**x)
b = b+(72*64**x)
c = c+(144*64**x)
d = d+(288*64**x)
e = e+(576*64**x)
f = f+(1152*64**x)
我解释一下这背后的原因：
The sequence A077947 is generated by 6 digital root preserving sequences stitched together; per the Python code these sequences initiate at the seed values a-f. The number of iterations required to calculate a given A077947 a(n) is ~n/6. The code when executed returns all the values for A077947 up to range(x), or ~x*6 terms of A077947. I find the repeated digital roots interesting as I look for periodic digital root preservation within sequences as a method to identify patterns within data. For example, digital root preserving sequences enable time series analysis of large datasets when estimating true-or-false status for alarms in large IT ecosystems that undergo maintenance (mod7 environments); such analysis is also related to predicting consumer demand / patterns of behavior. Appropriating those methods of analysis, carving A077947 into 6 digital root preserving sequences was meant to reduce complexity; the Python code reproduces A077947 across 6 "channels" with seed values a-f. This long paragraph boils down to statement, "The digital roots of the terms of the sequence repeat in the pattern (1, 1, 2, 5, 9, 9)." The bigger statement is that all sequences whose digital roots repeat with a pattern can be partitioned/separated into an equal number of distinct sequences and those sequences can be calculated independently. There was a bounty related to this sequence.
这段代码很难看，但如果不这样编码我似乎无法得到正确的答案；
我还没有想出如何将它写成一个函数，因为我似乎无法将重复值正确地存储在一个函数中。
当然，如果这产生了好的结果，我们希望 link 对 OEIS 参考文献进行讨论。
这是一个link序列： https://oeis.org/A077947
这与您的代码的行为相同，而且可以说更漂亮。您可能会看到使魔法常量不那么随意的方法。
factors = [ 1, 1, 2, 5, 9, 18 ]
cofactors = [ 36*(2**n) for n in range(6) ]
for x in range(10):
print(*factors)
for i in range(6):
factors[i] = factors[i] + cofactors[i] * 64**x
要仅计算其中一个子序列，在迭代时保持 i 固定就足够了。
这是一种无需第二个 for 循环的替代方法：
sequences = [ 1, 1, 2, 5, 9, 18 ]
multipliers = [ 36, 72, 144, 288, 576, 1152 ]
for x in range(100):
print(*sequences)
sequences = [ s + m*64**x for s,m in zip(sequences,multipliers) ]
[编辑] 查看这些值我注意到这个特定的序列也可以通过以下方式获得：
N[i+1] = 2 * N[i] + (-1,0,1 循环)
或
N[i+1] = 2 * N[i] + i mod 3 - 1 （假设索引从零开始）
i N[i] [-1,0,1] N[i+1]
0 1 -1 --> 2*1 - 1 --> 1
1 1 0 --> 2*1 + 0 --> 2
2 2 1 --> 2*2 + 1 --> 5
3 5 -1 --> 2*5 - 1 --> 9
4 9 0 --> 2*9 + 0 --> 18
...
因此，生成序列的更简单的循环可以是：
n = 1
for i in range(100):
print(n)
n = 2*n + i % 3 - 1
使用 functools 的 reduce 函数可以使这个更加简洁：
from functools import reduce
sequence = reduce(lambda s,i: s + [s[-1]*2 + i%3 - 1],range(20),[1])
print(sequence)
>>> [1, 1, 2, 5, 9, 18, 37, 73, 146, 293, 585, 1170, 2341, 4681, 9362, 18725, 37449, 74898, 149797, 299593, 599186]
使用您的多渠道方法和我建议的公式，这将得出：
sequences = [ 1, 1, 2, 5, 9, 18 ]
multipliers = [ 36, 72, 144, 288, 576, 1152 ]
allSequences = reduce(lambda ss,x: ss + [[ s + m*64**x for s,m in zip(ss[-1],multipliers) ]],range(100),[sequences])
for seq in allSequences: print(*seq) # print 6 by 6
[EDIT2] 如果您的所有序列都将具有相似的模式（即起始通道、乘数和计算公式），您可以将此类序列的打印概括为一个函数，因此每个序列只需要一行：
def printSeq(calcNext,sequence,multipliers,count):
for x in range(count):
print(*sequence)
sequence = [ calcNext(x,s,m) for s,m in zip(sequence,multipliers) ]
printSeq(lambda x,s,m:s*2+m*64**x,[1,1,2,5,9,18],multipliers=[36,72,144,288,576,1152],count=100)
[EDIT3] 改进 printSeq 函数。
我相信您不会总是需要一组乘数来计算每个通道中的下一个值。对该函数的改进是为 lambda 函数提供通道索引而不是乘数。这将允许您在需要时使用乘数数组，但也会让您使用更通用的计算。
def printSeq(name,count,calcNext,sequence):
p = len(sequence)
for x in range(count):
print(name, x,":","\t".join(str(s) for s in sequence))
sequence = [ calcNext(x,s,c,p) for c,s in enumerate(sequence) ]
lambda 函数有 4 个参数，预计 return 指定通道的下一个序列值：
s : current sequence value for the channel
x : iteration number
c : channel index (zero based)
p : number of channels
因此，在公式中使用数组可以这样表示：
printSeq("A077947",100,lambda x,s,c,p: s + [36,72,144,288,576,1152][c] * 64**x, [1,1,2,5,9,18])
但您也可以使用基于频道索引（和频道数）的更通用的公式：
printSeq("A077947",100,lambda x,s,c,p: s + 9 * 2**(p*x+c+2), [1,1,2,5,9,18])
或（基于 2*S + i%3 - 1 的 6 个通道）：
printSeq("A077947",100,lambda x,s,c,p: 64*s + 9*(c%3*2 - (c+2)%3 - 1) ,[1,1,2,5,9,18])
printSeq("A077947",100,lambda x,s,c,p: 64*s + 9*[-3,1,2][c%3],[1,1,2,5,9,18])
我的理由是，如果你有一个函数可以根据序列中的当前索引和值计算下一个值，你应该能够定义一个跨步函数来计算 N 个索引的值更远。
给定 F(i,S[i]) --> i+1,S[i+1]
F2(i,S[i]) --> i+2,S[i+2] = F(F(i,S[i]))
F3(i,S[i]) --> i+3,S[i+3] = F(F(F(i,S[i])))
...
F6(i,S[i]) --> i+6,S[i+6] = F(F(F(F(F(F(i,S[i]))))))
...
Fn(i,S[i]) --> i+n,S[i+n] = ...
这将始终有效，不需要乘法器数组。大多数情况下，应该可以仅使用代数来简化 Fn。
例如 A001045 : F(i,S) = i+1, 2*S + (-1)**i
printSeq("A001045",20,lambda x,s,c,p: 64*s + 21*(-1)**(x*p+c),[0,1,1,3,5,11])
请注意，从第 3 个值开始，可以在不知道索引的情况下计算该序列中的下一个值：
A001045: F(S) = 2*S + 1 - 2*0**((S+1)%4)
这是一个使用生成器的替代版本：
def Jacobsthal():
roots = [1, 1, 2, 5, 9, 18]
x = 0
while True:
yield roots
for i in range(6):
roots[i] += 36 * 2**i * 64**x
x += 1
这里有一个安全的使用方法：
j = Jacobsthal()
for _ in range(10):
print(j.send(None))
从OEIS注解可以看出，只需要3个初始值和1次递归3个前面的序列元素a(n-1)、a(n-2)和a(n-3)。然后您就可以轻松获得该系列。
# a(n-1)+a(n-2)+2*a(n-3)
# start values: 1, 1, 2
a1, a2, a3 = 1,1,2 #initial vales
m3, m2, m1 = 1,1,2 #multipliers for a(n-1):m3 a(n-2):m2 and a(n-3):m1
print a1, a2, a3,
for i in range(16):
a1,a2,a3=a2,a3,a3*m3+a2*m2+a1*m1
print a3,
这给你 1 1 2 5 9 18 37 73 146 293 585 1170 2341 4681 9362 18725 ...