-
Notifications
You must be signed in to change notification settings - Fork 2
/
Karatsuba.cpp
76 lines (64 loc) · 1.82 KB
/
Karatsuba.cpp
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
/**
* Description: Karatsuba Algorithm (Fast polynomial multiplication)
* Usage: multiply O(N^1.583)
* Source: https://github.com/dragonslayerx
*/
#include <iostream>
#include <cstdio>
using namespace std;
#define MOD 99991
//m should be a power of 2
class Karatsuba {
public:
static void multiply(int *A, int *B, int *C, int lA, int rA, int lB, int rB){
int m = rA-lA+1;
if (m == 1) {
C[0] = ((long long)A[lA]*B[lB]) % MOD;
return;
}
int z0[m], z1[m], z2[m];
int midA = (lA + rA) >> 1;
int midB = (lB + rB) >> 1;
multiply(A, B, z0, lA, midA, lB, midB);
multiply(A, B, z1, midA+1, rA, midB+1, rB);
int a[m], b[m];
int shift = m>>1;
int mid = m>>1;
for (int i = lA, j = 0; i <= midA; i++, j++) {
a[j] = A[i] + A[i+shift];
if (a[j] >= MOD) a[j] -= MOD;
}
for (int i = lB, j = 0; i <= midB; i++, j++) {
b[j] = B[i] + B[i+shift];
if (b[j] >= MOD) b[j] -= MOD;
}
multiply(a, b, z2, 0, mid-1, 0, mid-1);
for (int i = 0; i <= m-2; i++) {
C[i] = z0[i];
if (C[i] >= MOD) C[i] -= MOD;
}
C[m-1] = 0;
shift = m;
for (int i = 0; i <= m-2; i++) {
C[i+shift] = z1[i];
if (C[i+shift] >= MOD) C[i+shift] -= MOD;
}
shift = m>>1;
for (int i = 0; i <= m-2; i++) {
C[i+shift] += (z2[i] + (MOD-z1[i]) + (MOD-z0[i])) % MOD;
if (C[i+shift] >= MOD) {
C[i+shift] -= MOD;
}
}
}
};
int main()
{
int A[] = {1, 1, 1, 1};
int B[] = {1, 1, 0, 0};
int C[7] = {};
Karatsuba::multiply(A, B, C, 0, 3, 0, 3);
for (int i = 0; i < 7; i++) {
cerr << C[i] << endl;;
}
}