位运算优化素数环爆搜

前言

程序正在服务器上运行（我的服务器是出了名的腾讯云老爷机），估计没有个几年是跑不完了，不过我会逐步更新计算结果。

感谢HJQ巨佬提供的OEIS链接

OEIS 上给出的 $N/2\le 15$ 时的答案序列

当前计算结果

n	cnt	用时(ms)
4	1	0.005
6	1	0.001
8	2	0.001
10	48	0.009
12	512	0.075
14	1440	0.577
16	40512	10.085
18	385072	114.390
20	3154650	1523.373
22	106906168	45640.766
24	3197817022	1472444.354（24.53 分钟）
26	82924866213	49227017.707 (13.67 小时)
28	4025168862425	2237223393.049 (25.89 天)
30	Unknown	预计需要 858 天

代码

/// 采用位运算技巧来优化素数环问题的求解
#pragma GCC optimize(2) 

#include <cstdio>
#include <ctime>
#include <bitset>
using namespace std;

int isprime(int x) { /// 判断 x 是否是质数 
	if(x <= 1) return false;
	if(x == 2) return true;
	for(int i = 2; i*i <= x; i ++) {
		if(x % i == 0) return false;
	}
	return true; /// x 是质数 
}

long long cnt = 0; /// 统计答案总数 

long long avai[60 + 1];
int lg[(1 << 22) + 1]; /// 可用后继集合 

int Lg(long long x) { /// 计算对数 
	int ans = 0;
	while(x >= (1 << 20)) {
		ans += 20;
		x /= (1<<20);
	}
	return ans + lg[x];
}

void init(int n) {
	for(int i = 1; i <= n; i ++) {
		for(int j = 1; j <= n; j ++) {
			if(i == j) continue;
			if(isprime(i + j)) {
				avai[i] |= (1ll << (j-1)); /// avai[i] 表示能放在 i 后面的数的集合 
			}
		}
	}
	for(int j = 1, i = 1; i <= 20; i ++, j <<= 1) { /// 预处理对数数组 
		lg[j] = i;
	}
}

inline long long lowbit(long long x) { /// 得到最低二进制位 
	return x & (-x);
}

inline int inset(int x, long long S) { /// 判断 x 是否在集合 S 中 
	return ((1 << (x-1)) & S) != 0;
}

int arr[60 + 1], n;
void dfs(int pos, long long res) {
	if(pos == n + 1) {
		if(inset(arr[n], avai[1])) cnt ++;
	}else {
		long long S = res & avai[arr[pos - 1]]; /// 当前所有可用位置 
		while(S) {
			long long tmp = lowbit(S);
			int id = Lg(tmp);
			arr[pos] = id;
			dfs(pos + 1, res ^ tmp); /// 递归填数 
			S ^= tmp;
		}
	}
}

int main() {
	FILE* fpout = fopen("prime_ring.txt", "w"); /// 清空文件 
	fclose(fpout);
	init(60);
	for(n = 4; n <= 60; n += 2) {
		cnt = 0;
		FILE* fpout = fopen("prime_ring.txt", "a"); /// 文件尾部追加 
		int begin = clock(); /// 计时开始 
		if(!(n & 1)) {
			arr[1] = 1;
			dfs(2, (1 << n)-2);
		}
		fprintf(fpout, "n = %d, find %lld solutions.\n", n, cnt / 2); /// 考虑顺时针和逆时针，更正 %d -> %lld! 
		fprintf(fpout, "time consumed: %.3lf sec.\n\n", (clock() - begin)/(CLOCKS_PER_SEC * 1.0)); /// 输出计时结果 
		fflush(fpout);
		fclose(fpout); /// 之前版本的程序写错了，导致我在 linux 服务器上单位把毫秒标成了秒 
	}
	return 0;
}