Sum of Prod of Mod of Linear_abc402G

G - Sum of Prod of Mod of Linear

这个问题需要前置知识:下取整和(floor sum),广义下取整和(generalized floor sum)。这属于类欧几里得算法,因为他们的时间复杂度和欧几里得算法一样在log级别。

下取整和的推导和结论,结论在最后一张,可以直接去看结论:

广义下取整和的推导和下取整和思路类似,但需要引入Sp(x)和sp(x)两个函数:

 

错误指正:以上写的程序的终止条件是在归约步骤,a%m=0,则程序结束,这里有误,应该是在交换步骤,当a=0时,还要计算一个常数(第一项的值),计算完常数后不再继续递归,将其作为结果返回,故判断a=0是在交换步骤进行。 

其实这样存在大量的重复计算,程序的效率还比较低,以上我们已经有的优化效率的方法有:
1.当q=0时,直接返回S(p,N)作为答案

2.当n<=0时,返回0作为答案

针对大量重复计算,可以采用记忆化的方式,发现,针对同一组n,m,a,b,当递归层数相同的时候,当前的n,m,a,b是相同的,因此记忆化可以用递归层数代替前四个参数,同时这也是因为递归层数不会很大,用memo[dep][p][q]来记忆递归到第dep层,p,q值如何时,f(n,m,a,b,p,q)的结果。

递归层数可以通过实验试探出,当n,m,a,b上限是1e6时,大概是60层,当他们上限是1e9时,会变为85层左右。

接下来给出广义下取整和 (generalized floor sum,gfs) 的程序,可以在P5170 【模板】类欧几里得算法 - 洛谷验证程序是否正确。

#include<bits/stdc++.h>
using namespace std;
using ll=long long;

const ll maxn=1e6+5,maxpq=10,mod=998244353,maxdep=85;
ll fact[maxpq+5],inv[maxpq+5],invden[maxpq+5];
ll s[maxpq+5][maxpq+5]={
    {0,1},
    {0,-1,1},
    {0,1,-3,2},
    {0,0,1,-2,1},
    {0,-1,0,10,-15,6},
    {0,0,-1,0,5,-6,2},
    {0,1,0,-7,0,21,-21,6},
    {0,0,2,0,-7,0,14,-12,3},
    {0,-3,0,20,0,-42,0,60,-45,10},
    {0,0,-3,0,10,0,-14,0,15,-10,2},
    {0,5,0,-33,0,66,0,-66,0,55,-33,6}
};
ll den[maxpq+5]={1,2,6,4,30,12,42,24,90,20,66};
ll memo[maxdep][4][4];  //由于p+q+1最大是3,故将memo的p,q维度定义为4

inline void init_memo(){
    memset(memo,-1,sizeof(memo));
}

inline ll ipow(ll a,ll n){
    ll res=1;
    a=a%mod;
    while(n){
        if(n&1) res=res*a%mod;
        a=a*a%mod;
        n>>=1;
    }
    return res;
}


inline ll com(ll n,ll m){
    if(m<0 || m>n) return 0;
    if(m==n || m==0) return 1;
    ll res=fact[n]*inv[n-m]%mod*inv[m]%mod;
    return res;
}


inline ll get_S(ll p,ll x){
    ll res=0;
    for(ll i=0;i<=p+1;i++) res=(res+(s[p][i]*ipow(x,i)%mod))%mod;
    res=res*invden[p]%mod;
    return res;
}

inline ll getinv(ll x){
    return ipow(x,mod-2);
}

inline ll gfs(ll n,ll m,ll a,ll b,ll p,ll q,ll dep){
    if(memo[dep][p][q]!=-1) {
        return memo[dep][p][q];
    }
    if(n<=0) return memo[dep][p][q]=0;
    if(q==0) return memo[dep][p][q]=get_S(p,n);
    if(a>=m || b>=m || b<0) {  //归约步骤
        ll a1=a/m,b1=b/m;
        ll a2=a%m;
        ll b2=b%m;
        ll ans=0;
        for(ll j=0;j<=q;j++){
            for(ll t=0;t<=j;t++){
                ans=(ans+com(q,j)*com(j,t)%mod*ipow(a1,t)%mod*ipow(b1,j-t)%mod*gfs(n,m,a2,b2,p+t,q-j,dep+1)%mod)%mod;
            }
        }
        memo[dep][p][q]=ans;
        return ans;
    }else { //交换步骤
        ll K=(a*(n-1)+b)/m;
        ll ans=ipow(K,q)*get_S(p,n)%mod;
        if(a<=0) return memo[dep][p][q]=ans;
        ll delta=0;
        for(ll i=0;i<=q-1;i++){
            for(ll j=0;j<=p+1;j++){
                if(s[p][j]!=0){     //因很多s[x][y]=0,故进行此优化终止递归
                    ll t=gfs(K,a,m,m-b+a-1,i,j,dep+1);
                    delta=(delta+com(q,i)*t%mod*s[p][j]%mod)%mod;
                }
            }
        }
        delta=delta*invden[p]%mod;
        ans=(ans-delta)%mod;
        memo[dep][p][q]=ans;
        return ans;
    }
    return 0;       //guard
}

int main()
{
   // freopen("D:/in.txt","r",stdin);
   // freopen("D:/out.txt","w",stdout);
    ios::sync_with_stdio(0);cin.tie(0);
    //预计算一些值
    fact[0]=1;
    for(ll i=1;i<=maxpq;i++) fact[i]=fact[i-1]*i%mod;
    for(ll i=0;i<=maxpq;i++) {
        inv[i]=ipow(fact[i],mod-2); //阶乘和阶乘的逆元用于快速计算组合数
        invden[i]=ipow(den[i],mod-2);   //经测试,预处理小s公分母的逆元能很大程度上提高效率
    }

    ll T;
    cin>>T;
    while(T--)
    {
        ll n,m,a,b,p,q;
        cin>>n>>a>>b>>m;
        n++;
        init_memo();    //以下几组n,m,a,b相同,故他们共用memo
        ll ans2=gfs(n,m,a,b,0,2,0);
        ll ans3=gfs(n,m,a,b,1,1,0);
        ll ans1=gfs(n,m,a,b,0,1,0);

        cout<<(ans1+mod)%mod<<" ";
        cout<<(ans2+mod)%mod<<" ";
        cout<<(ans3+mod)%mod<<" ";
        cout<<"\n";
    }

    return 0;
}

标题题目解析

建议先去看官方题解,因为给出的代码实现是基于官方题解的,以下的题解思路与官方完全相同。

 

在实现过程中,需要注意:
1.对同一组(n,m,a,b)使用gfs的时候,共用memo以提高效率。

2.这里的b可能出现为负数的情况,需要略微修改归约步骤来应对b<0的情况,注意:C++负数取模的方法。

3.中间运算过程可能超过long long 范围,使用__int128进行计算,但最终结果在long long范围内,可以选择将__int128转成long long输出。

以下的实现是基于官方题解的,和上面所写实现方式略有不同,但思路一样。

#include<bits/stdc++.h>
using namespace std;
using ll=__int128;
using ull=unsigned long long ;

const ll maxn=1e6+5,maxpq=10,maxdep=60;
ll fact[maxpq+5];
ll s[maxpq+5][maxpq+5]={
    {0,1},
    {0,-1,1},
    {0,1,-3,2},
    {0,0,1,-2,1},
    {0,-1,0,10,-15,6},
    {0,0,-1,0,5,-6,2},
    {0,1,0,-7,0,21,-21,6},
    {0,0,2,0,-7,0,14,-12,3},
    {0,-3,0,20,0,-42,0,60,-45,10},
    {0,0,-3,0,10,0,-14,0,15,-10,2},
    {0,5,0,-33,0,66,0,-66,0,55,-33,6}
};
ll den[maxpq+5]={1,2,6,4,30,12,42,24,90,20,66};

ll memo[maxdep][4][4];

inline void init_memo(){
    memset(memo,-1,sizeof(memo));
}

inline ll ipow(ll a,ll n){
    ll res=1;
    while(n){
        if(n&1) res=res*a;
        a=a*a;
        n>>=1;
    }
    return res;
}

inline ll com(ll n,ll m){
    if(m<0 || m>n) return 0;
    if(m==n || m==0) return 1;
    ll res=fact[n]/(fact[n-m]*fact[m]);
    return res;
}


inline ll get_S(ll p,ll x){
    ll res=0;
    for(ll i=0;i<=p+1;i++) res+=s[p][i]*ipow(x,i);
    res/=den[p];
    return res;
}

inline ll gfs(ll n,ll m,ll a,ll b,ll p,ll q,ll dep){
    if(memo[dep][p][q]!=-1) {
        return memo[dep][p][q];
    }
    if(n<=0) return memo[dep][p][q]=0;
    if(q==0) return memo[dep][p][q]=get_S(p,n);
    if(a>=m || b>=m || b<0) {  //归约步骤
        ll a1=a/m,b1=b/m;
        if(b<0 && b%m!=0) b1--;
        ll a2=a%m;
        ll b2=(b%m+m)%m;
        ll ans=0;
        for(ll j=0;j<=q;j++){
            for(ll t=0;t<=j;t++){
                ans+=com(q,j)*com(j,t)*ipow(a1,t)*ipow(b1,j-t)*gfs(n,m,a2,b2,p+t,q-j,dep+1);
            }
        }
        memo[dep][p][q]=ans;
        return ans;
    }else { //交换步骤
        ll K=(a*(n-1)+b)/m;
        ll ans=ipow(K,q)*get_S(p,n);
        if(a<=0) return memo[dep][p][q]=ans;
        ll delta=0;
        for(ll i=0;i<=q-1;i++){
            for(ll j=0;j<=p+1;j++){
                if(s[p][j]!=0){
                    ll t=gfs(K,a,m,m-b+a-1,i,j,dep+1);
                    delta+=com(q,i)*t*s[p][j];
                }
            }
        }
        delta/=den[p];
        ans-=delta;
        memo[dep][p][q]=ans;
        return ans;
    }
    return 0;       //guard
}

istream& operator >> (istream &in,ll &x) {
    long long t;
    in>>t;
    x=t;
    return in;
}

ostream& operator << (ostream &out, const ll &x) {
        long long t=x;
        out<<t;
        return out;
}

int main()
{
    ios::sync_with_stdio(0);cin.tie(0);
    //预计算一些值
    fact[0]=1;
    for(ll i=1;i<=maxpq;i++) fact[i]=fact[i-1]*i;

    ll T;cin>>T;
    while(T--)
    {
        ll n,m,a,b1,b2;
        cin>>n>>m>>a>>b1>>b2;
        if(b1>b2) swap(b1,b2);
        if(a==0) {
            cout<<n*b1*b2<<"\n";
            continue;
        }
        ll res1=a*a*(n-1)*n/2*(2*n-1)/3+a*b2*n*(n-1)/2+a*b1*n*(n-1)/2+b1*b2*n;

        init_memo();
        ll res2=a*gfs(n,m,a,b2,1,1,0)+b1*gfs(n,m,a,b2,0,1,0);
        ll res3=gfs(n,m,a,b2,0,2,0);

        init_memo();
        res2+=a*gfs(n,m,a,b1,1,1,0)+b2*gfs(n,m,a,b1,0,1,0);
        res2*=m;
        ll delta;
        ll X=(a*(n-1)+b2)/m;
        ll dX_1=min((X*m-b1+a-1)/a,n);
        ll dX_2=min((X*m-b2+a-1)/a,n);

        init_memo();
        delta=-gfs(X,a,m,-b1+a-1,1,1,0);
        init_memo();
        delta+=gfs(X,a,m,-b2+a-1,1,1,0)-X*(dX_1-dX_2);
        res3+=delta;
        res3*=m*m;
        ll ans=res1-res2+res3;
        cout<<ans<<"\n";
    }

    return 0;
}

 

 

评论
成就一亿技术人!
拼手气红包6.0元
还能输入1000个字符
 
红包 添加红包
表情包 插入表情
 条评论被折叠 查看
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值