Create a gist now

Instantly share code, notes, and snippets.

@taroyabuki /magicsquare.cpp Secret
Last active Sep 6, 2016

What would you like to do?
Download ZIP
スパコンで約2時間36分かかったという、5×5の魔方陣の全解列挙を、パソコンで試す（C++）http://blog.unfindable.net/archives/7179
Raw
magicsquare.cpp
#include <iostream>
#include <fstream>
#include <string>
#include <algorithm>
#include <chrono>
#include <vector>
#include <omp.h>
#define MAXMIN(i) int ux##i, un##i; setUnusedMaxMin(used, ux##i, un##i)
using namespace std;
string outputfile = "result";
const int n = 5;
const int N = n * n;
const int S = N * (N + 1) / 2 / n;
inline void setUnusedMaxMin(const unsigned long used, int& ux, int& un) {
unsigned long tmp = used ^ ((1 << (N + 1)) - 2);
#ifdef _MSC_VER
_BitScanReverse((unsigned long*)&ux, tmp);
_BitScanForward((unsigned long*)&un, tmp);
#endif
#ifdef __GNUC__
ux = 31 - __builtin_clz(tmp);
un = __builtin_ctz(tmp);
#endif
}
int main() {
auto t0 = std::chrono::high_resolution_clock::now();
int threads = omp_get_max_threads();
cout << "num of threads: " << threads << endl;
auto files = vector<ofstream>();
for (int t = 0; t < threads; ++t) {
string filename = outputfile + to_string(t);
files.emplace_back(filename, ios_base::binary);
}
int sum = 0;
cout << "x7 = 1" << endl;
#pragma omp parallel for schedule(dynamic, 1024) reduction (+:sum)
for (int i = 0 * N * N * N * N + 1 * N * N * N + 2 * N * N + 3 * N + 4; i <= (N - 4) * N * N * N * N + (N - 1) * N * N * N + (N - 2) * N * N + (N - 3) * N + (N - 4); ++i) {
const int x07 = i / (N * N * N * N) + 1; if ((i % (N * N * N * N)) == 0) cout << "x7 = " << x07 << endl; unsigned long used = (1 << x07);
const int x09 = (i / (N * N * N)) % N + 1; if (x09 > x07 && x09 < N) { used += (1 << x09);
const int x17 = (i / (N * N)) % N + 1; if (x17 > x09) { used += (1 << x17);
const int x19 = (i / N) % N + 1; if (x19 > x07 && !(used & (1 << x19))) { used ^= (1 << x19);
const int x13 = i % N + 1; if (x13 <= 13 && !(used & (1 << x13))) { used ^= (1 << x13); MAXMIN(01);
for (int x01 = max(un01, S - x07 - x13 - x19 - ux01); x01 <= min(ux01, S - x07 - x13 - x19 - un01); ++x01) { if (!(used & (1 << x01))) { used ^= (1 << x01);
const int x25 = S - x01 - x07 - x13 - x19; if (!(used & (1 << x25))) { used ^= (1 << x25); MAXMIN(05);
for (int x05 = max(un05, S - x09 - x13 - x17 - ux05); x05 <= min(ux05, S - x09 - x13 - x17 - un05); ++x05) { if (!(used & (1 << x05))) { used ^= (1 << x05);
const int x21 = S - x05 - x09 - x13 - x17; if (!(used & (1 << x21))) { used ^= (1 << x21); MAXMIN(02);
for (int x02 = max(un02, S - min(x01 + x05, x07 + x17) - ux02 - ux02 - 1); x02 <= min(ux02, S - max(x01 + x05, x07 + x17) - un02 - un02 - 1); ++x02) { if (!(used & (1 << x02))) { used ^= (1 << x02); MAXMIN(03);
for (int x03 = max(un03, S - x01 - x02 - x05 - ux03); x03 <= min(ux03, S - x01 - x02 - x05 - un03); ++x03) { if (!(used & (1 << x03))) { used ^= (1 << x03);
const int x04 = S - x01 - x02 - x03 - x05; if (!(used & (1 << x04))) { used ^= (1 << x04); MAXMIN(12);
for (int x12 = max(un12, S - x02 - x07 - x17 - ux12); x12 <= min(ux12, S - x02 - x07 - x17 - un12); ++x12) { if (!(used & (1 << x12))) { used ^= (1 << x12);
const int x22 = S - x02 - x07 - x12 - x17; if (!(used & (1 << x22))) { used ^= (1 << x22); MAXMIN(24);
for (int x24 = max(un24, S - min(x04 + x09 + x19, x21 + x22 + x25) - ux24); x24 <= min(ux24, S - max(x04 + x09 + x19, x21 + x22 + x25) - un24); ++x24) { if (!(used & (1 << x24))) { used ^= (1 << x24);
const int x14 = S - x04 - x09 - x19 - x24; if (!(used & (1 << x14))) { used ^= (1 << x14);
const int x23 = S - x21 - x22 - x24 - x25; if (!(used & (1 << x23))) { used ^= (1 << x23); MAXMIN(08);
for (int x08 = max(un08, S - x03 - x13 - x23 - ux08); x08 <= min(ux08, S - x03 - x13 - x23 - un08); ++x08) { if (!(used & (1 << x08))) { used ^= (1 << x08);
const int x18 = S - x03 - x08 - x13 - x23; if (!(used & (1 << x18))) { used ^= (1 << x18); MAXMIN(06);
for (int x06 = max(un06, S - x07 - x08 - x09 - ux06); x06 <= min(ux06, S - x07 - x08 - x09 - un06); ++x06) { if (!(used & (1 << x06))) { used ^= (1 << x06);
const int x10 = S - x06 - x07 - x08 - x09; if (!(used & (1 << x10))) { used ^= (1 << x10); MAXMIN(11);
for (int x11 = max(un11, S - min(x01 + x06 + x21, x12 + x13 + x14) - ux11); x11 <= min(ux11, S - max(x01 + x06 + x21, x12 + x13 + x14) - un11); ++x11) { if (!(used & (1 << x11))) { used ^= (1 << x11);
const int x15 = S - x11 - x12 - x13 - x14; if (!(used & (1 << x15))) { used ^= (1 << x15);
const int x16 = S - x01 - x06 - x11 - x21; if (!(used & (1 << x16))) { used ^= (1 << x16);
const int x20 = S - x05 - x10 - x15 - x25; if (S == x16 + x17 + x18 + x19 + x20 && !(used & (1 << x20))) {
int t = omp_get_thread_num();
files[t] << (char)x01 << (char)x02 << (char)x03 << (char)x04 << (char)x05
<< (char)x06 << (char)x07 << (char)x08 << (char)x09 << (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;
if (x13 == 13) ++sum;
else {
sum += 2;
files[t] << (char)(26 - x01) << (char)(26 - x02) << (char)(26 - x03) << (char)(26 - x04) << (char)(26 - x05)
<< (char)(26 - x06) << (char)(26 - x07) << (char)(26 - x08) << (char)(26 - x09) << (char)(26 - x10)
<< (char)(26 - x11) << (char)(26 - x12) << (char)(26 - x13) << (char)(26 - x14) << (char)(26 - x15)
<< (char)(26 - x16) << (char)(26 - x17) << (char)(26 - x18) << (char)(26 - x19) << (char)(26 - x20)
<< (char)(26 - x21) << (char)(26 - x22) << (char)(26 - x23) << (char)(26 - x24) << (char)(26 - x25);
}}
used ^= (1 << x16); }
used ^= (1 << x15); }
used ^= (1 << x11); }}
used ^= (1 << x10); }
used ^= (1 << x06); }}
used ^= (1 << x18); }
used ^= (1 << x08); }}
used ^= (1 << x23); }
used ^= (1 << x14); }
used ^= (1 << x24); }}
used ^= (1 << x22); }
used ^= (1 << x12); }}
used ^= (1 << x04); }
used ^= (1 << x03); }}
used ^= (1 << x02); }}
used ^= (1 << x21); }
used ^= (1 << x05); }}
used ^= (1 << x25); }
used ^= (1 << x01); }}}}}}}
for (auto& f : files) f.close();
cout << "num of solution: " << sum << endl;
auto t1 = std::chrono::high_resolution_clock::now();
cout << std::chrono::duration_cast<std::chrono::milliseconds>(t1 - t0).count() / 1000. << " sec.\n";
return 0;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment