/magicsquare.cpp
Secret
Last active Mar 6, 2018
スパコンで約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