Skip to content

Instantly share code, notes, and snippets.

@taroyabuki

taroyabuki/magicsquare.cpp Secret

Last active December 20, 2018 15:21
Show Gist options
  • Save taroyabuki/b94ed5b24a385e460af7 to your computer and use it in GitHub Desktop.
Save taroyabuki/b94ed5b24a385e460af7 to your computer and use it in GitHub Desktop.
Download ZIP
スパコンで約2時間36分かかったという、5×5の魔方陣の全解列挙を、パソコンで試す（C++）http://blog.unfindable.net/archives/7179
Raw
magicsquare.cpp
#include <iostream>
#include <fstream>
#include <sstream>
#include <algorithm>
#include <ctime>
#include <omp.h>
using namespace std;
const char* outputfile = "result";
int main()
{
clock_t start = clock();
const int B = 5;
const int N = 25;
const int S = 65;
const int mask = (1 << B) - 1;
int threads = omp_get_max_threads();
cout << "num of threads: " << threads << endl;
int* counts = new int[threads];
ofstream** files = new ofstream*[threads];
for (int t = 0; t < threads; ++t){
counts[t] = 0;
stringstream ss;
ss << outputfile << t;
files[t] = new ofstream(ss.str().c_str(), ios_base::binary);
}
#pragma omp parallel
#pragma omp for schedule(dynamic, 1024)
for (int i = 0; i < (1 << 25); ++i) {
if ((i % (1 << 20)) == 0) cout << (i / (1 << 20)) << "/32" << endl;
const int x7 = ((i >> (B * 0)) & mask) + 1;
if (x7 > N - 3) continue;
const int x9 = x7 + ((i >> (B * 1)) & mask) + 1;
if (x9 > N) continue;
const int x17 = x9 + ((i >> (B * 2)) & mask) + 1;
if (x17 > N) continue;
int used = (1 << x7) + (1 << x9) + (1 << x17);
const int x19 = x7 + ((i >> (B * 3)) & mask) + 1;
if (x19 > N || (used & (1 << x19))) continue;
used ^= (1 << x19);
const int x13 = ((i >> (B * 4)) & mask) + 1;
if (x13 > N || (used & (1 << x13))) continue;
used ^= (1 << x13);
for (int x1 = max(1, S - x7 - x13 - x19 - N); x1 <= min(N, S - x7 - x13 - x19 - 1); ++x1) {
if (!(used & (1 << x1))) {
used ^= (1 << x1);
const int x25 = S - x1 - x7 - x13 - x19;
if (!(used & (1 << x25))) {
used ^= (1 << x25);
for (int x5 = max(1, S - x9 - x13 - x17 - N); x5 <= min(N, S - x9 - x13 - x17 - 1); ++x5) {
if (!(used & (1 << x5))) {
used ^= (1 << x5);
const int x21 = S - x5 - x9 - x13 - x17;
if (!(used & (1 << x21))) {
used ^= (1 << x21);
for (int x2 = max(1, max(S - x1 - x5, S - x7 - x17) - 49); x2 <= min(N, min(S - x1 - x5, S - x7 - x17) - 3); ++x2) {
if (!(used & (1 << x2))) {
used ^= (1 << x2);
for (int x3 = max(1, S - x1 - x2 - x5 - N); x3 <= min(N, S - x1 - x2 - x5 - 1); ++x3) {
if (!(used & (1 << x3))) {
used ^= (1 << x3);
const int x4 = S - x1 - x2 - x3 - x5;
if (!(used & (1 << x4))) {
used ^= (1 << x4);
for (int x12 = max(1, S - x2 - x7 - x17 - N); x12 <= min(N, S - x2 - x7 - x17 - 1); ++x12) {
if (!(used & (1 << x12))) {
used ^= (1 << x12);
const int x22 = S - x2 - x7 - x12 - x17;
if (!(used & (1 << x22))) {
used ^= (1 << x22);
for (int x24 = max(1, max(S - x4 - x9 - x19, S - x21 - x22 - x25) - N); x24 <= min(N, min(S - x4 - x9 - x19, S - x21 - x22 - x25) - 1); ++x24) {
if (!(used & (1 << x24))) {
used ^= (1 << x24);
const int x14 = S - x4 - x9 - x19 - x24;
if (1 <= x14 && x14 <= N && !(used & (1 << x14))) {
used ^= (1 << x14);
const int x23 = S - x21 - x22 - x24 - x25;
if (1 <= x23 && x23 <= N && !(used & (1 << x23))) {
used ^= (1 << x23);
for (int x8 = max(1, S - x3 - x13 - x23 - N); x8 <= min(N, S - x3 - x13 - x23 - 1); ++x8) {
if (!(used & (1 << x8))) {
used ^= (1 << x8);
const int x18 = S - x3 - x8 - x13 - x23;
if (!(used & (1 << x18))) {
used ^= (1 << x18);
for (int x6 = max(1, S - x7 - x8 - x9 - N); x6 <= min(N, S - x7 - x8 - x9 - 1); ++x6) {
if (!(used & (1 << x6))) {
used ^= (1 << x6);
const int x10 = S - x6 - x7 - x8 - x9;
if (!(used & (1 << x10))) {
used ^= (1 << x10);
for (int x11 = max(1, max(S - x1 - x6 - x21, S - x12 - x13 - x14) - N); x11 <= min(N, min(S - x1 - x6 - x21, S - x12 - x13 - x14) - 1); ++x11) {
if (!(used & (1 << x11))) {
used ^= (1 << x11);
const int x15 = S - x11 - x12 - x13 - x14;
if (1 <= x15 && x15 <= N && !(used & (1 << x15))) {
used ^= (1 << x15);
const int x16 = S - x1 - x6 - x11 - x21;
if (1 <= x16 && x16 <= N && !(used & (1 << x16))) {
used ^= (1 << x16);
const int x20 = S - x5 - x10 - x15 - x25;
if (S == x16 + x17 + x18 + x19 + x20 && 1 <= x20 && x20 <= N && !(used & (1 << x20))) {
int t = omp_get_thread_num();
++counts[t];
*files[t] << (char)x1 << (char)x2 << (char)x3 << (char)x4 << (char)x5 << (char)x6 << (char)x7 << (char)x8 << (char)x9 << (char)x10 << (char)x11 << (char)x12 << (char)x13 << (char)x14 << (char)x15 << (char)x16 << (char)x17 << (char)x18 << (char)x19 << (char)x20 << (char)x21 << (char)x22 << (char)x23 << (char)x24 << (char)x25;
}
used ^= (1 << x16);
}
used ^= (1 << x15);
}
used ^= (1 << x11);
}
}
used ^= (1 << x10);
}
used ^= (1 << x6);
}
}
used ^= (1 << x18);
}
used ^= (1 << x8);
}
}
used ^= (1 << x23);
}
used ^= (1 << x14);
}
used ^= (1 << x24);
}
}
used ^= (1 << x22);
}
used ^= (1 << x12);
}
}
used ^= (1 << x4);
}
used ^= (1 << x3);
}
}
used ^= (1 << x2);
}
}
used ^= (1 << x21);
}
used ^= (1 << x5);
}
}
used ^= (1 << x25);
}
used ^= (1 << x1);
}
}
}
int sum = 0;
for (int t = 0; t < threads; ++t) {
sum += counts[t];
files[t]->close();
}
cout << "num of solution: " << sum << endl;
clock_t now = clock();
double duration = double(now - start) / CLOCKS_PER_SEC;
cout << "time for count: " << duration << " sec." << endl;
return 0;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment