Create a gist now

Instantly share code, notes, and snippets.

@taroyabuki /magicsquare.cpp Secret
Last active Mar 6, 2018

What would you like to do?
スパコンで約2時間36分かかったという、5×5の魔方陣の全解列挙を、パソコンで試す(C++)http://blog.unfindable.net/archives/7179
#include <iostream>
#include <fstream>
#include <string>
#include <algorithm>
#include <chrono>
#include <vector>
#include <omp.h>
#ifdef _MSC_VER
#include <intrin.h>
#endif
#define MAXMIN(i) int ux##i, un##i; setUnusedMaxMin(used, ux##i, un##i)
#define FLIP(x) used ^= (1 << (x))
#define NOTUSED(x) !(used & (1 << (x)))
#define B(x) (char)(x)
#define C(x) (char)(26 - (x))
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; i <= 21 * N * N * N * N; ++i) {
const int x07 = i / (N * N * N * N) + 2; 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 (NOTUSED(x09) && x09 < N) { FLIP(x09);
const int x17 = (i / (N * N)) % N + 1; if (NOTUSED(x17)) { FLIP(x17);
const int x19 = (i / N) % N + 1; if (NOTUSED(x19)) { FLIP(x19);
const int x13 = i % N + 1; if (NOTUSED(x13) && x13 <= 13) { FLIP(x13); MAXMIN(01);
for (int x01 = max(un01, S - x07 - x13 - x19 - ux01); x01 <= min({ux01, x07, 22, S - x07 - x13 - x19 - un01}); ++x01) { if (NOTUSED(x01)) { FLIP(x01);
const int x25 = S - x01 - x07 - x13 - x19; if (NOTUSED(x25) && (x07 < x19 || x07 < x25)) { FLIP(x25); MAXMIN(05);
for (int x05 = max(un05, S - x09 - x13 - x17 - ux05); x05 <= min({ux05, 24, S - x09 - x13 - x17 - un05}); ++x05) { if (NOTUSED(x05) && (x07 < x09 || x07 < x05)) { FLIP(x05);
const int x21 = S - x05 - x09 - x13 - x17; if (NOTUSED(x21) && max(x05, x09) < max(x17, x21)) { FLIP(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 (NOTUSED(x02)) { FLIP(x02); MAXMIN(03);
for (int x03 = max(un03, S - x01 - x02 - x05 - ux03); x03 <= min(ux03, S - x01 - x02 - x05 - un03); ++x03) { if (NOTUSED(x03)) { FLIP(x03);
const int x04 = S - x01 - x02 - x03 - x05; if (NOTUSED(x04)) { FLIP(x04); MAXMIN(12);
for (int x12 = max(un12, S - x02 - x07 - x17 - ux12); x12 <= min(ux12, S - x02 - x07 - x17 - un12); ++x12) { if (NOTUSED(x12)) { FLIP(x12);
const int x22 = S - x02 - x07 - x12 - x17; if (NOTUSED(x22)) { FLIP(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 (NOTUSED(x24)) { FLIP(x24);
const int x14 = S - x04 - x09 - x19 - x24; if (NOTUSED(x14)) { FLIP(x14);
const int x23 = S - x21 - x22 - x24 - x25; if (NOTUSED(x23)) { FLIP(x23); MAXMIN(08);
for (int x08 = max(un08, S - x03 - x13 - x23 - ux08); x08 <= min(ux08, S - x03 - x13 - x23 - un08); ++x08) { if (NOTUSED(x08)) { FLIP(x08);
const int x18 = S - x03 - x08 - x13 - x23; if (NOTUSED(x18)) { FLIP(x18); MAXMIN(06);
for (int x06 = max(un06, S - x07 - x08 - x09 - ux06); x06 <= min(ux06, S - x07 - x08 - x09 - un06); ++x06) { if (NOTUSED(x06)) { FLIP(x06);
const int x10 = S - x06 - x07 - x08 - x09; if (NOTUSED(x10)) { FLIP(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 (NOTUSED(x11)) { FLIP(x11);
const int x15 = S - x11 - x12 - x13 - x14; if (NOTUSED(x15)) { FLIP(x15);
const int x16 = S - x01 - x06 - x11 - x21; if (NOTUSED(x16)) { FLIP(x16);
const int x20 = S - x05 - x10 - x15 - x25; if (NOTUSED(x20) && S == x16 + x17 + x18 + x19 + x20) {
int t = omp_get_thread_num();
sum += 2;
files[t] << B(x01) << B(x02) << B(x03) << B(x04) << B(x05)
<< B(x06) << B(x07) << B(x08) << B(x09) << B(x10)
<< B(x11) << B(x12) << B(x13) << B(x14) << B(x15)
<< B(x16) << B(x17) << B(x18) << B(x19) << B(x20)
<< B(x21) << B(x22) << B(x23) << B(x24) << B(x25);
files[t] << B(x07) << B(x06) << B(x08) << B(x10) << B(x09)
<< B(x02) << B(x01) << B(x03) << B(x05) << B(x04)
<< B(x12) << B(x11) << B(x13) << B(x15) << B(x14)
<< B(x22) << B(x21) << B(x23) << B(x25) << B(x24)
<< B(x17) << B(x16) << B(x18) << B(x20) << B(x19);
if (x13 != 13) {
sum += 2;
files[t] << C(x01) << C(x02) << C(x03) << C(x04) << C(x05)
<< C(x06) << C(x07) << C(x08) << C(x09) << C(x10)
<< C(x11) << C(x12) << C(x13) << C(x14) << C(x15)
<< C(x16) << C(x17) << C(x18) << C(x19) << C(x20)
<< C(x21) << C(x22) << C(x23) << C(x24) << C(x25);
files[t] << C(x07) << C(x06) << C(x08) << C(x10) << C(x09)
<< C(x02) << C(x01) << C(x03) << C(x05) << C(x04)
<< C(x12) << C(x11) << C(x13) << C(x15) << C(x14)
<< C(x22) << C(x21) << C(x23) << C(x25) << C(x24)
<< C(x17) << C(x16) << C(x18) << C(x20) << C(x19);
}}
FLIP(x16); }
FLIP(x15); }
FLIP(x11); }}
FLIP(x10); }
FLIP(x06); }}
FLIP(x18); }
FLIP(x08); }}
FLIP(x23); }
FLIP(x14); }
FLIP(x24); }}
FLIP(x22); }
FLIP(x12); }}
FLIP(x04); }
FLIP(x03); }}
FLIP(x02); }}
FLIP(x21); }
FLIP(x05); }}
FLIP(x25); }
FLIP(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