牛客竞赛数学专题班FFT与拉格朗日插值（FFT、快速多项式乘法、拉格朗日插值）

原创已于 2024-04-24 10:45:01 修改 · 938 阅读

17 ·

本内容遵循CC 4.0 BY-SA版权协议

标签

#算法

于 2024-04-24 10:16:51 首次发布

本文围绕牛客竞赛算法编程高难度练习赛展开，涵盖多道题目。如B题用生成函数构造式子并去重；Fuzzy Search题将字符拆分分析，转化为卷积形式；还涉及用拉格朗日插值法解幂和问题等，给出各题思路及代码。

(1条未读私信) 牛客竞赛_ACM/NOI/CSP/CCPC/ICPC算法编程高难度练习赛_牛客竞赛OJ (nowcoder.com)

A题 A * B Problem Plus

就是基本的fft板子题，就直接滤过了（要是要我解释FFT的话，我也听不懂你在讲什么，因为我只是个蒟蒻）

B题TSUM - Triple Sums
给定一个由 N 个不同整数组成的序列 s。
考虑序列中不同索引的三个整数的所有可能的和。
对于每个可能的和，输出生成它的不同索引三元组的数量。

我们下意识的会用生成函数构造一个这样的式子：

$(x^{s1}+x^{s2}+....x^{sn})^3$ ，嗯，似乎很合理，但是题目要求的是三个不同的索引三元组，我们这里明显有重复的比如说 $x^{2*s1+s2} or x^{3*s1}$ 之类的重复项，那么我们需要减去这些重复项；

可以构造 $(x^{2s1}+x^{2s2}+...x^{2sn})*(x^{s1}+x^{s2}+...x^{sn})$ 去重

但我们还要从三个中抽出两个去构造左边的式子，因此这个式子还要乘3；

明显又会有多减去的部分， $x^{3si}$ 就是被多减去的部分明显多减了2倍；

就有如下式子

$(x^{s1}+x^{s2}+....x^{sn})^3$ $-3*(x^{2s1}+x^{2s2}+...x^{2sn})*(x^{s1}+x^{s2}+...x^{sn})+2*(x^{3s1}+x^{3s2}+...+x^{3sn})$

但由于位置并不固定所以我们可能得到这样的三元组：

1 2 3 or 1 3 2 or 2 1 3等等。

但是毫无疑问这些三元组在题目中的意义是一样的

故整个式子还要除去 $A_{3}^{3}$ 也就是6；

令f1= $x^{s1}+x^{s2}+....x^{sn}$ ,f2= $x^{2s1}+x^{2s2}+...x^{2sn}$ ,f3= $x^{3s1}+x^{3s2}+...+x^{3sn}$

答案就是

$(f_{1}^3-3*f_{2}*f_{1}+2f_{3})/6$

代码如下：

#include<bits/stdc++.h>
using namespace std;
#define int long long
const int mod=998244353,gi[2]={3,332748118},N=1e6+10;
int q_pow(int a,int b)
{
    int t=1;
    while(b)
    {
        if(b&1)t=t*a%mod;
        b>>=1;
        a=a*a%mod;
    }
    return t;
}
void Read(int len,int fi[])
{
    int t=0;
    for(int i=0;i<len-1;i++)
    {
        int k=len/2;
        if(t>i)swap(fi[t],fi[i]);
        while(t>=k)t-=k,k/=2;
        t+=k;
    }
}
void NTT(int len,int fi[],int on)
{
    Read(len,fi);
    for(int i=2;i<=len;i<<=1)
    {
        int wi=q_pow(gi[on],(mod-1)/i);
        for(int j=0;j<len;j+=i)
        {
            int wn=1;
            for(int k=j;k<j+i/2;k++)
            {
                int t1=fi[k],t2=fi[k+i/2]*wn%mod;
                fi[k]=(t1+t2)%mod,fi[k+i/2]=(t1-t2+mod)%mod;
                wn=wn*wi%mod;
            }
        }
    }
    if(on==1)
    {
        int invlen=q_pow(len,mod-2);
        for(int i=0;i<len;i++)
        {
            fi[i]=fi[i]*invlen%mod;
        }
    }
}
int f1[N],f2[N],f3[N],n;
signed main()
{
    int k;
    int len=1;
    cin>>n;
    for(int i=0;i<n;i++)
    {
        cin>>k;
//         cout<<(k+40000)*3<<endl;
        f1[k+40000]=1;
        f2[(k+40000)*2]=1;
        f3[(k+40000)*3]=1;
       
    }
    while(len<240000)len<<=1;
    NTT(len,f1,0),NTT(len,f2,0),NTT(len,f3,0);
    for(int i=0;i<len;i++)
    {
        int t1=f1[i]*f1[i]%mod*f1[i]%mod,t2=3*f2[i]%mod*f1[i]%mod,t3=2*f3[i]%mod;
        f1[i]=(t1-t2+t3+mod)%mod*q_pow(6,mod-2)%mod;
    }
    NTT(len,f1,1);
    for(int i=0;i<N;i++)
    {
        if(f1[i])
        {
            cout<<i-120000<<" : "<<f1[i]<<endl;
        }
    }
}

第三题.Fuzzy Search
给出一个门限值 k和两个只包含 AGC 四种字符的基因串 S 和 T。现在你要找出在下列规则中 T 在 S 中出现了几次。

T 在 S的第 i 个位置中出现，当且仅当把 T 的首字符和 S 的第i个字符对齐后，T中的每一个字符能够在 S 中找到一个位置偏差不超过k的相同字符。
即对于所有的 j∈[1,∣T∣]，都存在一个 p∈[1,∣S∣] 使得 ∣(i+j−1)−p∣≤k 且 Sp=Tj 。

例如 k=1时，ACAT 出现在 AGCAATTCAT 的 2号,3号和 6号位置。 (编号从 1开始)

诈一看似乎没什么头绪，

但是如果我们把所有的字符拆开来看，一个一个分析是不是会好点。

分析题干给的 k=1 S="AGCAATTCAT",T="ACAT;

单独分析A；

对于S 我们设置两个函数f(x),g(x),若f(x)=1则表示S在第x个位置可以匹配到字符'A'；g(x)同理

那么AGCAATTCAT所对应的f(x)就是1111110111(因为S是可以模糊匹配的）

ACAT 对应的g(x)就是1010;

那么T在S在以第i个位置匹配上的个数就是 $\sum_{j=0}^{m-1}f(j+i)*g(j)$

但这并不是我们熟悉的卷积形式，如果把T转置以后那么相应的g也就倒过来了，就有一个新式子

$\sum_{j=0}^{m-1}f(j+i)*g(m-j-1)$

就变成了我们常见的卷积形式了；

要是我们把四个字母全部这样操作一次加起来的和恰好等于m，就意味着T在S的第i个位置匹配上了。

代码如下：

#include<bits/stdc++.h>
using namespace std;
#define int long long
const int mod=998244353,gi[2]={3,332748118},N=1e6+10;
int q_pow(int a,int b)
{
    int t=1;
    while(b)
    {
        if(b&1)t=t*a%mod;
        b>>=1;
        a=a*a%mod;
    }
    return t;
}
void Read(int len,int fi[])
{
    int t=0;
    for(int i=0;i<len-1;i++)
    {
        int k=len/2;
        if(t>i)swap(fi[t],fi[i]);
        while(t>=k)t-=k,k/=2;
        t+=k;
    }
}
void NTT(int len,int fi[],int on)
{
    Read(len,fi);
    for(int i=2;i<=len;i<<=1)
    {
        int wi=q_pow(gi[on],(mod-1)/i);
        for(int j=0;j<len;j+=i)
        {
            int wn=1;
            for(int k=j;k<j+i/2;k++)
            {
                int t1=fi[k],t2=fi[k+i/2]*wn%mod;
                fi[k]=(t1+t2)%mod,fi[k+i/2]=(t1-t2+mod)%mod;
                wn=wn*wi%mod;
            }
        }
    }
    if(on==1)
    {
        int invlen=q_pow(len,mod-2);
        for(int i=0;i<len;i++)
        {
            fi[i]=fi[i]*invlen%mod;
        }
    }
}
int f[4][N]={0},g[4][N]={0},sum[N],n,m,l;
string duibi="AGCT";
string s1,s2;
signed main()
{
    cin>>n>>m>>l;
    cin>>s1>>s2;
    for(int i=m-1;i>=0;i--)
    {
        if(s2[i]=='A')g[0][m-i-1]=1;
        else if(s2[i]=='G')g[1][m-i-1]=1;
        else if(s2[i]=='C')g[2][m-i-1]=1;
        else g[3][m-i-1]=1;
    }
    int k[4]={0};
    for(int i=0;i<n;i++)
    {
        for(int j=0;j<4;j++)
        {
            if(k[j])
            {
                f[j][i]=1;
                k[j]--;
            }
            if(s1[i]==duibi[j])
            {
                f[j][i]=1;
                k[j]=l;
            }
        }
    }
    for(int i=n-1;i>=0;i--)
    {
        for(int j=0;j<4;j++)
        {
            if(k[j])
            {
                f[j][i]=1;
                k[j]--;
            }
            if(s1[i]==duibi[j])
            {
                f[j][i]=1;
                k[j]=l;
            }
        }
    }
    int len=1;
    while(len<(n+m))len<<=1;
    for(int i=0;i<4;i++)
    {
        NTT(len,f[i],0);NTT(len,g[i],0);
        for(int j=0;j<len;j++)f[i][j]=f[i][j]*g[i][j]%mod;
        NTT(len,f[i],1);
        for(int j=0;j<len;j++)sum[j]+=f[i][j];
    }
    int res=0;
    for(int i=0;i<len;i++)
    {
        if(sum[i]==m)
        res++;
    }
    cout<<res<<endl;
}

第四题：The Sum of the k-th Powers

求 $\sum_{i=1}^{n}i^{k}mod(1e9+7)$

（要是大佬来看这一题，立马就是一个k+1次的多项式然后拉格朗日插值法a了。）

我之前一直挺纳闷的，为什么一定是一个k+1次的多项式呢。

我们先从k比较小的开始讨论（试图找规律的菜狗）

当k=1时就是（n+1)*n/2;

k=2时。。。。。；就要花点功夫了。。。

然后一个偶然的机会我看到一个式子

$n^2=2\binom{n}{2}+\binom{n}{1}$ ,而组合数又恰好有一个性质

$\binom{n}{m}=\binom{n-1}{m-1}+\binom{n-1}{m-2}$

那么 $\sum_{i=1}^{n}\binom{i}{k}=\binom{k}{k}+\binom{k+1}{k}+...\binom{n}{k}$ 又 $\binom{k}{k}=\binom{k+1}{k+1}$

这个式子就等于 $\binom{n+1}{k+1}$

而 $n^k$ 明显可以拆分成 $a_{0}\binom{n}{k}+a_{1}\binom{n}{k-1}+a_{2}\binom{n}{k-2}....+a_{k}\binom{n}{0}$

a0必然不为0；

求和式就变成了 $a_{0}\binom{n+1}{k+1}+a_{1}\binom{n+1}{k}+a_{2}\binom{n+1}{k-1}....+a_{k}\binom{n}{1}$

也就可以得出结论 $\sum_{i=1}^{n}i^{k}$

必然是一个n+1阶的多项式。(至于a0,a1,a2....我不知道有没有快捷的求法，实力有限）

然后就可以用拉格朗日插值法去做了，但这里也有点小技巧

就是可以直接求x=1,x=2,x=3,....x=k+1,x=k+2的值,要是要求的n小于k+1就直接输出就行；

M(i)=(-1)^(k+2-i)*(i-2)!*(k+2-i)!;

Mi(x)=(x-1)*(x-2)*....*(x-k-2)/(x-i);

求Mi(x)就是logx，所有M(i)就是O(n)

完全没问题。

答案就是 $f(n)=\sum_{i=1}^{k+2}f(i)*Mi(n)/M(i)$

代码如下：

#include<stdio.h>
#define ll long long
const int mod=1e9+7,N=1e6+10;
ll x[N],cnt[N],n,k;
ll q_pow(ll a,ll b)
{
    ll t=1;
    while(b)
    {
        if(b&1)
        t=t*a%mod;
        b>>=1;
        a=a*a%mod;
    }
    return t;
}
void slove()
{
    ll as=1,sum=0;
    cnt[0]=cnt[1]=1;
    for(int i=2;i<=k+2;i++)
        cnt[i]=cnt[i-1]*i%mod;
    for(int i=1;i<=k+2;i++)
        as=as*(n-i)%mod;
    for(int i=1;i<=k+2;i++)
    {
        int c=k-i+2,flag=1;
        if(c&1)flag=-1;
        sum=(sum+x[i]*as%mod*q_pow(cnt[i-1]*cnt[k-i+2]%mod*flag,mod-2)%mod*q_pow(n-i,mod-2)%mod)%mod;
    }
    printf("%lld",(sum+mod)%mod);
}
int main()
{
    scanf("%lld%lld",&n,&k);
    x[1]=1;
    for(int i=2;i<=k+2;i++)
    {
        x[i]=(x[i-1]+q_pow(i,k))%mod;
    }
    if(n<=k+2)
    {
        printf("%lld",x[n]);
    }
    else{
        slove();
    }
}

第五题：Polynomial

做过了第四题，第五题其实就很显然了，这个和式可以看成F(R)-F(L-1)

f(x)为n阶多项式

而F(R)明显是明显为n+1阶多项式，就要录入n+2个点，所以要先算出f(n+1)然后

求出n+2个F(i)后插值就可以了。

代码如下：

#include<bits/stdc++.h>
using namespace std;
#define int long long
const int N=1e3+10,M=1e7+10,mod=9999991;
int q_pow(int a,int b)
{
    int t=1;
    while(b)
    {
        if(b&1)t=t*a%mod;
        a=a*a%mod;
        b>>=1;
    }
    return t;
}
int inv[M];
int f[N],invfac[N];
int n,m;
void init()
{
    invfac[0]=1;
    inv[1]=inv[0]=1;
    for(int i=2;i<=mod+5;++i)
        inv[i]=(mod-mod/i)*inv[mod%i]%mod;
    for(int i=1;i<N;i++)invfac[i]=invfac[i-1]*inv[i]%mod;
}
int get(int x)
{
    if(x<=n)return f[x];
    int sum=0;
    int res=1;
    for(int i=0;i<=n;i++)
        res=res*(x-i)%mod;
    for(int i=0;i<=n;i++)
    {
        int k=res*inv[x-i]%mod;
        sum+=k*invfac[n-i]%mod*invfac[i]%mod*(((n-i)&1?-1:1)+mod)%mod*f[i]%mod;
        sum%=mod;
    }
    return sum;
}
void solve()
{
    cin>>n>>m;
    for(int i=0;i<=n;i++)
    {
        cin>>f[i];
        f[i]%=mod;
    }
    f[n+1]=get(n+1);
    n++;
    for(int i=1;i<=n;i++)
    {
        f[i]=f[i-1]+f[i];
        f[i]%=mod;
    }
    while(m--)
    {
        int l,r;
        cin>>l>>r;
        cout<<(get(r)-get(l-1)+mod)%mod<<endl;
    }
}
signed main()
{
    init();
    int t;
    cin>>t;
    while(t--)solve();
}

第六题：Assigning Prizes
构造一个长度为N的数组b(ai≤bi≤R)，并且满足bi≥bi+1(1≤i≤n−1)。
求数组b方案数，mod 1e9+7。

这个问题看起来不是很简单，实际上也不是很容易。

令f(i,d)，表示第i个位置选择bi=d并且[i,n]已经确定的方案数。

那么明显有 $f(i,d)=\sum_{j=b[i+1]}^{d}f(i+1,j)$

再分析在第n层的时候(a[n]<=x<=R)f(n,x)=1,第二层就是f(n-1,x)=i-a[n],第三层就是一次多项式的和式，那么显然(也不是那么明显)，f(i,d)其实是有关d的n-i阶多项式。

知道了，那么要怎么处理呢？

由于最后一定是一个n阶多项式，我们直接存n+1个点。

令F(i,d)= $\sum_{j=0}^{d}f(i+1,j)$

而f(i,d)=F(i,d)-F(i,b[i+1]-1)

我们直接对上一层做和，得到的就是F(i,d),通过插值的方式获取F(i,b[i+1]-1)的值，然后做差消除这些影响，就可以得到所有的f(i,d)了。

代码如下：

#include<bits/stdc++.h>
using namespace std;
#define ll long long
#define int long long
const int mod=1e9+7;
const int N=5e3+10;
ll q_pow(ll a,ll b)
{
	ll t=1;
	while(b)
	{
		if(b&1)t=t*a%mod;
		b>>=1;
		a=a*a%mod;
	}
    return t;
}
int inv[N],fnv[N];
int s1[N],s2[N];
void init()
{
	inv[0]=inv[1]=1;
	fnv[0]=1;
	for(int i=2;i<N;i++)
	{
		inv[i]=(mod-mod/i)*inv[mod%i]%mod;
	}
	for(int i=1;i<N;i++)
	{
		fnv[i]=fnv[i-1]*inv[i]%mod;
	}
}
int a[N];
int lagrange(int f[],int n,int x)
{
	int sum=0,s=1;
	s1[0]=x,s2[n+1]=1;
	for(int i=1;i<=n;i++)
	{
		s1[i]=s1[i-1]*(x-i)%mod;
	}
	for(int i=n;i>=0;i--)
	{
		s2[i]=s2[i+1]*(x-i)%mod;
	}
	for(int i=0;i<=n;i++)
	{
		int up = (ll)(i == 0 ? 1 : s1[i - 1]) * s2[i + 1] % mod;
		int down = (ll)fnv[i] * fnv[n - i] %mod* (n - i & 1 ? -1 : 1) % mod;
        sum = (sum + (ll)f[i] * down % mod* up % mod) % mod;
	
	}
    return (sum%mod+mod)%mod;
}
int l[N]={0};
void solve()
{
	int n,R;
	cin>>n>>R;
	l[0]=1;
	for(int i=1;i<=n;i++)cin>>a[i],l[i]=1;
	for(int i=n-1;i>=1;i--)
	{
		a[i]=max(a[i],a[i+1]);
	}
	for(int i=n-1;i>=0;i--)
	{
		for(int j=1;j<=n;j++)
		{
			l[j]=(l[j-1]+l[j])%mod;
		}
		int temp=lagrange(l,n-i,a[i+1]-1);
		for(int j=0;j<=n;j++)
		{
			l[j]=(l[j]-temp)%mod;
		}
	}
//	for(int i=0;i<=n;i++)cout<<l[i]<<' ';
	cout<<lagrange(l,n,R);
}
signed main()
{
	init();
	int t;
	t=1;
//	cin>>t;
	while(t--)solve();
}