ACM►学习笔记►DP

树形dp

  1. 给你一棵n个节点的树，要你求出包含每个点的连通点集的数量。答案对1e9+7取模。

两次dfs，第一次处理每个节点的子孙方向的联通点集数量，
$$
dp[i]=\prod _{all..son}{(dp[son]+1)}
$$

第二次处理父亲方向的点集数量
$$
\begin{aligned}
&dp2[i] \
=&(\prod _{all..brother}{(dp1[brother]+1)})*dp2[father]+1 \
=&\frac{dp1[father]}{dp1[i]+1}*dp2[father]+1 \
\end{aligned}
$$

特别注意这里小心没有逆元

Others

排序算法比较 

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
#include<iostream>
#include<iomanip>
#include<vector>
#include<cstring>
using namespace std;

//快速排序
//unstable
//time:   sita(n*lgn)   omiga(n*lgn)   O(n^2)
//memery  O(lgn)
void quicksort(int *a, int *b){  //[a,b)
    int *l = a , *r = b-1;
    if(l>=r)return;
    int *mid = l; //suppose that the little some is in [l,mid-1] and bound is in mid
    while(r>mid){
        if(*r<*l)swap(*(++mid),*r);
        else r--;
    }
    swap(*l,*mid);
    quicksort(a,mid);
    quicksort(mid+1,b);
}


//插入排序
//stable
//time:   sita(n^2) omiga(n)  O(n^2)
//memery: O(1)
void insertionsort(int *a,int *b){
    for(int *i = a; i<b; i++){
        int *j = i;
        while(j>a&&*j<*(j-1)){//可能每次都直接跳出去，结果成了omiga(n)
            swap(*j,*(j-1));
            j--;
        }
    }
}


//选择排序
//unstable
//time:   sita(n^2) omiga(n^2)  O(n^2)
//memery: O(1)
void selectsort(int *a,int *b){
    for(int *i=a;i<b;i++){
        int*min = i;
        for(int *j=i+1;j<b;j++){
            if(*min>*j)min=j;
        }
        swap(*min,*i);
    }
}


//基数排序??
//unstable
//time:   omige(log(k,maxvalue)*n)  O(log(k,maxvalue)*n)
//memery:
void radixsort(int*a,int*b){
    const int k =15;
    vector<int>bin[2][k+1];// 滚动数组
    int t=0;
    
    for(int*i=a;i<b;i++)bin[t][(*i)&k].push_back(*i);
    t^=1;
    
    for(int times=1;times<8;times++){
        for(int i=0;i<=k;i++)bin[t][i].clear();
        for(int i=0;i<=k;i++){
            for(int x:bin[t^1][i]){
                bin[t][(x>>(4*times))&k].push_back(x);
            }
        }
        t^=1;
    }
    for(int i=0;i<k;i++){
        for(int x:bin[t^1][i]){
            *(a++)=x;
        }
    }
}


//计数排序??
//unstable
//time:   sita(maxvalue-minvalue) omiga(maxvalue-minvalue) O(maxvalue-minvalue)
//memery: O(maxvalue-minvalue)
void countingsort(int*a,int*b){
    int maxvalue = (1<<31); // 最小的int
    int minvalue = -1^(1<<31);// 最大的int
    for(int*i=a;i<b;i++){
        maxvalue=max(maxvalue,*i);
        minvalue=min(minvalue,*i);
    }
    int*data = new int[maxvalue-minvalue+1];
    memset(data,0,sizeof(int)*(maxvalue-minvalue+1));
    for(int*i=a;i<b;i++)data[(*i)-minvalue]++;
    for(int i=minvalue;i<=maxvalue;i++){
        while(data[i-minvalue]--)*(a++)=i;
    }
    delete[] data;
}


//桶排序??
//time:   sita(k*O(sort(n/k))  omiga(k*O(sort(n/k))  O(k*O(sort(n/k))
//memery:
void bucketsort(int*a,int*b){
    int maxvalue = (1<<31); // 最小的int
    int minvalue = -1^(1<<31);// 最大的int
    for(int*i=a;i<b;i++){
        maxvalue=max(maxvalue,*i);
        minvalue=min(minvalue,*i);
    }
    if(maxvalue==minvalue)return;
    const int k = 10;
    vector<int>data[k];
    for(int*i=a;i<b;i++){
        int belong = ((*i)-minvalue)/((maxvalue-minvalue)/10+1);
        data[belong].push_back(*i);
    }
    for(int i=0;i<k;i++){
        sort(data[i].begin(),data[i].end());
        for(int x:data[i])*(a++)=x;
    }
}


//冒泡排序
//stable
//time： sita(n^2) omiga（n） O(n^2)
//memery: O(1)
void bubblesort(int*a,int *b){
    for(int i=0;i<b-a;i++){
        bool noswap = true;
        for(int*j=a+1;j<b-i;j++){
            if(*(j-1)>*j){
                swap(*(j-1),*j);
                noswap = false;
            }
        }
        if(noswap)break;
    }
}


//归并排序
//stable
//time:   sita(n*lgn) omiga(n*lgn) O(n*lgn)
//memery: O(n)
void mergesort(int*a,int*b){
    int*mid = a+ (b-a)/2; // =(a+b)/2
    if(a>=b-1)return;
    
    mergesort(a,mid);
    mergesort(mid,b);
    
    int*c=new int[b-a];
    int*pa=a,*pm=mid,*pc=c;
    while(pa<mid&&pm<b){
        *(pc++)=*pa<=*pm?*(pa++):*(pm++);
    }
    while(pa<mid)*(pc++)=*(pa++);
    while( pm<b )*(pc++)=*(pm++);
    memcpy(a,c,sizeof(int)*(b-a));
}


//希尔排序
//unstable
//time:   sita(n^1.3) omiga(??) O(n^2)
//memery: O(1)
void shellsort(int*a,int*b){
    for(int gap=(b-a)>>1;gap;gap>>=1){//枚举gap
        for(int*i=a+gap;i<b;i++){
            int x=*i,*j=i;
            while(j-gap>=a&&x<*(j-gap)){
                *j=*(j-gap);
                j-=gap;
            }
            *j=x;
        }
    }
}


//堆排序
//unstable
//time:   sita(n*lgn)   omiga(n*lgn)   O(n*lgn)
//memery: O(1)
void heapsort(int*a,int*b){
    const int n=b-a;
    a--;
    
    //up
    for(int i=1;i<=n;i++){//建立大根堆
        int j=i;
        while(j!=1&&a[j]>a[j>>1]){//不是根且比父亲大
            swap(a[j],a[j>>1]);
            j>>=1;
        }
    }
   
    //down
    for(int i=n;i>=1;i--){
        swap(a[1],a[i]);
        int cur = 1;
        while((cur<<1)<i){//左儿子存在
            int son = (cur<<1|1)<i && a[cur<<1|1]>a[cur<<1];//若满足此条件则son为1，代表着右儿子优先
            if(a[cur]<a[cur<<1|son])swap(a[cur],a[cur<<1|son]);
            else break;//若优先的儿子都不能换
            cur = cur<<1|son;
        }
    }
}


const int n = 1e4;
int ans[n];

inline void compare(int *a,int *b,void sort(int*,int*),string name){
    int data[n];
    memcpy(data,a,sizeof(int)*(b-a));
    time_t begin = clock();
    sort(data, data+n);
    time_t spend = clock()-begin;
    cout<<setw(16)<<std::left<<name<<" time:"<<setw(13)<<spend;
    //check it
    for(int i=0;i<n;i++){
        if(ans[i]==data[i]&&i==n-1)cout<<"everything is right"<<endl;
        else if(ans[i]!=data[i]){
            cout<<"error"<<endl;
            break;
        }
    }
}

int main(){
    srand(int(time(NULL)));

    int a[n];
    for(int i=0;i<n;i++)a[i]=rand()%9000000+1000000;
    memcpy(ans,a,sizeof(a));
    sort(ans,ans+n);
    
    compare(a,a+n,sort,"* STL(fastest)");
    compare(a,a+n,quicksort,"quicksort");
    compare(a,a+n,insertionsort,"insertionsort");
    compare(a,a+n,selectsort,"selectsort");
    compare(a,a+n,radixsort,"radixsort");
    compare(a,a+n,countingsort,"countingsort");
    compare(a,a+n,bucketsort,"bucketsort");
    compare(a,a+n,bubblesort,"bubblesort");
    compare(a,a+n,mergesort,"mergesort");
    compare(a,a+n,shellsort,"shellsort");
    compare(a,a+n,heapsort,"heapsort");
    
    puts("ok");
}



/*
 
 输出:
 
 * STL(fastest)   time:468          everything is right
 quicksort        time:1265         everything is right
 insertionsort    time:178130       everything is right
 selectsort       time:126705       everything is right
 radixsort        time:3128         everything is right
 countingsort     time:50445        everything is right
 bucketsort       time:1002         everything is right
 bubblesort       time:355497       everything is right
 mergesort        time:3295         everything is right
 shellsort        time:3960         everything is right
 heapsort         time:2347         everything is right
 ok
 
 
 */
ACM►学习笔记►字符串

Manacher算法 

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
struct Manacher {//鉴于马拉车算法较复杂，此处有少量修改，
    //s[i]=ma[i<<1]
    //mp[i]表示以i为中心的最长回文串的半径，且mp[i]-1恰好为此回文串包含原字符串的字符的数量
    //可以证明ma字符串所包含的回文串总数=原字符串b所包含的回文串总数+2n+2
    static const int maxn = 1e6 + 666;
    char ma[maxn << 1];
    int mp[maxn << 1], begat[maxn];//begta[i]-> 以i开头的回文串的数量  begin at

    void build(char *str) {
        int len = strlen(str + 1), l = 0;
        ma[l++] = '$';//$#.#.#.#.#.#.#.#
        ma[l++] = '#';
        for (int i = 1; i <= len; i++) {
            ma[l++] = str[i];
            ma[l++] = '#';
        }
        ma[l] = mp[l] = 0;
        int mx = 0, id = 0;
        for (int i = 0; i < l; i++) {
            mp[i] = mx > i ? min(mp[2 * id - i], mx - i) : 1;
            while (ma[i + mp[i]] == ma[i - mp[i]])mp[i]++;
            if (i + mp[i] > mx) {
                mx = i + mp[i];
                id = i;
            }
        }
        //for(int i=2;i<=l;i++)palindrome+=mp[i]>>1;//回文串个数

        //若不用dalt数组，此后可删掉
        for (int i = 1; i <= len; i++)begat[i] = 0;
        for (int i = 2; i < l; i++) {
            int s = i - mp[i] + 1;//ma串最长回文左端点s
            s = (s + 1) / 2;//变为str串最长回文左端点，向上取整,因为str[i]对应smp[i<<1]
            int t = s + mp[i] / 2 - 1;//右端点
            begat[s]++;
            begat[t + 1]--;
        }
        for (int i = 1; i <= len + 1; i++)begat[i] += begat[i - 1];//+1是为了还原
    }
};
ACM►学习笔记►字符串

AC自动机

AC自动机

所谓AC自动机，其实是kmp算法与字典树的结合,不懂这两个，是无法学会的。

自动机

自动机，是一个五元组，包括了状态的非空有穷集合，符号的有限集合，状态转移函 数， 开始状态，终止状态集合，而在字典树上，增加了两个新的东西，一个标记了终止状态集合，另一个辅助了状态转移函数。 我们利用字典树上 的节点来表示状态，而边则用来表示状态转移函数的一部分。

Read More

ACM►学习笔记►字符串

回文自动机

回文自动机，就像AC自动机一样，他采取字典树来储存状态集合，也像AC自动机是KMP 算法与字典树结合一样，回文自动机就是manacher算法和字典树结合的新算法。

在回文自动机里面，状态指的是，以字典树上所表示的字符串的逆序串的以末端字符镜像对称后得到的新串,简单说，baab在字典树上为root->a>b,cacaacac在字典树上为root->a->c->a->c,想像根就是对称中心，往儿子走就意味着，在对称中心两端加上同样的数字。当然这样子是无法表示奇数回文的，所以我们设立两个根就可以了。一个串的回文自动机，储存了这个串的所有回文子串。

回文自动机的状态转移的依据是回文，举个例子，如下图，如果黑色串代表以黄色字符为结尾 的最长的回文串，红色串代表黑色串的最长真后缀回文串，（定义：若回文串a为某串b的真后缀，则a为b的真后缀回文串，定义：若后缀a为某串b的后缀且 a！=b，则a为b的真后缀）当我们在黑串后面加上一个绿字符形成新串的时候，（姑且叫他黑绿串）回文自动机中的节点会发生什么样的变化呢？显然增加了 某些新的回文串，但是我们考虑这样一个事实，如果我们把黑绿串的最长后缀回文串加入到自动机中以后，我们就把黑绿串新增的所有回文串都可以被回文自动机所表示，证明类似于manacher算法，请自行证明。

下面来讨论如何把它加进去。我们假设粉色字符刚好是是红色串前的一个字符，如果粉色字符和绿色 字符为同一个字符的时候，我们可以肯定，黑绿串的最长真后缀回文串就是粉色字符+红串+绿色字符。此刻我们只需要新增一个节点了。如果他们不相等的话，重复 我们对黑串的算法与红串之上即可解决。

然后我们来考虑新增节点的所表示的回文后缀的最长真后缀回文串，我们重复使用上一段的算法，即可找到。

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
struct palindrome_tree{
    static const int maxn=1.2e5+5;
    int trans[maxn][26],len[maxn],suf[maxn],num[maxn];
    //len代表回文长度,suf指向回文后缀,类似于fail指针，num是最长后缀的数量，经过calc之后是后缀数量
    int last,cnt;//last是上一个回文后缀，cnt是所有节点数
    
    int new_node(int _len,int _suf,int _num){//长度，后缀，数量
        for(int i=0;i<26;i++)trans[cnt][i]=0;
        len[cnt]=_len;
        suf[cnt]=_suf;
        num[cnt]=_num;
        return cnt++;
    }
    
    void ini(){
        cnt=0;
        int root_even=new_node(0,1,0);//=1
        int root_odd=new_node(-1,1,0);//=0
        last=root_odd;
    }
    
    int query_longest_palindrom(){
        int ret=1;
        for(int i=0;i<cnt;i++){
            ret=max(ret,len[i]);
        }
        return ret;
    }
    
    void build(char*s){//s是要建立回文自动机的字符串,下标从1开始
        int len=(int)strlen(s+1);
        for(int i=1;i<=len;i++)extend(s,i);
        calc();
    }
    
    void extend(char*s,int cur){
        int w=s[cur]-'a';//当前结点的值
        int p=last;//上一次匹配到的回文后缀
    while( cur-len[p]-1<1 || s[cur-len[p]-1] != s[cur]) p=suf[p]; // BUG 数组越界
        //现在p结点的cur儿子,要么是匹配成功的最长非自身回文后缀,要么是自身这一个字符
        
        if(!trans[p][w]){//如果此回文后缀未出现过，要创建节点
            int v=suf[p];//我们来找他的suffix link回边，
            while( s[cur-len[v]-1] != s[cur] )v=suf[v];
            //此时意味着找到了suffix link 是v的儿子
            trans[p][w]=new_node(len[p]+2,trans[v][w],0);
        }
      
        last=trans[p][w];//这一次匹配到的回文后缀
        num[last]++;
    }
    
    void calc(){
        for( int i=cnt-1;~i;i-- )num[suf[i]] += num[i];//结点创建顺序保证了suf[i]<i
    }
};
ACM►学习笔记►字符串

错位比较

codeforces 1043B

本质上让你求一个字符串可能的循环节，题目的第二个样例可以看作求序列 1，2，2，1，2，的循环节显然可以是3或5， 于是算法出来了，我们只要求出最小的循环节就可以了，最小的循环节怎么来求呢，我们hash预处理原字符串以便后面O1查询区间hash值，对于字符串s如下图，我们对原字符串复制一份，错位比较红色竖线之间的子串，如果相等，意味着什么读者可以自行思考，于是我们已经找到了On的做法了，我们只要从小到大枚举 蓝色区域即可

ACM►学习笔记►字符串

后缀数组

定义：

先定义一个字符串s，假设它的长度为n，s[i]表示第i个元素 ，s[i…]代表以s[i]开头且包含s[i]的后缀。我们定义新的数组 sa[i]为一个0-n的排列，且sa[i]为后缀s[i…]在所有后缀中按 照从小到大排序的排名。最后定义rank是sa的反函数。

sa数组也称为后缀数组

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
//定义以s[i]开头的后缀是suf(i)
//后缀数组性质：suf(sa[i])<suf(sa[i+1])
//模版从下标为1的地方开始，引入的字符串下标也要从下标为1开始

struct SA{
    static const int MAXN=2e5+555;
    static int lg[MAXN];
    int h[MAXN][20],rank[MAXN],sa[MAXN],n;

    void init(char*s,int len){//第一个参数要是从下标为1开始的字符串，第二个参数是要从下标为
        static int x[MAXN],y[MAXN],c[MAXN];//全部是辅助数组
        n=len;//初始化后缀数组内部的n->s后缀数组长度
        int m=1000;//桶的大小
        for(int i=1;i<=n;i++)sa[i]=y[i]=0;//初始化sa,y
        for(int i=1;i<=n;i++)x[i]=s[i];//把s复制到x
        for(int i=1;i<=m;i++)c[i]=0;//初始化c
        for(int i=1;i<=n;i++)c[x[i]]++;//对x计数
        for(int i=1;i<=m;i++)c[i]+=c[i-1];//计数前缀和
        for(int i=n;i>=1;i--)sa[c[x[i]]--]=i;//按照计数排序后的结果
        for(int k=1;k<=n;k<<=1){
            int num=0;
            for(int i=n-k+1;i<=n;i++)y[++num]=i;
            for(int i=1;i<=n;i++)if(sa[i]>k)y[++num]=sa[i]-k;
            for(int i=1;i<=m;i++)c[i]=0;//初始化c
            for(int i=1;i<=n;i++)c[x[i]]++;//对x计数
            for(int i=1;i<=m;i++)c[i]+=c[i-1];//计数前缀和
            for(int i=n;i>=1;i--)sa[c[x[y[i]]]--]=y[i];
            for(int i=1;i<=n;i++)y[i]=x[i];
            for(int i=1;i<=n;i++)x[i]=0;
            num=1;
            x[sa[1]]=1;
            for(int i=2;i<=n;i++){
                if((y[sa[i]]!=y[sa[i-1]])||(y[sa[i]+k]!=y[sa[i-1]+k])){
                    x[sa[i]]=++num;
                }
                else x[sa[i]]=num;
            }
            if(num>=n)break;
            m=num;
        }
        for(int i=1;i<=n;i++)rank[i]=x[i];
        //获取高度数组
        int k=0;
        for(int i=1;i<=n;i++){
            if(k)k--;
            int j=sa[rank[i]-1];
            while((i+k<=n)&&(j+k<=n)&&(s[i+k]==s[j+k]))k++;
            h[rank[i]][0]=k;
        }
        //对高度数组做rmq
        for(int j=1;j<20;j++){
            int d=1<<j;
            for(int i=1;i+2*d-1<=n;i++)h[i][j]=min(h[i][j-1],h[i+d][j-1]);
        }
        if(lg[1]!=0)for(int i=1;i<MAXN;i++)lg[i]=trunc(log(i+0.5)/log(2));
    }

    int lcp(int x,int y){//注意没有下标检查，如果访问越界的话，会错。
        int L=min(rank[x],rank[y])+1;
        int R=max(rank[x],rank[y]);
        int k=lg[R-L+1];
        return min(h[L][k],h[R-(1<<k)+1][k]);
    }
}sa;
int SA::lg[SA::MAXN];
ACM►学习笔记►数学

FFT

题目1

2018宁夏网络赛H题

给你类似剪刀石头布游戏的五种手势和十种克制关系。每种手势克制其他两种手势并被另外两种手势克制。给你两个字符串分别表示A，B的 手势顺序，A长B短，你可以随意挪动B相对于A的位置，求B最多能赢多少次？

1
2
3
4
5
6
7
8
9
10
11
若两个串为abcde和ede
则共有三种方法来决斗，
abcde
eda

abcde
 eda

abcde
  eda
为每一个字符标上指数，将一个串逆置后，发现是fft

题目2

BZOJ4259

给你两个字符串原串和匹配串，串里只包含小写字母和∗,∗可以表示任意字符，问匹配串在原串不同位置中出现多少次，起始位置不同即不同位置。
我们还是利用上一题逆置短串构造卷积的方法，我们把*的值当作赋为0，pow(str1[i]−str2[i],2)*str1[i]*str2[i] 化简消去负数项即可

FFT

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
const int maxn=4e6//开四倍所有数组。准确的数值是乘法结果数组的2进制对齐那么大。这个值小于四倍，
const double pi=acos(-1.0);

struct cp{
    double x,y;
    cp(double X=0,double Y=0){x=X,y=Y;}
    cp operator + (cp b){return cp(x+b.x,y+b.y);}
    cp operator - (cp b){return cp(x-b.x,y-b.y);}
    cp operator * (cp b){return cp(x*b.x-y*b.y,x*b.y+y*b.x);}
};

// opt = 1 把系数表达式变为点值表达式
// opt =-1 把点值表达式变成系数表达式
// 传参的时候必须要求n为2的bits次方
int r[maxn];
void fft(cp *a,int n,int bits,int opt){
    for (int i=0; i<n; ++i){
        r[i]=(r[i>>1]>>1)|((i&1)<<(bits-1));
        if (i<r[i]) swap(a[i],a[r[i]]);
    }
    for (int k=1; k<n; k<<=1){
        cp wn=cp(cos(pi/k),opt*sin(pi/k));
        for (int i=0; i<n; i+=(k<<1)){
            cp w=cp(1,0);
            for (int j=0; j<k; ++j,w=w*wn){
                cp x=a[i+j],y=w*a[i+j+k];
                a[i+j]=x+y,a[i+j+k]=x-y;
            }
        }
    }
}

// 没有必要分x和y的长度,没有意义，对时间复杂度影响不大
// 传数组使用闭区间[0,n-1],n-1为多项式最高次数
// xy甚至可以传一个x或y一样的参数 例如mul(f,g,f,n,m)等价于f*=g
cp cpa[maxn],cpb[maxn];//经测试, 两个1e6的多项式乘起来大概的时间700ms, 
void mul(int*x,int*y,int*xy,int xn,int yn){
    int exn=1,bits=0;
    while(exn-1 < xn-1+yn-1)exn<<=1,bits++;//exn为xy数组拓展的长度，bits为lg(exn)
    for(int i=0 ;i<xn ;i++)cpa[i]=cp(x[i],0);
    for(int i=xn;i<exn;i++)cpa[i]=cp(0,0);//按0补齐
    for(int i=0 ;i<yn ;i++)cpb[i]=cp(y[i],0);
    for(int i=yn;i<exn;i++)cpb[i]=cp(0,0);//按0补齐
    fft(cpa,exn,bits,1);fft(cpb,exn,bits,1);//化为点值表示法
    for(int i=0; i<exn;i++)cpa[i]=cpa[i]*cpb[i];//乘
    fft(cpa,exn,bits,-1);//化为系数表示法
    for(int i=0; i<exn;i++)xy[i]=cpa[i].x/exn+0.5;//赋值结果
}

NTT

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
typedef long long ll;
// ll mod=1004535809,g=3;// 479*pow(2,21)+1
// ll mod=998244353,g=3;//  119*pow(2,23)+1
// ll mod=2281701377,g=3;//  17*pow(2,27)+1

ll qpow(ll a,ll b,ll mod){
    ll ret=1;
    while(b){
        if(b&1) ret=ret*a%mod;
        a=a*a%mod;
        b>>=1;
    }
    return ret;
}

const ll maxn=1e6;
ll r[maxn];
void ntt(ll*a,ll n,ll bits,ll opt,ll g,ll mod){
    for(ll i=0;i<n;i++) {
        r[i]=(r[i>>1]>>1)|((i&1)<<(bits-1));
        if(i<r[i]) swap(a[i],a[r[i]]);
    }
    for(ll k=1;k<n;k<<=1){
        ll wn=qpow(g,(mod-1)/(k<<1),mod);//用原根代替单位根
        if(opt==-1) wn=qpow(wn,mod-2,mod);//逆变换则乘上逆元
        for (ll i=0;i<n;i+=(k<<1)){
            for (ll j=0,w=1; j<k; j++,w=w*wn%mod){
                ll x=a[i+j], y=w*a[i+j+k]%mod;
                a[i+j]=(x+y)%mod, a[i+j+k]=(x-y+mod)%mod;
            }
        }
    }
    if(opt==-1) {
        ll inv=qpow(n,mod-2,mod);
        for(ll i=0;i<n;i++) a[i]=a[i]*inv%mod;
    }
}

ll cpx[maxn],cpy[maxn];
void mul(ll*x,ll*y,ll*xy,ll xn,ll yn,ll g,ll mod){
    ll exn=1,bits=0;
    while(exn-1 < xn-1+yn-1)exn<<=1,bits++;//exn为xy数组拓展的长度，bits为lg(exn)
    for(ll i=0 ;i<xn ;i++)cpx[i]=x[i];
    for(ll i=xn;i<exn;i++)cpx[i]=0;//按0补齐
    for(ll i=0 ;i<yn ;i++)cpy[i]=y[i];
    for(ll i=yn;i<exn;i++)cpy[i]=0;//按0补齐
    ntt(cpx,exn,bits,1,g,mod); ntt(cpy,exn,bits,1,g,mod);//化为点值表示法
    for(ll i=0; i<exn;i++)cpx[i]=cpx[i]*cpy[i]%mod;//乘
    ntt(cpx,exn,bits,-1,g,mod);//化为系数表示法
    for(ll i=0; i<exn;i++)xy[i]=cpx[i];//赋值结果
}
ACM►学习笔记►图论

二分图最小费用固定流

这是我自己给这类图取的名字

给出定义， 有五类边

  • 第一类为原点到左边的点，容量无穷大,有费用

  • 第二类为左边的点到右边的点，一对一，容量为任意常数，费用0

  • 第三类为右边的点到汇点，容量无穷大，费用0

  • 第四类为左边的点依次连i->i+1，容量无穷大，费用0

  • 第五类从右边连向左边，容量无穷大，费用0

如下图左边的图

给出问题,求当有容量的边的都满流时的最小费用

我们发现左图的最大流一定满足有容量的边是满流的，我们尝试将其转化为最小费用最大流 我们发现，如果我们对左边的图跑最大流，当左边的图有流量进入y点集的时候，他一定会进 入源点，不会流回x点集，这就很烦，我们要的不是最大流，而是指定边的满流，我们尝试“阻 止”这个过程，怎么阻止呢，我们在考虑一个z点集，从源点出发，每当y点集有流量流入汇点，就从源点流等量流量到对应的z点集，我们最后用z点集完全替代y点集的回流功能，这样以后， 最小费用最大流就成了我们要求的最小费用固定流。

再仔细想想，我们发现y点集已经没有作用了，删掉他们，最后成了上图的右边的图。 于是解决了。

对于流入汇点的流量，我们可以构造一个新的点集，每当有流量流入汇点，我们就从源点流入 新的点集，实现回流机制。