-
-
Save taroyabuki/b94ed5b24a385e460af7 to your computer and use it in GitHub Desktop.
スパコンで約2時間36分かかったという、5×5の魔方陣の全解列挙を、パソコンで試す（C++）http://blog.unfindable.net/archives/7179
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| #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