Strassen算法代码

本文仅代码,无理论解释

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

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

于是我就手写了一个……我尽可能的减少使用的空间同时加快速度了,而且是通过递归实现 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 || other.col & 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;
}

Strassen算法代码
https://blog.mauve.icu/2021/10/13/cpp/Strassen/
作者
Shiroha
发布于
2021年10月13日
许可协议