-
Notifications
You must be signed in to change notification settings - Fork 2
/
karatsuba.cpp
43 lines (35 loc) · 1.62 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
#include <string.h>
int aux[16][3][100010]; // First dimension should be log_2(MAX_DEGREE)
// third should be (2 * MAX_DEGREE)
int MOD; // All operations are done % MOD. Define the MOD
const int LIM = 35; // Magical constant, in case of TLE try playing with it
// It assumes that r[] is initialized with 0, p1[] and p2[] are normalized
// (i.e.: %MOD), g1 = degree(p1), g2 = degree(p2). Don't call with _lvl
// coefficient of x^i is in index i of the array
// Complexity: O(max(g1, g2) ^ lg_2 3) = O(max(g1,g2) ^ 1.58)
void karatsuba(int r[], int g1, int p1[], int g2, int p2[], int _lvl = 0) {
if (g1 <= LIM || g2 <= LIM) { // if one is small, do O(g1g2)
for (int i = 0; i <= g1; i++)
for (int j = 0; j <= g2; j++)
r[i+j] = (r[i+j] + p1[i]*p2[j])%MOD;
return;
}
int degree = max(g1, g2);
int degree2 = degree/2;
// A0*B0
karatsuba(r, degree2, p1, degree2, p2, _lvl+1);
// A1*B1
karatsuba(r+2*degree2+2, g1-degree2-1, p1+degree2+1, g2-degree2-1, p2+degree2+1, _lvl+1);
// (A0+A1)*(B0+B1)
for (int i = 0; i <= degree2; i++)
aux[_lvl][0][i] = (p1[i] + ((i + degree2 + 1 <= g1)?p1[i+degree2+1]:0)) % MOD;
for (int i = 0; i <= degree2; i++)
aux[_lvl][1][i] = (p2[i] + ((i + degree2 + 1 <= g2)?p2[i+degree2+1]:0)) % MOD;
memset(aux[_lvl][2], 0, (degree + 2)*sizeof(int));
karatsuba(aux[_lvl][2], degree2, aux[_lvl][0], degree2, aux[_lvl][1], _lvl+1);
// puts result in r[]
for (int i = 0; i <= 2*degree2; i++)
aux[_lvl][2][i] = (aux[_lvl][2][i] + 2 * MOD - r[i] - r[i+2*degree2+2]) % MOD;
for (int i = 0; i <= 2*degree2; i++)
r[degree2 + i + 1] = (r[degree2 + i + 1] + aux[_lvl][2][i]) % MOD;
}