斯特恩双原子序列的高效生成

Question

Stern的双原子序列可以更详细地阅读over here；但是，为了我的目的，我现在将定义它。在

---

Stern双原子序列的定义

让n是生成fusc函数的数字。表示为fusc(n)。在

如果n为0，则返回值为0。
如果n为1，则返回值为1。在

如果n为偶数，则返回值为fusc(n / 2)。
如果n为奇数，则返回值为fusc((n - 1) / 2) + fusc((n + 1) / 2)。在

---

目前，我的Python代码在大部分的生成过程中都是蛮力的，除了除以两部分之外，因为它永远不会产生任何变化。在

def fusc (n):
    if n <= 1:
        return n

    while n > 2 and n % 2 == 0:
        n /= 2

    return fusc((n - 1) / 2) + fusc((n + 1) / 2)

但是，我的代码必须能够处理大小为1000s百万位的数字，并且递归地运行函数数千数百万次似乎不是非常有效或实用的。在

有没有什么方法可以在算法上改进我的代码，使大量的数字可以通过而不必递归调用函数那么多次？在

Answer 1

对于一百万位的记忆，递归堆栈将非常大。我们可以先看看一个足够大的数字，我们可以用手来计算，fusc(71)在这种情况下：

fusc(71) = fusc(35) + fusc(36)

fusc(35) = fusc(17) + fusc(18)
fusc(36) = fusc(18)

fusc（71）=1*fusc（17）+2*fusc（18）

fusc(17) = fusc(8) + fusc(9)
fusc(18) = fusc(9)

fusc（71）=1*fusc（8）+3*fusc（9）

fusc(8) = fusc(4)
fusc(9) = fusc(4) + fusc(5)

fusc（71）=4*fusc（4）+3*fusc（5）

fusc(4) = fusc(2)
fusc(3) = fusc(1) + fusc(2)

fusc（71）=7*fusc（2）+3*fusc（3）

fusc(2) = fusc(1)
fusc(3) = fusc(1) + fusc(2)

fusc（71）=11*fusc（1）+3*fusc（2）

fusc(2) = fusc(1)

fusc（71）=14*fusc（1）=14

我们意识到在这种情况下，我们可以完全避免递归，因为我们总是可以用a * fusc(m) + b * fusc(m+1)的形式来表达fusc(n)，同时将m的值减少到0。从上面的示例中，您可能会发现以下模式：

if m is odd:
a * fusc(m) + b * fusc(m+1) = a * fusc((m-1)/2) + (b+a) * fusc((m+1)/2)
if m is even:
a * fusc(m) + b * fusc(m+1) = (a+b) * fusc(m/2) + b * fusc((m/2)+1)

因此，您可以使用一个简单的循环函数在O（lg（n））时间内解决问题

def fusc(n):
  if n == 0: return 0
  a = 1
  b = 0
  while n > 0:
    if n%2:
      b = b + a
      n = (n-1)/2
    else:
      a = a + b
      n = n/2
  return b

Answer 2

lru_cache对你的情况很有帮助。确保maxsize是2的幂次方。你的应用程序可能需要调整一下这个尺寸。cache_info()会有帮助的。在

也可以使用//代替/进行整数除法。在

from functools import lru_cache


@lru_cache(maxsize=512, typed=False)
def fusc(n):

    if n <= 1:
        return n

    while n > 2 and n % 2 == 0:
        n //= 2

    return fusc((n - 1) // 2) + fusc((n + 1) // 2)

print(fusc(1000000000078093254329870980000043298))
print(fusc.cache_info())

是的，这只是菲利普·马尔查克提出的一个建议。在

您可以在while循环中使用位操作获得额外的tiny加速：

^{pr2}$

---

更新：

这里有一个简单的“手工”方法：

def fusc(n, _mem={}):  # _mem will be the cache of the values 
                       # that have been calculated before
    if n in _mem:      # if we know that one: just return the value
        return _mem[n]

    if n <= 1:
        return n

    while not n & 1:
        n >>= 1
    if n == 1:
        return 1    

    ret = fusc((n - 1) // 2) + fusc((n + 1) // 2)
    _mem[n] = ret  # store the value for next time
    return ret

---

更新

在阅读了short article by dijkstra himself的小更新之后。在

文章指出，f(n) = f(m)如果m的第一位和最后一位与{}相同，并且两者之间的位是颠倒的。这样做的目的是使n越小越好。在

这就是位掩码(1<<n.bit_length()-1)-2的用途（第一位和最后一位是0；中间的那些1；xoringn，如上所述）。在

我只能做一些小的基准测试；我感兴趣的是，如果这对您的输入的魔力有任何帮助。。。这将减少缓存的内存，并有望带来一些加速。在

def fusc_ed(n, _mem={}):

    if n <= 1:
        return n

    while not n & 1:
        n >>= 1
    if n == 1:
        return 1

    # https://www.cs.utexas.edu/users/EWD/transcriptions/EWD05xx/EWD578.html
    # bit invert the middle bits and check if this is smaller than n
    m = n ^ (1<<n.bit_length()-1)-2
    n = m if m < n else n

    if n in _mem:
        return _mem[n]

    ret = fusc(n >> 1) + fusc((n >> 1) + 1)
    _mem[n] = ret
    return ret

我不得不增加递归限制：

import sys
sys.setrecursionlimit(10000)  # default limit was 1000

基准测试给出了奇怪的结果；使用下面的代码并确保我总是启动一个新的interperter（有一个空的_mem），我有时会得到更好的运行时；在其他情况下，新代码会慢一些。。。在

基准代码：

print(n.bit_length())

ti = timeit('fusc(n)', setup='from __main__ import fusc, n', number=1)
print(ti)

ti = timeit('fusc_ed(n)', setup='from __main__ import fusc_ed, n', number=1)
print(ti)

我得到了三个随机结果：

6959
24.117448464001427
0.013900151001507766

6989
23.92404893300045
0.013844672999766772

7038
24.33894686200074
24.685758719999285

那就是我停下来的地方。。。在