数论-FTT 和 NTT

FTT

 

NKOJ3071 高精度乘（输入 a, b ，输出两个数的积）

 

 1 #include <stdio.h>
 2 #include <complex>
 3 #include <math.h>
 4 
 5 using namespace std;
 6 
 7 typedef complex<double> CP;
 8 typedef long long LL;
 9 
10 const int _N = 300005;
11 const double PI = asin(1)*2;
12 
13 LL rev[_N], Ans[_N];
14 char str1[_N], str2[_N];
15 CP X[_N], Y[_N];
16 
17 void GetRev(LL bit)
18 {
19     for (LL i = 0; i < (1<<bit); ++i)
20         rev[i] = (rev[i>>1]>>1) | ((i&1)<<(bit-1));
21     return;
22 }
23 
24 void FFT(CP *A, LL n, LL ty)
25 {
26     LL i, k, len;
27     for (i = 0; i < n; ++i)
28         if (i < rev[i]) swap(A[i], A[rev[i]]);
29         
30     for (len = 1; len < n; len <<= 1) {
31         CP wn = exp(CP(0, ty*PI/len));
32         for (i = 0; i < n; i += len<<1) {
33             CP wi(1, 0);
34             for (k = i; k < i+len; ++k) {
35                 CP t0 = A[k], t1 = wi*A[k+len];
36                 A[k] = t0+t1, A[k+len] = t0-t1;
37                 wi *= wn;
38             }
39         }
40     }
41     if (ty == -1)
42         for (i = 0; i < n; ++i) A[i] /= n;
43     return;
44 }
45 
46 bool Input(LL &len, char *str)
47 {
48     char tt; bool flag = false; int i;
49     while (((tt = getchar()) < '0' || tt > '9') && tt != '-');
50     if (tt == '-') flag = true, i = -1;
51     else str[0] = tt, i = 0;
52     while ((tt = getchar()) >= '0' && tt <= '9') str[++i] = tt;
53     len = i+1;
54     return flag;
55 }
56 
57 int main()
58 {
59     LL flag1, flag2, l1, l2, i;
60     flag1 = Input(l1, str1), flag2 = Input(l2, str2);
61 //    printf("\n\n%s %s\n\n %lld %lld\n", str1, str2, l1, l2);
62     LL a = 1, x = 0;
63     while (a < l1+l2-1) a <<= 1, ++x;
64 //    printf("\n\na = %lld, x = %lld\n\n", a, x);
65     for (i = 0; i < l1; ++i) X[i] = (double)(str1[l1-1-i]-'0');
66     for (i = 0; i < l2; ++i) Y[i] = (double)(str2[l2-1-i]-'0');
67     GetRev(x), FFT(X, a, 1), FFT(Y, a, 1);
68     for (i = 0; i < a; ++i) X[i] *= Y[i];
69     FFT(X, a, -1);
70     for (i = 0; i < l1+l2; ++i)
71         Ans[i] += (long long)(X[i].real() + 0.5), Ans[i+1] += Ans[i]/10, Ans[i] %= 10;
72     for (i = l1+l2; i >= 0 && !Ans[i]; --i);
73     if (i == -1) { printf("0\n"); return 0; }
74     if (flag1 ^ flag2) putchar('-');
75     while (i >= 0) putchar(Ans[i] + '0'), --i;
76     putchar('\n');
77     return 0;
78 }

 

 

NTT

1.数论阶和原根

 

相关文章：数论之原根 https://blog.csdn.net/fuyukai/article/details/50894609

 

ordn a 表示当 a, n 互素时，在模 n 意义下，a 的数论阶，即满足 ax Ξ 1 (mod n) 的最小正整数 x .

结论：(ordn a) | x

推论：(ordn a) | phi(x)

 

当 n > 0 且 ordn a = phi(n) 时，a 为 n 的一个原根.

素数一定有原根，整数不一定. 当一个数 n 有原根时，它的原根数量是 phi(phi(n)) .

在模 n 意义下，n 的原根为 g，gi 互不相同 , 0 <= i <= phi(n) - 1

这一点性质类似 FFT 中的单位根 w ，也就是 NTT 的基础。

 

求素数 n 的原根

从小到大枚举每一个 x | phi(n) , 判断是否满足 ax Ξ 1 (mod n)，找到一个合理的 x 即为 n 的原根.

 

2.NTT操作

 

 1 #include <stdio.h>
 2 #include <algorithm>
 3 
 4 using namespace std;
 5 
 6 typedef long long LL;
 7 
 8 const int _N = 300005;
 9 const long long P = 998244353LL;
10 const long long G = 3LL;
11 
12 LL rev[_N], Ans[_N];
13 char str1[_N], str2[_N];
14 LL X[_N], Y[_N];
15 
16 void GetRev(LL bit)
17 {
18     for (LL i = 0; i < (1<<bit); ++i)
19         rev[i] = (rev[i>>1]>>1) | ((i&1)<<(bit-1));
20     return;
21 }
22 
23 LL Mont(LL t1, LL t2)
24 {
25     t1 %= P;
26     LL ans = 1;
27     while (t2) {
28         if (t2 & 1) ans = ans*t1%P;
29         t2 >>= 1, t1 = t1*t1%P;
30     }
31     return ans;
32 }
33 
34 void NTT(LL *A, LL n, LL ty)
35 {
36     LL i, k, len;
37     for (i = 0; i < n; ++i)
38         if (i < rev[i]) swap(A[i], A[rev[i]]);
39         
40     for (len = 1; len < n; len <<= 1) {
41         LL wn = Mont(G, (P-1)+ty*(P-1)/(len<<1));
42         for (i = 0; i < n; i += len<<1) {
43             LL wi = 1;
44             for (k = i; k < i+len; ++k) {
45                 LL t0 = A[k], t1 = wi*A[k+len]%P;
46                 A[k] = t0+t1, A[k+len] = t0-t1;
47                 if ((A[k] = t0+t1) >= P) A[k] -= P;
48                 if ((A[k+len] = t0-t1) < 0) A[k+len] += P;
49                 wi = wi*wn%P;
50             }
51         }
52     }
53     if (ty == -1) {
54         LL tmp = Mont(n, P-2);
55         for (i = 0; i < n; ++i) A[i] = A[i]*tmp%P;
56     }
57     return;
58 }
59 
60 bool Input(LL &len, char *str)
61 {
62     char tt; bool flag = false; int i;
63     while (((tt = getchar()) < '0' || tt > '9') && tt != '-');
64     if (tt == '-') flag = true, i = -1;
65     else str[0] = tt, i = 0;
66     while ((tt = getchar()) >= '0' && tt <= '9') str[++i] = tt;
67     len = i+1;
68     return flag;
69 }
70 
71 int main()
72 {
73     LL flag1, flag2, l1, l2, i;
74     flag1 = Input(l1, str1), flag2 = Input(l2, str2);
75 //    printf("\n\n%s %s\n\n %lld %lld\n", str1, str2, l1, l2);
76     LL a = 1, x = 0;
77     while (a < l1+l2-1) a <<= 1, ++x;
78 //    printf("\n\na = %lld, x = %lld\n\n", a, x);
79     for (i = 0; i < l1; ++i) X[i] = str1[l1-1-i]-'0';
80     for (i = 0; i < l2; ++i) Y[i] = str2[l2-1-i]-'0';
81     GetRev(x), NTT(X, a, 1), NTT(Y, a, 1);
82     for (i = 0; i < a; ++i) X[i] = X[i]*Y[i]%P;
83     NTT(X, a, -1);
84     for (i = 0; i < l1+l2; ++i)
85         Ans[i] += X[i], Ans[i+1] += Ans[i]/10, Ans[i] %= 10;
86     for (i = l1+l2; i >= 0 && !Ans[i]; --i);
87     if (i == -1) { printf("0\n"); return 0; }
88     if (flag1 ^ flag2) putchar('-');
89     while (i >= 0) putchar(Ans[i] + '0'), --i;
90     putchar('\n');
91     return 0;
92 }

 

posted @ 2018-05-02 17:34  derchg  阅读(1034)  评论(0)  编辑  收藏  举报