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
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
| /*********************************************
* 多倍長乗算 ( by Toom-Cook 法 (3-way) )
* - 最下位の桁を配列の先頭とする考え方
*********************************************/
#include <cstdlib> // for rand()
#include <iostream> // for cout
#include <math.h> // for pow()
#include <stdio.h> // for printf()
//#define TEST // テスト ( 乗算回数・処理時間 ) するならコメント解除
#define D_MAX 729 // 計算可能な最大桁数 ( 3 のべき乗 )
#define D 729 // 実際に計算する桁数 ( D_MAX 以下 )
using namespace std;
/*
* 計算クラス
*/
class Calc
{
int A[D]; // 被乗数配列
int B[D]; // 乗数配列
#ifdef TEST
int cnt_mul; // 乗算回数
clock_t t1, t2; // 計算開始CPU時刻、計算終了CPU時刻
double tt; // 計算時間
#endif
public:
Calc(); // コンストラクタ
void calcToomCook(); // 計算 ( Toom-Cook 法 )
private:
void multiplyNormal(int *, int *, int, int *); // 乗算 ( 標準(筆算)法 )
void multiplyToomCook3(int *, int *, int , int *); // 乗算 ( Toom-Cook 法 (3-way) )
void doCarry(int *, int); // 繰り上がり・借り処理
void display(int *, int *, int *); // 結果出力
};
/*
* コンストラクタ
*/
Calc::Calc()
{
/* ====================================== *
* テストなので、被乗数・乗数は乱数を使用 *
* ====================================== */
int i; // LOOP インデックス
// 被乗数・乗数桁数設定
for (i = 0; i < D; i++) {
A[i] = rand() % 10;
B[i] = rand() % 10;
}
}
/*
* 計算 ( Toom-Cook 法 )
*/
void Calc::calcToomCook()
{
int a[D_MAX]; // 被乗数配列
int b[D_MAX]; // 乗数配列
int z[D_MAX * 2]; // 計算結果用配列
int i; // LOOPインデックス
#ifdef TEST
t1 = clock(); // 計算開始時刻
for (int l = 0; l < 1000; l++) {
cnt_mul = 0; // 乗算回数リセット
#endif
// 配列初期設定 ( コンストラクタで設定した配列を設定 )
for (i = 0; i < D; i++) {
a[i] = A[i];
b[i] = B[i];
}
// 最大桁に満たない部分は 0 を設定
for (i = D; i < D_MAX; i++) {
a[i] = 0;
b[i] = 0;
}
// 乗算 ( Toom-Cook 法 (3-way) )
multiplyToomCook3(a, b, D_MAX, z);
// 繰り上がり・借り処理
doCarry(z, D_MAX * 2);
#ifdef TEST
}
t2 = clock(); // 計算終了時刻
tt = (double)(t2 - t1) / CLOCKS_PER_SEC; // ==== 計算時間
#endif
// 結果出力
display(a, b, z);
}
/*
* 乗算 ( 標準(筆算)法 )
*/
void Calc::multiplyNormal(int *a, int *b, int tLen, int *z)
{
int i, j; // ループインデックス
// 計算結果初期化
for(i = 0; i < tLen * 2; i++) z[i] = 0;
// 各配列を各桁とみなして乗算
for (j = 0; j < tLen; j++) {
for (i = 0; i < tLen; i++) {
z[j + i] += a[i] * b[j];
#ifdef TEST
cnt_mul++; // 乗算カウント
#endif
}
}
}
/*
* 乗算 ( Toom-Cook 法 (3-way) )
* 結果用配列は以下のように配置し、
* +----+----+----+----+----+----+----+----+----+----+
* | c0 | c2 | c4 | c1 | c3 |
* +----+----+----+----+----+----+----+----+----+----+
* 最後に、c1, c3 を所定の位置に加算する。
* +----+----+----+----+----+----+
* | c0 | c2 | c4 |
* +----+----+----+----+----+----+
* +----+----+----+----+
* | c1 | c3 |
* +----+----+----+----+
* その他、計算に必要な配列ポインタを詳細に設定。
*/
void Calc::multiplyToomCook3(int *a, int *b, int tLen, int *z)
{
// ==== 変数宣言
int *a0 = &a[0]; // 被乗数/右側配列ポインタ
int *a1 = &a[tLen / 3]; // 被乗数/中央配列ポインタ
int *a2 = &a[tLen * 2/ 3]; // 被乗数/左側配列ポインタ
int *b0 = &b[0]; // 乗数 /右側配列ポインタ
int *b1 = &b[tLen / 3]; // 乗数 /中央配列ポインタ
int *b2 = &b[tLen * 2/ 3]; // 乗数 /左側配列ポインタ
int *c0 = &z[(tLen / 3) * 0]; // c0 用配列ポインタ
int *c2 = &z[(tLen / 3) * 2]; // c2 用配列ポインタ
int *c4 = &z[(tLen / 3) * 4]; // c4 用配列ポインタ
int c1 [(tLen / 3) * 2]; // c1 用配列
int c3 [(tLen / 3) * 2]; // c3 用配列
int a_m2 [tLen / 3]; // a( -2) 用配列
int a_m1 [tLen / 3]; // a( -1) 用配列
int a_0 [tLen / 3]; // a( 0) 用配列
int a_1 [tLen / 3]; // a( 1) 用配列
int a_inf [tLen / 3]; // a(inf) 用配列
int b_m2 [tLen / 3]; // b( -2) 用配列
int b_m1 [tLen / 3]; // b( -1) 用配列
int b_0 [tLen / 3]; // b( 0) 用配列
int b_1 [tLen / 3]; // b( 1) 用配列
int b_inf [tLen / 3]; // b(inf) 用配列
int c_m2 [(tLen / 3) * 2]; // c( -2) 用配列
int c_m1 [(tLen / 3) * 2]; // c( -1) 用配列
int c_0 [(tLen / 3) * 2]; // c( 0) 用配列
int c_1 [(tLen / 3) * 2]; // c( 1) 用配列
int c_inf [(tLen / 3) * 2]; // c(inf) 用配列
int i; // LOOPインデックス
// ==== 9 桁（配列 9 個）になった場合は標準乗算
if (tLen <= 9) {
multiplyNormal(a, b, tLen, z);
return;
}
// ==== a(-2) = 4 * a2 - 2 * a1 + a0, b(1) = 4 * b2 - 2 * b1 + b0 (by シフト演算)
for(i = 0; i < tLen / 3; i++) {
a_m2[i] = (a2[i] << 2) - (a1[i] << 1) + a0[i];
b_m2[i] = (b2[i] << 2) - (b1[i] << 1) + b0[i];
}
// ==== a(-1) = a2 - a1 + a0, b(1) = b2 - b1 + b0
for(i = 0; i < tLen / 3; i++) {
a_m1[i] = a2[i] - a1[i] + a0[i];
b_m1[i] = b2[i] - b1[i] + b0[i];
}
// ==== a(0) = a0, b(0) = b0
for(i = 0; i < tLen / 3; i++) {
a_0[i] = a0[i];
b_0[i] = b0[i];
}
// ==== a(1) = a2 + a1 + a0, b(1) = b2 + b1 + b0
for(i = 0; i < tLen / 3; i++) {
a_1[i] = a2[i] + a1[i] + a0[i];
b_1[i] = b2[i] + b1[i] + b0[i];
}
// ==== a(inf) = a2, b(inf) = b2
for(i = 0; i < tLen / 3; i++) {
a_inf[i] = a2[i];
b_inf[i] = b2[i];
}
// ==== c(-2) = a(-2) * b(-2)
multiplyToomCook3(a_m2, b_m2, tLen / 3, c_m2 );
// ==== c(-1) = a(-1) * b(-1)
multiplyToomCook3(a_m1, b_m1, tLen / 3, c_m1 );
// ==== c(0) = a(0) * b(0)
multiplyToomCook3(a_0, b_0, tLen / 3, c_0 );
// ==== c(1) = a(1) * b(1)
multiplyToomCook3(a_1, b_1, tLen / 3, c_1 );
// ==== c(inf) = a(inf) * b(inf)
multiplyToomCook3(a_inf, b_inf, tLen / 3, c_inf);
// ==== c4 = 6 * c(inf) / 6
for(i = 0; i < (tLen / 3) * 2; i++)
c4[i] = c_inf[i];
// ==== c3 = -c(-2) + 3 * c(-1) - 3 * c(0) + c(1) + 12 * c(inf) / 6
for(i = 0; i < (tLen / 3) * 2; i++) {
c3[i] = -c_m2[i];
c3[i] += (c_m1[i] << 1) + c_m1[i];
c3[i] -= (c_0[i] << 1) + c_0[i];
c3[i] += c_1[i];
c3[i] += (c_inf[i] << 3) + (c_inf[i] << 2);
c3[i] /= 6;
}
// ==== c2 = 3 * c(-1) - 6 * c(0) + 3 * c(1) - 6 * c(inf) / 6
for(i = 0; i < (tLen / 3) * 2; i++) {
c2[i] = (c_m1[i] << 1) + c_m1[i];
c2[i] -= (c_0[i] << 2) + (c_0[i] << 1);
c2[i] += (c_1[i] << 1) + c_1[i];
c2[i] -= (c_inf[i] << 2) + (c_inf[i] << 1);
c2[i] /= 6;
}
// ==== c1 = c(-2) - 6 * c(-1) + 3 * c(0) + 2 * c(1) - 12 * c(inf) / 6
for(i = 0; i < (tLen / 3) * 2; i++) {
c1[i] = c_m2[i];
c1[i] -= (c_m1[i] << 2) + (c_m1[i] << 1);
c1[i] += (c_0[i] << 1) + c_0[i];
c1[i] += (c_1[i] << 1);
c1[i] -= (c_inf[i] << 3) + (c_inf[i] << 2);
c1[i] /= 6;
}
// ==== c0 = 6 * c(0) / 6
for(i = 0; i < (tLen / 3) * 2; i++)
c0[i] = c_0[i];
// ==== z = c4 * x^4 + c3 * x^3 + c2 * x^2 + c1 * x + c0
// (c0, c2, c4 は最初から所定の位置に格納されているので、
// c1, c3 のみ加算 )
for(i = 0; i < (tLen / 3) * 2; i++) z[i + tLen / 3] += c1[i];
for(i = 0; i < (tLen / 3) * 2; i++) z[i + (tLen / 3) * 3] += c3[i];
}
/*
* 繰り上がり・借り処理
*/
void Calc::doCarry(int *a, int tLen) {
int cr; // 繰り上がり
int i; // ループインデックス
cr = 0;
for(i = 0; i < tLen; i++) {
a[i] += cr;
if(a[i] < 0) {
cr = -(-(a[i] + 1) / 10 + 1);
} else {
cr = a[i] / 10;
}
a[i] -= cr * 10;
}
// オーバーフロー時
if (cr != 0) printf("[ OVERFLOW!! ] %d\n", cr);
}
/*
* 結果出力
*/
void Calc::display(int *a, int *b, int *z)
{
int i; // LOOPインデックス
// 上位桁の不要な 0 を削除するために、配列サイズを取得
int aLen = D_MAX, bLen = D_MAX, zLen = D_MAX * 2;
while (a[aLen - 1] == 0) if (a[aLen - 1] == 0) aLen--;
while (b[bLen - 1] == 0) if (b[bLen - 1] == 0) bLen--;
while (z[zLen - 1] == 0) if (z[zLen - 1] == 0) zLen--;
// a 値
printf("a =\n");
for (i = aLen - 1; i >= 0; i--) {
printf("%d", a[i]);
if ((aLen - i) % 10 == 0) printf(" ");
if ((aLen - i) % 50 == 0) printf("\n");
}
printf("\n");
// b 値
printf("b =\n");
for (i = bLen - 1; i >= 0; i--) {
printf("%d", b[i]);
if ((bLen - i) % 10 == 0) printf(" ");
if ((bLen - i) % 50 == 0) printf("\n");
}
printf("\n");
// z 値
printf("z =\n");
for (i = zLen - 1; i >= 0; i--) {
printf("%d", z[i]);
if ((zLen - i) % 10 == 0) printf(" ");
if ((zLen - i) % 50 == 0) printf("\n");
}
printf("\n\n");
#ifdef TEST
printf("Counts of multiply / 1 loop = %d\n", cnt_mul); // 乗算回数
printf("Total time of all loops = %f seconds\n", tt); // 処理時間
#endif
}
/*
* メイン処理
*/
int main()
{
try
{
// 計算クラスインスタンス化
Calc objCalc;
// 乗算 ( Toom-Cook 法 )
objCalc.calcToomCook();
}
catch (...) {
cout << "例外発生！" << endl;
return -1;
}
// 正常終了
return 0;
}
|
