在 C++ 中：

int main()
{
    for(int n = 1; n <= 50; n++)
        cout << setw(2) << n << "  "
            << setw(20) << ((n < 2) ? 0 : (unsigned long long)((1llu<<(n-2))/7))
            << endl;
}

据我了解，在 Python 中，

for n in range(1,1001):
    print(n,"  ",(1<<(n-2)) // 7)

https://oeis.org/A033138

适用于 n 到一千（Python 中为 15 秒）的方法。

基于基于三元组的交替总和的 7 整除的符号（此处）。我们把数字分成3位数的面，把它们加起来，如果结果能被7整除，那么也是原来的数字。

77744744: 77-744+744 = 77 делится на 7

对于初级三合会，只有两个选项，立即设置，我们整理第二个三合会。我们得到所有组合的可能数量的计数器，选择那些可以被 7 整除的数量。

from collections import Counter

def num74(n):
    if n < 5: return 0
    triads = [[444,447,474,477,744,747,774,777],[4,7],[44,47,74,77]]
    sums = Counter({444:1, 744:1})
    mul = 1
    for k in range(n//3-1):
        mul = - mul
        temp = Counter()
        for v in triads[0]:
            for ky in sums:
                temp[ky + v*mul] += sums[ky]
        sums = temp
    if n%3:
        mul = - mul
        temp = Counter()
        for v in triads[n%3]:
            for ky in sums:
                temp[ky + v*mul] += sums[ky]
        sums = temp

    res = sum([sums[ky] for ky in sums if (ky%7) == 0])
    return(res)

我没有@Harry 猜想的那么漂亮的解决方案，但是

• 它具有复杂性O(log n)，其中n是数字中的位数，
• 它是可解释的、可证明的和可推广的。

首先，请注意所有数字必须看起来像a = <b44>，即以两个四位结尾。否则，a它不能被 4 整除。

让我们找到余数a % 7：(b*100 + 44) % 7 = (b*100)%7 + 44%7 = (b%7)*(100%7) + 2 = (b%7)*2 + 2

因此，当且仅当它以 44 结尾时，a它可以被 7 和 4 整除，并且通过丢弃最后两位数字获得的数字必须能被 7 整除，余数为 6。为什么是 6？因为aba6*2 +2 = 14 == 0 (mod 7)

要找出有多少这样的数字很长，m = n-2您需要稍微概括一下这个问题。

m让我们知道由长度为:的 4 和 7 组成的数字的残差直方图r_m = [ r0_m, r1_m, ..., r6_m]，其中r0_m是长度m可被 7 整除且余数为 0 的数字数，r1_m- 余数为 1 等。

r_m又是如何关联的r_{m+1}？

设b一些长度m。让我们把右边的数字 4 和 7 添加到它上面。我们得到两个数字b1 = <b4>, b2 = <b7>。它们除以 7 的余数是多少？

b1 % 7 = ((10*b) + 4)%7 = (b%7)*(10%7) + 4 = 3*(b%7) + 4

b2 % 7 = ((10*b) + 7)%7 = (b%7)*(10%7) + 0 = 3*(b%7)

也就是说，每个余数r生成两个余数3r和3r+4。我们明确地写出来。左边是长度数的余数，m右边是长度数的两个余数m+1

0 -> 0, 4
1 -> 3, 0
2 -> 6, 3
3 -> 2, 6
4 -> 5, 2
5 -> 1, 5
6 -> 4, 1

在矩阵形式中，它看起来像这样：

因为m=0只有一个数，0。因此，初始残差向量r_0 = [1,0,0,0,0,0,0]

r_m = (A^m)@r_0

这是从到A的转移矩阵r_mr_{m+1}

要回答原始问题，您需要计算r_{n-2}并从中获取第六个分量（从 0 开始编号），对应于n-2除以 7 时余数为 6的长度数

import numpy as np

A = np.array([[1,1,0,0,0,0,0],
              [0,0,0,0,0,1,1],
              [0,0,0,1,1,0,0],
              [0,1,1,0,0,0,0],
              [1,0,0,0,0,0,1],
              [0,0,0,0,1,1,0],
              [0,0,1,1,0,0,0]
             ], dtype=np.object)

r0 = np.array([1,0,0,0,0,0,0])

for i in range(3, 101):
    B = np.linalg.matrix_power(A, i-2)
    r = B@r0
    print(i, r[6])

类型在构造函数A中明确指定，object因此元素是任意长度的 Pythonic 整数。如果未指定类型，则它将numpy优化整数并将其存储为 32 位整数。

对于求幂，使用numpy.linalg.matrix_power，因为内置函数按元素求幂pow。

结果

3 0
4 0
5 1
6 2
7 4
8 9
9 18
10 36
11 73
12 146
13 292
14 585
15 1170
16 2340
17 4681
18 9362
19 18724
20 37449
21 74898
22 149796
23 299593
24 599186
25 1198372
26 2396745
27 4793490
28 9586980
29 19173961
30 38347922
31 76695844
32 153391689
33 306783378
34 613566756
35 1227133513
36 2454267026
37 4908534052
38 9817068105
39 19634136210
40 39268272420
41 78536544841
42 157073089682
43 314146179364
44 628292358729
45 1256584717458
46 2513169434916
47 5026338869833
48 10052677739666
49 20105355479332
50 40210710958665
51 80421421917330
52 160842843834660
53 321685687669321
54 643371375338642
55 1286742750677284
56 2573485501354569
57 5146971002709138
58 10293942005418276
59 20587884010836553
60 41175768021673106
61 82351536043346212
62 164703072086692425
63 329406144173384850
64 658812288346769700
65 1317624576693539401
66 2635249153387078802
67 5270498306774157604
68 10540996613548315209
69 21081993227096630418
70 42163986454193260836
71 84327972908386521673
72 168655945816773043346
73 337311891633546086692
74 674623783267092173385
75 1349247566534184346770
76 2698495133068368693540
77 5396990266136737387081
78 10793980532273474774162
79 21587961064546949548324
80 43175922129093899096649
81 86351844258187798193298
82 172703688516375596386596
83 345407377032751192773193
84 690814754065502385546386
85 1381629508131004771092772
86 2763259016262009542185545
87 5526518032524019084371090
88 11053036065048038168742180
89 22106072130096076337484361
90 44212144260192152674968722
91 88424288520384305349937444
92 176848577040768610699874889
93 353697154081537221399749778
94 707394308163074442799499556
95 1414788616326148885598999113
96 2829577232652297771197998226
97 5659154465304595542395996452
98 11318308930609191084791992905
99 22636617861218382169583985810
100 45273235722436764339167971620

由于求幂m使用具有复杂度的算法O(log m)，所以 的复杂度n等于O(log(n-2))，即O(log n)

该算法也可用于其他条件。例如，对于由数字 5 和 9 组成的数字，它们可以被 5 和 9 整除。余数将有不同的矩阵和不同的条件，但解决方案将保持不变。

甚至可以不使用两个四人组进行优化。在这种情况下，矩阵是 28 x 28 。

A28 = np.zeros((28,28), dtype=np.object)
# подготовка матрицы переходов
for r in range(28):
    i1 = (r*10+4)%28
    i2 = (r*10+7)%28
    A28[i1, r] = 1
    A28[i2, r] = 1

# начальный вектор остатков
r028 = np.zeros(28, dtype=np.object)
r028[0] = 1

for i in range(1, 101):
    B = np.linalg.matrix_power(A28, i)
    r = B@r028
    print(i, r[0])

结果是一样的。

@Harry 的解决方案

我不会承诺推导出确切的公式，但2^(n-2)//7它很容易证明是正确的。加减一）

n只有合适的长度数字2^(n-2)，因为最后两位数字是固定的。除以 7 时的余数是均匀分布的，相差不超过 1。因此，具有固定余数的数字的数量——在我们的例子中，余数 6——是2^(n-2) / 7正负 1。

我不知道如何严格证明关于残差均匀分布的论文。从一般考虑，矩阵的结构说明了这一点A——每行正好有两个，行是混合的，所以矩阵A应该很好地对齐残差的直方图，混合最大值和最小值。

最令人惊讶的是，该公式不仅适用于大数，也适用于小数，其中波动加负1与函数的值相当相称

我会使用这个算法：

parts = [[444, 447, 474, 477, 744, 747, 774, 777], [4], [44, 47, 74, 77]]
ends = [[444, 744], [4], [44]]


def calculate(count, size, pos, value):
    res = 0

    if pos < count - 1:
        for part in parts[size]:
            res += calculate(count, 0, pos + 1, value + (part if pos % 2 else -part))
    else:
        for part in ends[size]:
            if (value + (part if pos % 2 else -part)) % 7 == 0:
                return 1

    return res


n = 10
print(calculate(n // 3 + (n % 3 != 0), n % 3, 0, 0))

1. 基本原则 - 搜索通过三倍数
2. 最后一个三元组是根据可被 4 整除的标准选择的 - （最后 2 位数字必须能被 4 整除）
3. 根据该数字的偶数和奇数三倍（千）之和必须能被 7 整除的标准，检查数字本身是否能被 7 整除