strassen算法c语言（strassen矩阵乘法算法C语言）

本文仅代码，无理论解释

实话实说，我觉得这个算法在C系列的语言下，简直垃圾到爆炸……毕竟是一群完全不懂程序数学家对着纸弄出来的，看起来好像非常的有用，实际上耗时是非常爆炸的。

但是《算法导论》里有啊……然后上课又要求手写一个

于是我就手写了一个……我尽可能的减少使用的空间同时加快速度了，当 n = 512 的时候，内存使用量峰值没有超过 10mb，而且是通过递归实现 Strassen 算法

其中，in.txt 已经预先准备了 3000000 个范围在 0-100 随机数，避免程序在运算过程中爆 int（虽然完全可以取1000）

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 /** * Created by Mauve on 3/29/2020. * Copyright © 2020 Mauve, All Rights Reserved */ #include <bits/stdc++.h> using namespace std; /** * 矩阵相乘 * 最终结果耗时结果保存至 * https://www.desmos.com/calculator/gl4tm5i1zu */ struct mat { unsigned row, col; mat(unsigned r, unsigned c) : row(r), col(c) {} virtual int &pos_ref(unsigned i, unsigned j) = 0; virtual int pos(unsigned i, unsigned j) const = 0; }; struct base_mat; struct sub_mat; stack<sub_mat *> sub_data; struct base_mat : mat { int *data; base_mat(unsigned r, unsigned c) : mat(r, c), data(new int[row * col]) {} ~base_mat() { delete[] data; } inline int &pos_ref(unsigned i, unsigned j) override { return *(data + i * col + j); } inline int pos(unsigned i, unsigned j) const override { return *(data + i * col + j); } }; unsigned min_mul; struct sub_mat : mat { mat *a, *b; bool is_add; unsigned offset_ai, offset_aj, offset_bi, offset_bj; explicit sub_mat(mat *data) : mat(data->row, data->col), a(data), b(nullptr), is_add(false), offset_ai(0), offset_aj(0), offset_bi(0), offset_bj(0) { sub_data.push(this); } sub_mat(mat *data, bool of_i, bool of_j) : mat(data->row >> 1u, data->col >> 1u), a(data), b(nullptr), is_add(false), offset_ai(of_i ? data->row >> 1u : 0), offset_aj(of_j ? data->col >> 1u : 0), offset_bi(0), offset_bj(0) { sub_data.push(this); } inline int &pos_ref(unsigned i, unsigned j) override { assert(b == nullptr); return a->pos_ref(i + offset_ai, j + offset_aj); } inline int pos(unsigned i, unsigned j) const override { if (b == nullptr) return a->pos(i + offset_ai, j + offset_aj); return a->pos(i + offset_ai, j + offset_aj) + (is_add ? 1 : -1) * b->pos(i + offset_bi, j + offset_bj); } inline sub_mat *operator+(sub_mat &other) { auto res = new sub_mat(this); res->b = &other; res->is_add = true; return res; } inline sub_mat *operator-(sub_mat &other) { auto res = new sub_mat(this); res->b = &other; res->is_add = false; return res; } mat *operator*(sub_mat &other) { assert(col == other.row); auto res = new base_mat(row, other.col); if (col & 1u || row & 1u || col <= min_mul || row <= min_mul || other.col <= min_mul) { memset(res->data, 0, sizeof(int) * res->row * res->col); for (int k = 0; k < col; k++) for (int i = 0; i < row; ++i) for (int j = 0; j < other.col; ++j) res->pos_ref(i, j) += pos(i, k) * other.pos(k, j); } else { size_t sub_data_size = sub_data.size(); #define a(i, j) (*new sub_mat(this, i == 2 , j == 2)) #define b(i, j) (*new sub_mat(&other, i == 2 , j == 2)) auto m1 = *(a(1, 1) + a(2, 2)) * *(b(1, 1) + b (2, 2)); auto m2 = *(a(2, 1) + a(2, 2)) * b(1, 1); auto m3 = a(1, 1) * *(b(1, 2) - b(2, 2)); auto m4 = a(2, 2) * *(b(2, 1) - b(1, 1)); auto m5 = *(a(1, 1) + a(1, 2)) * b(2, 2); auto m6 = *(a(2, 1) - a(1, 1)) * *(b(1, 1) + b(1, 2)); auto m7 = *(a(1, 2) - a(2, 2)) * *(b(2, 1) + b(2, 2)); #undef a #undef b unsigned half_row = row >> 1u, half_col = col >> 1u; #define m(t) (m##t->pos(i, j)) // C11 for (unsigned i = 0; i < half_row; ++i) for (unsigned j = 0; j < half_col; ++j) res->pos_ref(i, j) = m(1) + m(4) - m(5) + m(7); // C12 for (unsigned i = 0; i < half_row; ++i) for (unsigned j = 0; j < half_col; ++j) res->pos_ref(i, j + half_col) = m(3) + m(5); // C21 for (unsigned i = 0; i < half_row; ++i) for (unsigned j = 0; j < half_col; ++j) res->pos_ref(i + half_row, j) = m(2) + m(4); // C22 for (unsigned i = 0; i < half_row; ++i) for (unsigned j = 0; j < half_col; ++j) res->pos_ref(i + half_row, j + half_col) = m(1) - m(2) + m(3) + m(6); #undef m delete dynamic_cast<base_mat *>(m1); delete dynamic_cast<base_mat *>(m2); delete dynamic_cast<base_mat *>(m3); delete dynamic_cast<base_mat *>(m4); delete dynamic_cast<base_mat *>(m5); delete dynamic_cast<base_mat *>(m6); delete dynamic_cast<base_mat *>(m7); while (sub_data.size() > sub_data_size) { delete sub_data.top(); sub_data.pop(); } } return res; } }; unsigned N = 2; void solve() { cerr << "N = " << N << endl; base_mat a(N, N), b(N, N); for (int i = 0; i < N; ++i) for (int j = 0; j < N; ++j) cin >> a.pos_ref(i, j); for (int i = 0; i < N; ++i) for (int j = 0; j < N; ++j) cin >> b.pos_ref(i, j); for (int t = 1; t < min(10u, N); t += 3) { auto x = new sub_mat(&a), y = new sub_mat(&b); min_mul = t; auto time_1 = clock(); auto z = *x * *y; auto time_2 = clock(); cerr << "t = " << t << " time: " << double(time_2 - time_1) / CLOCKS_PER_SEC << endl; delete dynamic_cast<base_mat *>(z); while (!sub_data.empty()) { delete sub_data.top(); sub_data.pop(); } } auto x = new sub_mat(&a), y = new sub_mat(&b); min_mul = 10000; auto time_1 = clock(); auto z = *x * *y; auto time_2 = clock(); cerr << "tradition: " << double(time_2 - time_1) / CLOCKS_PER_SEC << endl; delete dynamic_cast<base_mat *>(z); while (!sub_data.empty()) { delete sub_data.top(); sub_data.pop(); } N *= 2; if (N >= 1000) exit(0); } signed main() { ios_base::sync_with_stdio(false); cin.tie(nullptr); cout.tie(nullptr); #ifdef ACM_LOCAL freopen("in.txt", "r", stdin); freopen("out.txt", "w", stdout); long long test_index_for_debug = 1; char acm_local_for_debug; while (cin >> acm_local_for_debug && acm_local_for_debug != '~') { cin.putback(acm_local_for_debug); if (test_index_for_debug > 20) { throw runtime_error("Check the stdin!!!"); } auto start_clock_for_debug = clock(); solve(); auto end_clock_for_debug = clock(); cout << "Test " << test_index_for_debug << " successful" << endl; cerr << "Test " << test_index_for_debug++ << " Run Time: " << double(end_clock_for_debug - start_clock_for_debug) / CLOCKS_PER_SEC << "s" << endl; cout << "--------------------------------------------------" << endl; } #else solve(); #endif return 0; }