我有一些简单的代码可以执行以下操作。
它遍历所有可能长度为 n 的列表 F
与 +-1 个条目。对于每一个,它遍历所有可能的长度 2n
列表 S
与 +-1 个条目,其中 $S$ 的前半部分只是后半部分的副本。该代码计算 F
与长度为 n
的每个 S
子列表的内积。对于每个 F,S,它计算在第一个非零内积之前为零的内积。
这里是代码。
#!/usr/bin/python
from __future__ import division
import itertools
import operator
import math
n=14
m=n+1
def innerproduct(A, B):
assert (len(A) == len(B))
s = 0
for k in xrange(0,n):
s+=A[k]*B[k]
return s
leadingzerocounts = [0]*m
for S in itertools.product([-1,1], repeat = n):
S1 = S + S
for F in itertools.product([-1,1], repeat = n):
i = 0
while (i<m):
ip = innerproduct(F, S1[i:i+n])
if (ip == 0):
leadingzerocounts[i] +=1
i+=1
else:
break
print leadingzerocounts
n=14
的正确输出是
[56229888, 23557248, 9903104, 4160640, 1758240, 755392, 344800, 172320, 101312, 75776, 65696, 61216, 59200, 59200, 59200]
使用 pypy,对于 n = 14,这需要 1 分 18 秒。不幸的是,我真的很想在 16、18、20、22、24、26 上运行它。我不介意使用 numba 或 cython,但如果可能的话,我想靠近 python。
非常感谢任何帮助加快这一进程。
我会在这里记录最快的解决方案。 (如果我错过了更新的答案,请告诉我。)
- n = 22,时间为 9 分 35.081 秒,由 Eisenstat (C) 提供
- Eisenstat (pypy) 在 1 分 16.344 秒时的 n = 18
- n = 18,2 分 54.998 秒,Tupteq (pypy)
- n = 14 at 26s by Neil (numpy)
- n - 14 点 11 分 59 秒 192 秒 by kslote1 (pypy)
最佳答案
这个新代码利用问题的循环对称性获得了另一个数量级的加速。这个 Python 版本使用 Duval 算法枚举项链; C 版本使用蛮力。两者都包含下面描述的加速。 在我的机器上,C 版本在 100 秒内解决了 n = 20! 粗略的计算表明,如果你让它在单核上运行一周,它可以做到 n = 26,并且如下所述,它适合并行性。
import itertools
def necklaces_with_multiplicity(n):
assert isinstance(n, int)
assert n > 0
w = [1] * n
i = 1
while True:
if n % i == 0:
s = sum(w)
if s > 0:
yield (tuple(w), i * 2)
elif s == 0:
yield (tuple(w), i)
i = n - 1
while w[i] == -1:
if i == 0:
return
i -= 1
w[i] = -1
i += 1
for j in range(n - i):
w[i + j] = w[j]
def leading_zero_counts(n):
assert isinstance(n, int)
assert n > 0
assert n % 2 == 0
counts = [0] * n
necklaces = list(necklaces_with_multiplicity(n))
for combo in itertools.combinations(range(n - 1), n // 2):
for v, multiplicity in necklaces:
w = list(v)
for j in combo:
w[j] *= -1
for i in range(n):
counts[i] += multiplicity * 2
product = 0
for j in range(n):
product += v[j - (i + 1)] * w[j]
if product != 0:
break
return counts
if __name__ == '__main__':
print(leading_zero_counts(12))
C 版:
#include <stdio.h>
enum {
N = 14
};
struct Necklace {
unsigned int v;
int multiplicity;
};
static struct Necklace g_necklace[1 << (N - 1)];
static int g_necklace_count;
static void initialize_necklace(void) {
g_necklace_count = 0;
for (unsigned int v = 0; v < (1U << (N - 1)); v++) {
int multiplicity;
unsigned int w = v;
for (multiplicity = 2; multiplicity < 2 * N; multiplicity += 2) {
w = ((w & 1) << (N - 1)) | (w >> 1);
unsigned int x = w ^ ((1U << N) - 1);
if (w < v || x < v) goto nope;
if (w == v || x == v) break;
}
g_necklace[g_necklace_count].v = v;
g_necklace[g_necklace_count].multiplicity = multiplicity;
g_necklace_count++;
nope:
;
}
}
int main(void) {
initialize_necklace();
long long leading_zero_count[N + 1];
for (int i = 0; i < N + 1; i++) leading_zero_count[i] = 0;
for (unsigned int v_xor_w = 0; v_xor_w < (1U << (N - 1)); v_xor_w++) {
if (__builtin_popcount(v_xor_w) != N / 2) continue;
for (int k = 0; k < g_necklace_count; k++) {
unsigned int v = g_necklace[k].v;
unsigned int w = v ^ v_xor_w;
for (int i = 0; i < N + 1; i++) {
leading_zero_count[i] += g_necklace[k].multiplicity;
w = ((w & 1) << (N - 1)) | (w >> 1);
if (__builtin_popcount(v ^ w) != N / 2) break;
}
}
}
for (int i = 0; i < N + 1; i++) {
printf(" %lld", 2 * leading_zero_count[i]);
}
putchar('\n');
return 0;
}
您可以通过利用符号对称性 (4x) 并仅迭代那些通过第一个内积测试的向量(渐近地,O(sqrt(n))x)来获得一点加速。
import itertools
n = 10
m = n + 1
def innerproduct(A, B):
s = 0
for k in range(n):
s += A[k] * B[k]
return s
leadingzerocounts = [0] * m
for S in itertools.product([-1, 1], repeat=n - 1):
S1 = S + (1,)
S1S1 = S1 * 2
for C in itertools.combinations(range(n - 1), n // 2):
F = list(S1)
for i in C:
F[i] *= -1
leadingzerocounts[0] += 4
for i in range(1, m):
if innerproduct(F, S1S1[i:i + n]):
break
leadingzerocounts[i] += 4
print(leadingzerocounts)
C 版本,以了解我们在 PyPy 中损失了多少性能(PyPy 的 16 大致相当于 C 的 18):
#include <stdio.h>
enum {
HALFN = 9,
N = 2 * HALFN
};
int main(void) {
long long lzc[N + 1];
for (int i = 0; i < N + 1; i++) lzc[i] = 0;
unsigned int xor = 1 << (N - 1);
while (xor-- > 0) {
if (__builtin_popcount(xor) != HALFN) continue;
unsigned int s = 1 << (N - 1);
while (s-- > 0) {
lzc[0]++;
unsigned int f = xor ^ s;
for (int i = 1; i < N + 1; i++) {
f = ((f & 1) << (N - 1)) | (f >> 1);
if (__builtin_popcount(f ^ s) != HALFN) break;
lzc[i]++;
}
}
}
for (int i = 0; i < N + 1; i++) printf(" %lld", 4 * lzc[i]);
putchar('\n');
return 0;
}
这个算法是令人尴尬的并行,因为它只是对 xor
的所有值进行累加。对于 C 版本,粗略计算表明,几千小时的 CPU 时间足以计算出 n = 26
,按目前的费率计算,这相当于几百美元EC2。毫无疑问,需要进行一些优化(例如矢量化),但对于这样的一次性优化,我不确定程序员需要付出多少努力。
关于python - 如何在python中加速多个内积,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/26960749/