POJ3415 Common Substrings

题意: 给出两个字符串S1, S2,并给出数K, 求两个字符串有多少对公共子串的长度大于等于K. S1, S2长度 <= 100000, K <= 100000, S1和S2中都是拉丁文字母. (哪些是拉丁文字母?...)

所谓多少"对"公共子串, 就是一个子串出现在S1的1次出现, 与S2中的1次出现对应, 就是1对. 比如AABABABB和ABABAAAA, 则子串AB在S1出现了3次, 在S2出现了2次, 因此由乘法原理, 子串"AB"带来的公共子串对数就是2*3= 6. 当然, 这是在K <= 2的情况下.

 

         这题比较麻烦, 首先还是考虑用后缀数组, 把两个串拼接到一块, 然后求出后缀数组和height数组. 然后怎么统计有多少对公共子串呢? 先从简单的想起...

         后缀数组sa可以认为是一个后缀的排行榜,名次->位置的映射. 每个子串都是某个后缀的前缀. 重复的前缀都会在排行榜中连续的位置出现, 因此可以从前往后遍历height数组(排行榜上相邻两后缀的公共前缀)来计算出总共有多少对公共前缀. 当然, 我们还要根据位置判定每个前缀属于S1还是S2.

         我们可以定义一个数组cnt表示"有效前缀数量",cnt[j][len] = n表示Sj有n个长度为len的子串还可以与后面的字符串匹配. 如果对于某个长为len的前缀L, 到某时, 它后面还从未出现过height值小于len的情况, 说明此时的当前后缀的前len位仍然是L, 说明L所带来的计数应该保留; 如果某个时刻出现了height值小于len, 那么后面的字符串的前len位再也不会是L了, 此时L的计数就应该去掉. 因此, 我们从前往后遍历height数组的时候, 如果height[i]不小于height[i-1], 则直接计数, 让cnt[j][height[i-1]]++ (height[i-1]表示前一位与后一位的LCP, 这个长度的前缀还可以再往后匹配); 如果height[i]小于height[i-1], 说明有些长度的前缀应该被删掉了, 此时还是先把cnt[j][height[i-1]]++, 然后把所有长度大于height[i]的计数都处理掉(肯定不能直接删掉, 因为人家对公共前缀的数量有贡献, 需要算好贡献了多少之后再删除), 因为后面的字符串和前面的至多匹配height[i]这么长.

         比如height数组局部是: (之前各个长度的有效前缀数量都是0) (K = 1, 相当于没限制)

         0, 1, 2, 2, 3, 2, 0, ...

         这里怎么计数呢? 第一个0不管, 第二个是1, 说明第一个串有1长度的有效前缀, 就把cnt[j][1]++,变成1. (j是1还是2要看这个名次对应位置在S1还是S2, 这个不难) 然后第三个是2, 比前一个大, 说明第二个串有2长度的有效前缀, 就把cnt[j][2]++; 第四, 第五个同理. 到了第6个, 比前一个少了, 此时先把前一个3统计到cnt中(因为height数组里一个3说明有两个串前缀3是相同的), 然后删除长度大于2的, 也就是长度3的. 假设那两个前缀3相同的串分别在S1与S2中, 那么删除3的时候先在最终计数里面加上1*1=1, 然后把3的计数加到2上面, 因为长度3的有效前缀的也包含长度为2的有效前缀, 可以和其他长度2的有效前缀一起匹配, 并记入最终计数当中. (最终计数就是公共前缀的对数). 最后读到0, 那么就先把前面的2统计到cnt当中, 然后全部删除+计数.

        

         唉, 解释清楚比想清楚要麻烦许多呀...

        

         然后我们就有了更为确切的方法了. 从前往后遍历height,到下标i时, 先读出sa[i-1]对应的位置属于S1还是S2, 然后把对应的cnt[][height[i-1]]++, 然后判断, 如果height[i]小于height[i-1], 那么就把cnt数组中大于height[i]的都处理掉, 处理方法就是, 对于每一个长度len, 在最终计数tot中加上cnt[1][len] * cnt[2][len], 然后把这两个值都加给len-1, 然后把这两个值清0. 一直滚雪球地清理下去, 直到达到height[i]为止.

         但还有个小问题, 就是len可能很大, 比如A和B都是abcde...一个重复程度极低的串, 然后排行榜就会长这样...

         abcd...

         abcd...$...

         bcd...

         bcd...$...

         cd...

         cd...$...

         .......

         然后height数组大概就是0, 100000,0, 99999, 0, 99998...然后如果每次删除都把len一个一个的减少, 最后减掉0的话, 那就搞成N^2的了, 那就死了...咋办呢? 注意到len的值可能很大, 但不同的len的个数不会很大, 因此可以设计某个数据结构来维护哪些len是出现过的, 以及能够很方便地进行插入, 以及删除大于某个阈值的所有len. 那就用一个堆就行了呗. 每次新得到一个len的时候, 先去cnt数组查看是否两个cnt值都是0, 不是的话说明目前还没有统计过这个len, 于是把这个len加入到堆里面, 并统计cnt. 如果不都是0, 说明已经有统计过这个len, 那就直接统计cnt. 如果要删除, 那就直接弹堆顶, 然后把cnt相乘之后, 再乘以这个堆顶与弹完之后的堆顶的差, 因为这中间的所有len值的cnt都是0, 所以如果滚雪球下来cnt值都不会变. 然后把这个乘积加到tot上面, 再把旧堆顶的cnt加给新堆顶即可. 这样就可以处理len值过大造成的问题了.

#include <iostream>
#include <algorithm>
#include <cstring>
#include <cstdio>
#include <queue>
#define N 200007
using namespace std;

int sa[N + 1], ranks[N + 1], height[N + 1];
int sz; // 后缀数组长度
int cnt[N + 1]; //基数排序的桶
char input[N + 1]; //原始数组
int K;

// int efflen[N]={};
priority_queue<int> pq;
long long cnt1[N]={};
long long cnt2[N]={};
int delim;
bool vis[N]={};

struct node {
    int v[2]; // v[0],v[1]表示二元对。v[0]是高位，v[1]是低位。
    int p; // p表示这个后缀的起始位置。
    bool operator == (const node & n) const {
        return v[0] == n.v[0] && v[1] == n.v[1];
    }
    bool operator < (const node &n) const {
        if (v[0] != n.v[0]) return v[0] < n.v[0];
        return v[1] < n.v[1];
    }
}arr[N + 1], temp[N + 1];

void print() {
    cout << "sa" << endl;
    for(int i = 1; i <= sz; i++) cout << sa[i] << " ";
    cout << endl << "rank" << endl;
    for(int i = 1; i <= sz; i++) cout << ranks[i] << " ";
    cout << endl << "height" << endl;
    for(int i = 1; i <= sz; i++) cout << height[i] << " ";
    cout << endl;
}

// 填入ranks数组。注意相同的后缀的排名应该相同
void ra() {
    for(int i = 1, j = 1, k = 1; i <= sz; i = j, k++) {
        while(j <= sz && arr[i] == arr[j]) {
            ranks[arr[j++].p] = k;
        }
    }
}

// 基数排序 桶的范围是[0, base]
void sort(int base) {
    for(int i = 1; i >= 0; i--) {
        memset(cnt, 0, sizeof(int) * (base + 1));
        for(int j = 1; j <= sz; j++) cnt[arr[j].v[i]]++;
        for(int j = 1; j <= base; j++) cnt[j] += cnt[j - 1];
        for(int j = sz; j > 0; j--) temp[cnt[arr[j].v[i]]--] = arr[j]; // 从后往前遍历，这样是稳定的
        memcpy(arr, temp, sizeof(node) * (sz + 1));
    }
    ra();
}

void buildSA() {
    for(int i = 1; i <= sz; i++) {
        arr[i].p = i;
        arr[i].v[0] = input[i];
        arr[i].v[1] = 0;
    }
    sort(arr + 1, arr + sz + 1); //第一次有负数直接排序，不影响最终复杂度
    ra();
    for(int i = 1; i < sz; i <<= 1) {
        for(int j = 1; j <= sz; j++) {
            arr[j].v[0] = ranks[j];
            arr[j].v[1] = j + i <= sz ? ranks[j + i] : 0;
            arr[j].p = j;
        }
        sort(sz);
    }
    for(int i = 1; i <= sz; i++) {
        sa[ranks[i]] = i;
    }
}

void lcp() {
    int len = 0;
    for(int i = 1; i <= sz; i++) {
        if (len > 0) len--;
        if (ranks[i] > 1) {
            while(input[i + len] == input[sa[ranks[i] - 1] + len]) len++;
        }
        height[ranks[i]] = len;
    }
}

void insHeap(int clen, int pos)	//插入长度为clen, 位置为pos的记录 
{
	if(clen < K)		return;
	if(pos > 0 & vis[pos])	return;
	vis[pos] = true;
	if(cnt1[clen] == 0 && cnt2[clen] == 0)	pq.push(clen);	//如果原来没有这个长度则插入, 避免插入两个一样的长度值 

	if(pos < delim)		cnt1[clen]++;	//用pos判断在S1还是在S2, 统计出现次数 
    if(pos >= delim)
	{
		cnt2[clen]++;
	}
}

long long delHeap(int clen)		//删除所有长度len以上的记录, 并返回计数增加的数量 
{
	long long tot = 0;
	if(clen >= K)			//删除一个堆里没有的长度clen可能造成有些长度较大的记录直接落到clen以下, 统计错误 
	{
		if(cnt1[clen] == 0 && cnt2[clen] == 0)
			pq.push(clen);
	}
	else
	{
		clen = K-1;
	}

	while(clen < pq.top())		//反复删除大于clen的 
    {
    	int oTop = pq.top();
    	pq.pop();
        int nTop = pq.top();
        tot += cnt1[oTop]*cnt2[oTop]*(oTop - nTop);		//统计该次删除长度对应的公共子串计数 
      	cnt1[nTop] += cnt1[oTop];				//把删掉的长度的计数落到下一个最大长度上 
        cnt2[nTop] += cnt2[oTop];
        cnt1[oTop] = cnt2[oTop] = 0;				//清除该次删掉长度的计数 
    }
    return tot;
}

int main(int argc, const char * argv[])
{
    while(true)
	{
        scanf("%d", &K);
        if (K == 0) break;
        
        scanf("%s", input+1);		//构造工作串: S1$S2 
        sz = strlen(input+1);
        delim = sz+1;
        input[delim] = '$';
        scanf("%s", input+sz+2);
        sz = strlen(input+1);
        
        buildSA();
        lcp();
        //print();
        
        while(!pq.empty())	pq.pop();
        pq.push(K-1);		//监视哨, 避免弹空栈RE 
        memset(cnt1, 0, sizeof(cnt1));
        memset(cnt2, 0, sizeof(cnt2));
        memset(vis, 0, sizeof(vis));	//(vis数组最后并没有用...) 
        
        long long tot = 0;
        for(int i = 1; i <= sz; i++)
        {
        	if(height[i-1] <= height[i])	//height值上升, 则当前height值就是i-1位置的有效前缀长度 
        	{
        		insHeap(height[i], sa[i-1]);
        	}
        	else			//height值下降, 则之前height值超出的部分就没用了,
        	{			//先把i-1还没有统计上的有效前缀统计上, 然后把超出当前height的长度值删掉  
        		insHeap(height[i-1], sa[i-1]);
        		int low = max(height[i], K-1);
        		tot += delHeap(low);
        	}
        }
        insHeap(height[sz], sa[sz]);	//统计上最后一位的有效前缀长度 
        tot += delHeap(K-1);		//最后清空栈, 得到最终公共子串计数 
        
        printf("%lld\n", tot);
    }
    return 0;
}