蜜月日记 Part9

littlez_meow · 2025-01-14 14:37:03 · 休闲·娱乐

前言

争取写到 Day 100.

Day 81

数论已经好久没有在蜜月日记里出现了。

2025/01/14 [AT Xmas Contest 2019] Sum of (-1)^f(n)

题意

定义 f(x) 是 x 的质因数分解后各个质因子的次数和，求 \sum\limits_{i=1}^n (-1)^{f(i)}。

思路

显然 f(x) 是完全和性函数，故 (-1)^{f(x)} 是（完全）积性函数。下面设 F(x)=(-1)^{f(x)}。

可以直接 min_25，但是我不会。又懒得去推这个函数如何杜教筛。

所以考虑一些科技飞升做法，比如 PN 筛。下面会给出 PN 筛的流程，但不会有定义和证明。

首先寻找质数 p 处取值相同的积性函数 G。我们有 (-1)^{f(p)}=-1=\mu(p)，所以取 G(x)=\mu(x)。如无特殊说明，下面 p 均为质数。

然后设积性函数 H(x) 满足 F=G*H，这里 * 是狄利克雷卷积，则 H 只在 PN 处不为 0。

则我们有 F(p^k)=\sum\limits_{i=0}^k \mu(p^i)H(p^{k-i})。简单化简得到 (-1)^k=H(p^k)-H(p^{k-1})。结合 H(p)=0，有 H(p^k)=[2|k]。

我们有 \sum\limits_{i=1}^n F(i)=\sum\limits_{i=1}^n\sum\limits_{d|i}H(d)\mu(\dfrac{i}{d})。

枚举因数变枚举倍数得到 \sum\limits_{i=1}^n F(i)=\sum\limits_{d=1}^n\sum\limits_{i=1}^{\lfloor\frac{n}{d}\rfloor}\mu(i)H(d)。

设 S(x)=\sum\limits_{i=1}^x\mu(i)，化简得到 \sum\limits_{i=1}^n F(i)=\sum\limits_{d=1}^n H(d)S(\lfloor\dfrac{n}{d}\rfloor)。

枚举令 H(d)\neq 0 的 O(\sqrt n) 个 PN，用一次杜教筛算出 O(\sqrt n) 个本质不同的 S(\lfloor\dfrac{n}{d}\rfloor) 即可。

瓶颈在于一次杜教筛，时间复杂度 O(n^{\frac 2 3})。

代码

#include<bits/stdc++.h>
#include<ext/pb_ds/assoc_container.hpp>
#include<ext/pb_ds/hash_policy.hpp>
#define F(i,a,b) for(int i(a),i##i##end(b);i<=i##i##end;++i)
#define R(i,a,b) for(int i(a),i##i##end(b);i>=i##i##end;--i)
#define ll long long
#define File(a) freopen(#a".in","r",stdin);freopen(#a".out","w",stdout)
using namespace std;
const int MAXN=5e7+1,LIM=1e6+1;
bool vis[MAXN];
ll mu[MAXN];
int pr[3001135],cnt;
inline void init(){
    mu[1]=1;
    F(i,2,MAXN-1){
        if(!vis[i]) pr[++cnt]=i,mu[i]=-1;
        F(j,1,cnt){
            ll k=i*pr[j];
            if(k>=MAXN) break;
            vis[k]=1;
            if(i%pr[j]) mu[k]=-mu[i];
            else break;
        }
    }
    F(i,2,MAXN-1) mu[i]+=mu[i-1];
    return;
}
__gnu_pbds::gp_hash_table<ll,int>res;
ll n,ans;
ll S(ll x){
    if(x<MAXN&&x!=n) return mu[x];
    if(res.find(x)!=res.end()) return res[x];
    ll ans=1;
    for(ll l=2,r;l<=x;l=r+1){
        ll qwq=x/l;
        r=x/qwq;
        ans=ans-(r-l+1)*S(qwq);
    }
    return res[x]=ans;
}
void dfs(int step,ll now,ll H){
    ll pri=pr[step],nxt=pri*pri;
    if(step>cnt||now>n/nxt){
        ans+=H*S(n/now);
        return;
    }
    dfs(step+1,now,H);
    for(int expo(2);now*nxt<=n;++expo,nxt*=pri) dfs(step+1,now*nxt,H*((expo&1)==0));
    return;
}
signed main(){
    ios::sync_with_stdio(0);
    cin.tie(0),cout.tie(0);
    cin>>n;
    init(),S(n);
    dfs(1,1,1);
    cout<<ans;
    return 0;
}

Day 82

快要去 WC 了，浅卷一卷。

2025/01/16 [WC2021] 斐波那契

题意

定义数列 F_0=a,F_1=b,F_n=(F_{n-1}+F_{n-2})\bmod m，N 次询问 a,b 求满足 F_p=0 的最小非负整数 p。

思路

注意区分下文的 N,n。

设斐波那契数列为 \{f_n\}，则 F_n=af_{n-1}+bf_n。相当于要求最小的 n 满足 af_{n-1}+bf_n\equiv 0\pmod m。

关于模意义下的斐波那契数列，有经典结论：

【皮萨诺周期】 模 m 意义下斐波那契数列有长度 \le 6m 的最小正周期，称为皮萨诺周期。具体地，有如下几个结论计算：

• 若 \gcd(a,b)=1，则 ab 的皮萨诺周期是 a,b 的皮萨诺周期的最小公倍数。

• 对于奇质数 p\equiv 1,4\pmod 5，皮萨诺周期为 p-1。

• 对于奇质数 p\equiv 2,3\pmod 5，皮萨诺周期为 2p+2。

• 若 p 为质数，则 p^k 的皮萨诺周期为 p 的皮萨诺周期乘 p^{k-1}。

本题中只需要用到皮萨诺周期是 O(m) 的这个结论。

先特判 a=0 或 b=0，下面认为 a,b 全不为 0。

令 b\leftarrow m-b，则要求最小的 n 满足 af_{n-1}\equiv bf_n\pmod m。

设 a'=\dfrac{a}{\gcd(a,b,m)},b'=\dfrac{b}{\gcd(a,b,m)},m'=\dfrac{m}{\gcd(a,b,m)}，根据同余的基本性质，所求式变为 a'f_{n-1}\equiv b'f_n\pmod{m'}。注意这里 a',b',m' 并不是两两互质，不一定存在逆元。

根据另一个经典结论，\gcd(f_n,f_{n-1})=f_{\gcd(n,n-1)}=1，因此有 \gcd(f_{n-1},m')=\gcd(b',m')=r,\gcd(f_n,m')=\gcd(a',m')=s。

再设 a''=\dfrac{a'}{r},b''=\dfrac{b'}{s},m''=\dfrac{m'}{rs}，则 a''f_{n-1}\equiv b''f_n\pmod{m''}。此时它们两两互质，故可以求逆元。对于每个 m'|m 预处理三元组 (r,s,\dfrac{f_{n}}{f_{n-1}}\bmod{\dfrac{m'}{rs}})，用哈希表存一下就行。

懒了，用的 map，时间复杂度 O(\sigma(m)\log m+N\log(\sigma(m)))，其中 \sigma 是约数和函数。

代码

#include<bits/stdc++.h>
#define File(a) freopen(#a".in","r",stdin);freopen(#a".out","w",stdout)
#define ll long long
#define F(i,a,b) for(int i(a),i##i##end(b);i<=i##i##end;++i)
#define R(i,a,b) for(int i(a),i##i##end(b);i>=i##i##end;--i)
using namespace std;
const int MAXN=1e5+1;
void exgcd(int a,int b,int&x,int&y){
    if(b==0) return x=1,y=0,void();
    else return exgcd(b,a%b,y,x),y-=a/b*x,void();
}
ll inv(int x,int MOD){
    int inv,y;
    exgcd(x,MOD,inv,y);
    return (inv%MOD+MOD)%MOD;
}
int N,m;
struct Triple{
    int r,s,val;
    bool operator<(const Triple&x)const{
        return r==x.r?(s==x.s?val<x.val:s<x.s):r<x.r;
    }
};
map<Triple,int>ans[MAXN];
signed main(){
    ios::sync_with_stdio(0);
    cin.tie(0),cout.tie(0);
    cin>>N>>m;
    F(mm,2,m) if(m%mm==0){
        int fn1=0,fn=1;
        for(int n=1;;++n){
            int r=__gcd(fn1,mm),s=__gcd(fn,mm),MOD=mm/r/s;
            Triple qwq={r,s,fn/s*inv(fn1/r,MOD)%MOD};
            if(ans[mm].find(qwq)==ans[mm].end()) ans[mm][qwq]=n;
            int nxt=fn+fn1;
            nxt>=mm&&(nxt-=mm);
            fn1=fn,fn=nxt;
            if(fn1==0&&fn==1) break;
        }
    }
    for(;N;--N){
        int a,b,mm,d;
        cin>>a>>b,b=(m-b)%m;
        if(a==0){
            cout<<"0\n";
            continue;
        }
        if(b==0){
            cout<<"1\n";
            continue;
        }
        d=__gcd(__gcd(a,b),m);
        a/=d,b/=d,mm=m/d;
        int r=__gcd(a,mm),s=__gcd(b,mm),MOD=mm/r/s;
        a=a/r,b=b/s;
        auto it=ans[mm].find((Triple){s,r,a*inv(b,MOD)%MOD});
        if(it==ans[mm].end()) cout<<"-1\n";
        else cout<<it->second<<"\n";
    }
    return 0;
}

Day 83

WC 讲课的题，老早就会了，一直没写。

2025/01/18 [THUPC2019] 找树

题意

定义 \otimes_1, \otimes_2, \otimes_3 分别为按位与、按位或、按位异或运算。给定长为 o_1,o_2,\cdots,o_w，定义位运算 x \oplus y=z 表示二进制下 z 的从低到高第 i 位等于 x,y 的从低到高第 i 位进行 \otimes_{o_i} 运算。给定一个 n 个点 m 条边的无向图，边权 \in[0,2^w)\cap\mathbb Z，求运算 \oplus 意义下的最大生成树。

思路

数据范围超级小，那就应该不是传统最优化。

这个运算没有什么性质，但是众所周知的是我们可以用生成函数求边权和为某个数的生成树个数，这里同样可以。

设一条边权为 v 的边的生成函数为 x^v，多项式乘法定义为 \oplus 卷积，加法就是普通多项式加法，然后就可以跑矩阵树定理了。

又因为 FWT 是线性变换，可以在求行列式之前对矩阵中每个元素求 FWT，再对每一位求行列式，最后 IFWT 回来。

这里 \oplus 的 FWT 可以视为对每一维做不同变换，即对按位与的位做长度为 2 的 \min 卷积的变换（后缀和），按位或的为做长度为 2 的 \max 卷积的变换（前缀和），按位异或做长度为 2 的加法卷积的变换（FFT）。

方案数可以对大质数取模。

时间复杂度 O(n^3 2^w)。

代码

#include<bits/stdc++.h>
#define File(a) freopen(#a".in","r",stdin);freopen(#a".out","w",stdout)
#define ll long long
#define F(i,a,b) for(int i(a),i##i##end(b);i<=i##i##end;++i)
#define R(i,a,b) for(int i(a),i##i##end(b);i>=i##i##end;--i)
using namespace std;
const int MAXN=71,MAXL=4096,MOD=1e9+9,INV2=(MOD+1)/2;
inline ll qpow(ll base,int expo){
    ll res(1);
    while(expo){
        (expo&1)&&(res=res*base%MOD);
        base=base*base%MOD,expo>>=1;
    }
    return res;
}
int n,m,w,l;
char o[MAXN];
ll kir[MAXN][MAXN][MAXL];
inline void FWT(ll*poly,bool inv){
    for(int len=2,i=0;len<=l;len<<=1,++i){
        int qwq=len>>1;
        switch(o[i]){
            case '&':{
                for(int i(0);i<l;i+=len) F(j,i,i+qwq-1){
                    if(inv) (poly[j]-=poly[j+qwq])<0&&(poly[j]+=MOD);
                    else (poly[j]+=poly[j+qwq])>=MOD&&(poly[j]-=MOD);
                }
                break;
            }
            case '|':{
                for(int i(0);i<l;i+=len) F(j,i,i+qwq-1){
                    if(inv) (poly[j+qwq]-=poly[j])<0&&(poly[j+qwq]+=MOD);
                    else (poly[j+qwq]+=poly[j])>=MOD&&(poly[j+qwq]-=MOD);
                }
                break;
            }
            case '^':{
                for(int i(0);i<l;i+=len) F(j,i,i+qwq-1){
                    int x=poly[j],y=poly[j+qwq];
                    (poly[j]=x+y)>=MOD&&(poly[j]-=MOD);
                    (poly[j+qwq]=x-y)<0&&(poly[j+qwq]+=MOD);
                    if(inv) poly[j]=poly[j]*INV2%MOD,poly[j+qwq]=poly[j+qwq]*INV2%MOD;
                }
                break;
            }
        }
    }
    return;
}
ll mat[MAXN][MAXN],ans[MAXL];
inline ll det(){
    bool neg=0;
    F(i,1,n-1){
        if(!mat[i][i]){
            bool flag=1;
            F(j,i+1,n-1) if(mat[j][i]){
                flag=0;
                swap(mat[j],mat[i]);
                neg^=1;
                break;
            }
            if(flag) return 0;
        }
        ll inv=qpow(mat[i][i],MOD-2);
        F(j,i+1,n-1){
            ll factor=inv*mat[j][i]%MOD;
            F(k,i,n-1) (mat[j][k]=mat[j][k]-factor*mat[i][k]%MOD)<0&&(mat[j][k]+=MOD);
        }
    }
    ll res=(neg?MOD-1:1);
    F(i,1,n-1) res=res*mat[i][i]%MOD;
    return res;
}
signed main(){
    ios::sync_with_stdio(0);
    cin.tie(0),cout.tie(0);
    cin>>n>>m>>o;
    w=strlen(o),l=1<<w;
    F(i,1,m){
        int x,y,v;
        cin>>x>>y>>v;
        --kir[x][y][v],--kir[y][x][v];
        ++kir[x][x][v],++kir[y][y][v];
    }
    F(x,1,n-1) F(y,1,n-1){
        F(v,0,l-1) kir[x][y][v]<0&&(kir[x][y][v]+=MOD);
        FWT(kir[x][y],0);
    }
    F(i,0,l-1){
        F(x,1,n-1) F(y,1,n-1) mat[x][y]=kir[x][y][i];
        ans[i]=det();
    }
    FWT(ans,1);
    R(i,l-1,0) if(ans[i]) return cout<<i,0;
    cout<<"-1";
    return 0;
}

Day 84

考 WC 前紧急打一打多项式板子。

2025/01/19 [ABC385G] Counting Buildings

题意

求长为 n 的前缀最大值个数减后缀最大值个数为 k 的排列个数模 998244353。

思路

把排列从最大值处分开，右侧没有前缀最大值，左侧没有后缀最大值。

设前后缀缀最大值分别有 x,y 个。先考虑最大值左边，发现可以当成是把左边的元素分成了 x-1 个排列，强制要求最大值在第一个。这又相当于把左边的元素分成了 x-1 个轮换。

右边同理，是 y-1 个轮换。方案数就是把 n-1 个元素分成 x+y-2 个轮换，选出 x-1 个分配到左边的方案数。前半部分是第一类斯特林数，后半部分是组合数。即此时的方案数为 s(n-1,x+y-2)\dbinom{x+y-2}{x-1}。枚举 y，有 x=y+k 即可计算答案。

问题变成求一行第一类斯特林数。众所周知，第 n 行第一类斯特林数的生成函数是 x^{\bar n}=x(x+1)\cdots(x+n-1)，可以用倍增做到 O(n\log n)。不过我这里是打板子，写的分治 NTT，时间复杂度 O(n\log^2 n)。

代码

#include<bits/stdc++.h>
#define File(a) freopen(#a".in","r",stdin);freopen(#a".out","w",stdout)
#define ll long long
#define F(i,a,b) for(int i(a),i##i##end(b);i<=i##i##end;++i)
#define R(i,a,b) for(int i(a),i##i##end(b);i>=i##i##end;--i)
#define vi vector<int>
using namespace std;
const int MAXN=(1<<19)+5,MOD=998244353,G=3;
inline ll qpow(ll base,int expo){
    ll res(1);
    while(expo){
        (expo&1)&&(res=res*base%MOD);
        base=base*base%MOD,expo>>=1;
    }
    return res;
}
const int INVG=qpow(G,MOD-2);
ll powg[2][20][MAXN];
inline void init(){
    F(i,0,1) F(j,1,19){
        powg[i][j][0]=1;
        ll qwq=qpow(i?INVG:G,(MOD-1)>>j);
        F(k,1,MAXN-1) powg[i][j][k]=powg[i][j][k-1]*qwq%MOD;
    }
    return;
}
int rev[MAXN];
inline void NTT(int*poly,int len,bool inv){
    F(i,0,len-1) if(i<rev[i]) swap(poly[i],poly[rev[i]]);
    for(int i=1,expo=1;i<len;i<<=1,++expo) for(int j=0;j<len;j+=(i<<1)) F(k,0,i-1){
        int&x(poly[j|k]),&y(poly[i|j|k]);
        ll qwq=y*powg[inv][expo][k]%MOD;
        (y=x-qwq)<0&&(y+=MOD);
        (x+=qwq)>=MOD&&(x-=MOD);
    }
    if(inv){
        ll invl=qpow(len,MOD-2);
        F(i,0,len-1) poly[i]=poly[i]*invl%MOD;
    } 
    return;
}
int px[MAXN],py[MAXN];
inline vi mul(const vi&x,const vi&y){
    int n=x.size()+y.size()-1,l=1;
    while(l<n) l<<=1;
    F(i,0,l-1) rev[i]=(rev[i>>1]>>1)|((i&1)?(l>>1):0);
    F(i,0,x.size()-1) px[i]=x[i];
    F(i,x.size(),l-1) px[i]=0;
    F(i,0,y.size()-1) py[i]=y[i];
    F(i,y.size(),l-1) py[i]=0;
    NTT(px,l,0),NTT(py,l,0);
    F(i,0,l-1) px[i]=px[i]*1ll*py[i]%MOD;
    NTT(px,l,1);
    vi res(n);
    F(i,0,n-1) res[i]=px[i];
    return res;
}
vi solve(int l,int r){
    if(l>r) return vi({1});
    if(l==r) return vi({l,1});
    int mid=(l+r)>>1;
    return mul(solve(l,mid),solve(mid+1,r));
}
int n,k;
vi s;
ll fact[MAXN],inv[MAXN];
inline ll C(int x,int y){
    return x>=y&&y>=0?fact[x]*inv[y]%MOD*inv[x-y]%MOD:0;
}
signed main(){
    ios::sync_with_stdio(0);
    cin.tie(0),cout.tie(0);
    init();
    cin>>n>>k;
    s=solve(0,n-2);
    fact[0]=1;
    F(i,1,n) fact[i]=fact[i-1]*i%MOD;
    inv[n]=qpow(fact[n],MOD-2);
    R(i,n,1) inv[i-1]=inv[i]*i%MOD;
    ll ans=0;
    F(y,1,n){
        int x=y+k;
        if(x<=0) continue;
        if(x+y-2>n-1) break;
        ans=(ans+s[x+y-2]*C(x+y-2,x-1))%MOD;
    }
    cout<<ans;
    return 0;
}

Day 85

被 WC T2 位运算创飞了。感觉自己不会位运算相关 trick，尤其是既有异或又有加法的。

2025/01/21 [HNOI/AHOI2018] 寻宝游戏

题意

给定 n+q 个 m 位二进制数 a_1,a_2,\cdots,a_n,r_1,r_2,\cdots,r_q，对每个 t=1,2,\cdots,q，求出运算符序列 \otimes_1,\cdots,\otimes_n 的个数，满足 \otimes_i 为按位且和按位或之一，且 0\otimes_1 a_1\otimes_2\cdots\otimes_n a_n=r_t。答案对 10^9+7 取模。

思路

首先这个肯定是按位考虑。下面暂时认为 m=1。

先研究 \otimes_n 应该怎么取。如果是按位且 1，没有任何效果，就要以相同方式研究 \otimes_{n-1}；如果是按位且 0，无论前面如何，结果都是 0。按位或是类似的。

发现这个过程很像字典序比较。把按位或当成 0，按位且当成 1，以第 n 位为最高位，得到的二进制数设为 x。把 a_n 当成最高位所有 a_i 组成的二进制数设为 y。则有 x<y\Leftrightarrow 结果是 1。

然后是 m\neq 1 的情况，把每一位对应的不等式交起来就好了。

时间复杂度 O(nm\log n+nq)。

代码

#include<bits/stdc++.h>
#define File(a) freopen(#a".in","r",stdin);freopen(#a".out","w",stdout)
#define ll long long
#define F(i,a,b) for(int i(a),i##i##end(b);i<=i##i##end;++i)
#define R(i,a,b) for(int i(a),i##i##end(b);i>=i##i##end;--i)
#define fi first
#define se second
using namespace std;
const int MAXN=5002,MOD=1e9+7;
int n,m,q,val[MAXN];
pair<string,int>id[MAXN];
char ch[MAXN][MAXN],ask[MAXN];
signed main(){
    ios::sync_with_stdio(0);
    cin.tie(0),cout.tie(0);
    cin>>n>>m>>q;
    F(i,1,n) cin>>(ch[i]+1);
    F(i,1,m){
        R(j,n,1) id[i].fi+=ch[j][i],val[i]=((val[i]<<1)|(ch[j][i]^'0'))%MOD;
        id[i].se=i;
    }
    sort(id+1,id+m+1);
    F(i,1,n) id[0].fi+='0',id[m+1].fi+='1',val[m+1]=((val[m+1]<<1)|1)%MOD;
    id[0].se=0,id[m+1].se=m+1,val[0]=0,++val[m+1];
    for(;q;--q){
        cin>>(ask+1);
        int l=m+1,r=0;
        F(i,1,m){
            if(ask[id[i].se]=='0') r=max(i,r);
            else l=min(i,l); 
        }
        if(l<r) cout<<"0\n";
        else cout<<(val[id[l].se]-val[id[r].se]+MOD)%MOD<<"\n";
    }
    return 0;
}

Day 86

寒假摆会烂，写点典题。

2025/01/25 无标号无根树计数

题意

求 n 个点无标号无根树的数量。

思路

首先求出 n 个点无标号有根树的数量 f_n，再钦定重心为根。

设 f_n 的 OGF 为 F(x)。由于把根删掉之后是无标号有根树森林，即对无标号有根树进行 \operatorname{MSET} 构造。

在无标号体系下的 \operatorname{MSET} 构造定义为：

【欧拉变换】 \mathcal E(F(x))=\prod\limits_{i=1}^{\infty}\dfrac{1}{(1-x^i)^{[x^i]F(x)}}=\exp(\sum\limits_{i=1}^\infty\dfrac{F(x^i)}{i})，其中 \mathcal E 被称为欧拉变换或 Pólya exp，即无标号体系下的 \operatorname{MSET} 构造。

第一个等号的组合意义是背包，第二个等号是先求 \ln 再求 \exp。

于是有生成函数方程：

下面就是要化简成好处理的方式了。先展开成 F(x)=xe^{\sum\limits_{i=1}^{\infty}\frac{F(x^i)}{i}}。

两边同时求导得到 F'(x)=\mathcal E(F(x))+x\left(e^{\sum\limits_{i=1}^{\infty}\frac{F(x^i)}{i}}\right)'

使用链式求导法则，得到 F'(x)=\mathcal E(F(x))+x\exp(\sum\limits_{i=1}^{\infty}\dfrac{F(x^i)}{i})\sum\limits_{i=1}^{\infty}F'(x^i)x^{i-1}。

整理得到 F'(x)=\dfrac{F(x)}{x}+F(x)\sum\limits_{i=1}^{\infty}F'(x^i)x^{i-1}

两边同时乘 x 有 xF'(x)=F(x)\sum\limits_{i=1}^{\infty}F'(x^i)x^i。

设 G(x)=\sum\limits_{i=1}^{\infty}F'(x^i)x^i,g_n=[x^n]G(x)，不难推出 g_n=\sum\limits_{d|n}df_d。则有 nf_n=f_n+\sum\limits_{i=1}^{n-1}f_ig_{n-i}。

整理得到：

其中 f_1=1。

这个式子就可以分治 FFT 了。注意 l=1 时直接计算贡献有问题，需要特殊处理。

然后容斥掉根不是重心的方案数。

当 2\nmid n 时，只有一个重心，枚举大小 \ge \lfloor\dfrac n 2\rfloor 的子树大小，得到答案为 f_n-\sum\limits_{i=\lfloor\frac n 2\rfloor}^{n-1}f_if_{n-i}。

否则，还要减去以两个重心为根不同构的情况 \dbinom{f_{\frac n 2}}{2}。

最后答案为：

瓶颈在于分治 FFT，时间复杂度 O(n\log^2 n)。其实可以牛顿迭代但是常数过大。

代码

#include<bits/stdc++.h>
#define F(i,a,b) for(int i(a),i##i##end(b);i<=i##i##end;++i)
#define R(i,a,b) for(int i(a),i##i##end(b);i>=i##i##end;--i)
#define ll long long
#define File(a) freopen(#a".in","r",stdin);freopen(#a".out","w",stdout)
using namespace std;
const int MAXN=(1<<19)+5,MOD=998244353,G=3;
inline ll qpow(ll base,int expo){
    ll res(1);
    while(expo){
        (expo&1)&&(res=res*base%MOD);
        base=base*base%MOD,expo>>=1;
    }
    return res;
}
const int INVG=qpow(G,MOD-2);
ll powg[2][20][MAXN];
inline void init(){
    F(i,0,1) F(j,1,19){
        powg[i][j][0]=1;
        ll qwq=qpow(i?INVG:G,(MOD-1)>>j);
        F(k,1,MAXN-1) powg[i][j][k]=powg[i][j][k-1]*qwq%MOD;
    }
    return;
}
int rev[MAXN];
inline void NTT(int*poly,int len,bool inv){
    F(i,0,len-1) if(i<rev[i]) swap(poly[i],poly[rev[i]]);
    for(int i=1,expo=1;i<len;i<<=1,++expo) for(int j=0;j<len;j+=(i<<1)) F(k,0,i-1){
        int&x(poly[j|k]),&y(poly[i|j|k]);
        ll qwq=y*powg[inv][expo][k]%MOD;
        (y=x-qwq)<0&&(y+=MOD);
        (x+=qwq)>=MOD&&(x-=MOD);
    }
    if(inv){
        ll invl=qpow(len,MOD-2);
        F(i,0,len-1) poly[i]=poly[i]*invl%MOD;
    } 
    return;
}
inline int mul(int l,int*x,int*y,int*res){
    l=1<<(__lg(l)+1);
    F(i,0,l-1) rev[i]=(rev[i>>1]>>1)|((i&1)?(l>>1):0);
    static int px[MAXN],py[MAXN];
    memcpy(px,x,sizeof(int)*l);
    memcpy(py,y,sizeof(int)*l);
    NTT(px,l,0),NTT(py,l,0);
    F(i,0,l-1) res[i]=px[i]*1ll*py[i]%MOD,px[i]=py[i]=0;
    NTT(res,l,1);
    return l;
}
int n,f[MAXN],g[MAXN];
void DandQ_FFT(int l,int r){
    int mid=(l+r)>>1;
    if(l==r){
        if(l==1) f[l]=1;
        else f[l]=f[l]*qpow(l-1,MOD-2)%MOD;
        int qwq=f[l]*1ll*l%MOD;
        for(int i=l;i<=n;i+=l) (g[i]+=qwq)>=MOD&&(g[i]-=MOD);
        return;
    }
    DandQ_FFT(l,mid);
    if(l==1){
        static int p1[MAXN],p2[MAXN],res[MAXN];
        F(i,l,mid) p1[i-l]=f[i],p2[i-l]=g[i];
        int len=mul((mid-l)<<1,p1,p2,res);
        F(i,mid+1,r) (f[i]+=res[i-l-1])>=MOD&&(f[i]-=MOD);
        F(i,0,len) p1[i]=p2[i]=res[i]=0;
    }else{
        static int p1[MAXN],p2[MAXN],res[MAXN];
        F(i,l,mid) p1[i-l]=f[i];
        F(i,1,r-l) p2[i-1]=g[i];
        int len=mul(r+mid-l-l+1,p1,p2,res);
        F(i,mid+1,r) (f[i]+=res[i-l-1])>=MOD&&(f[i]-=MOD);
        F(i,0,len) p1[i]=p2[i]=res[i]=0;
        F(i,l,mid) p1[i-l]=g[i];
        F(i,1,r-l) p2[i-1]=f[i];
        len=mul(r-l+mid-l+1,p1,p2,res);
        F(i,mid+1,r) (f[i]+=res[i-l-1])>=MOD&&(f[i]-=MOD);
        F(i,0,len) p1[i]=p2[i]=res[i]=0;
    }
    DandQ_FFT(mid+1,r);
    return;
}
signed main(){
    ios::sync_with_stdio(0);
    cin.tie(0),cout.tie(0);
    init();
    cin>>n;
    DandQ_FFT(1,n);
    ll ans=f[n];
    F(i,(n>>1)+1,n-1) ans=(ans-f[i]*1ll*f[n-i])%MOD;
    if(!(n&1)) ans=(ans-f[n>>1]*(f[n>>1]-1ll)/2)%MOD;
    ans<0&&(ans+=MOD);
    cout<<ans;
    return 0;
}

Day 87

很神秘的题目。

2025/01/26 [ARC055D] 隠された等差数列

题意

给定一个长为 n 的数列 d_0,\cdots,d_{n-1}，满足 d_0\neq 0,d_i\in[0,9]\cap\mathbb N。求最小的 A 满足存在首项为 A 公差为 B 的等差数列 a_0,a_1,\cdots 和正整数 x 使得 \forall i\in[0,n-1],d_i=\lfloor\dfrac{a_i}{10^{x-1}}\rfloor，或报告无解。

思路

先从无解情况开始思考。发现 d 的增量 d_{i+1}-d_i 可以看做是对位加再加上进位的贡献。进位最多 +1，所以 d_{i+1}-d_i 在模 10 意义下只有两种取值，且相差 1。如果不满足这个条件，必然无解。

然后只有一种取值也是好处理的，此时这种取值必然是 0，取 A=d_i,B=0,x=1 即可。

接下来我们假设 d_{i+1}-d_i 有差为 1 的两种取值。

首先你会发现 A,B<10^x，因为超过的部分不会对 d 造成影响。

我们设 l_i=\lfloor\dfrac{a_i}{10^x}\rfloor。根据上面的观察 l_0 最低位是 d_0 其他位是 0。

然后你发现 d_i<d_{i+1} 代表第 x 位不进位，d_i>d_{i+1} 代表第 x 位进位，相等的情况可以用另一种 d_{i+1}-d_i 的取值是 1 还是 9 来判断 是否进位，1 代表不进位 9 代表进位。根据这些第 x 位是否进位的信息，我们可以求出所有的 l_i。

于是我们就有不等式组 l_i\times10^x\le A+Bi<(l_i+1)\times10^x。设 A'=\dfrac{A}{10^x},B'=\dfrac{B}{10^x}，同时除以 10^x 得到 l_i\le A'+B'i<l_i+1。

从小到大枚举 x，因为 d_0\neq 0，故如果此时存在 A 必有 10^x\le A<10^{x+1}。枚举到第一个存在 A 的 x 就可以直接输出答案了。

问题转化为给定平面上若干线段 x=i,y\in[l_i,l_i+1)，求是否存在一条斜率为 B' 截距为 A' 的直线穿过所有线段，若存在求 A' 最小值。

发现若 B' 存在，必然存在某条直线穿过某个 (i,l_i) 使得 A' 也存在，且 B' 越大 A' 越小。所以我们只需要求 B' 的取值范围即可，而且这个范围和 x 无关。

记点集 P=\{(i,l_i)|i=0,1,\cdots n\}，Q=\{(i,l_i+1)|i=0,1,\cdots n\}，限制改写成 P 都在直线上或在直线下方，Q 都在直线上方。

进一步的，你发现只有 P 的上凸壳和 Q 的下凸壳中的点会做出贡献。设凸包大小为 k，

求出凸包后，斜率具有单调性，双指针的时间复杂度 O(k)。直接枚举点对是 O(k^2)。

官方题解说 k=O(n^{\frac 2 3})，所以上面两种方法都是可行的。

答案的上界是 O(n^2)，所以枚举 x 的复杂度是 O(\log n)，不是瓶颈。

总的时间复杂度是 O(n) 或 O(n^{\frac 4 3})，取决于找 B' 的范围如何实现。

好吧，最后懒了，考虑到题目年代久远，n\le 10^4，O(n^2) 都能跑。

代码

#include<bits/stdc++.h>
#define F(i,a,b) for(int i(a),i##i##end(b);i<=i##i##end;++i)
#define R(i,a,b) for(int i(a),i##i##end(b);i>=i##i##end;--i)
#define ll long long
#define File(a) freopen(#a".in","r",stdin);freopen(#a".out","w",stdout)
using namespace std;
const int MAXN=1e4+2;
const double EPS=1e-9;
char s[MAXN];
int n,d[MAXN],l[MAXN];
set<int>delta;
signed main(){
    ios::sync_with_stdio(0);
    cin.tie(0),cout.tie(0);
    cin>>s;
    n=strlen(s);
    F(i,0,n-1) d[i]=s[i]-'0';
    F(i,0,n-2) delta.insert((d[i+1]-d[i]+10)%10);
    if(delta.size()>2) return cout<<"-1",0;
    if(delta.size()==1||n==1) return cout<<d[0],0;
    int x=*delta.begin(),y=*delta.rbegin();
    if((x-y+10)%10!=1&&(y-x+10)%10!=1) return cout<<"-1",0;
    l[0]=d[0];
    if((y-x+10)%10==9) x=9,y=10;
    F(i,1,n-1){
        if((l[i-1]+x)%10==d[i]) l[i]=l[i-1]+x;
        else l[i]=l[i-1]+y;
    }
    double bl=0,br=1e9;
    F(i,0,n-1) F(j,i+1,n-1){
        bl=max(bl,(l[j]-l[i]-1.0+EPS)/(j-i));
        br=min(br,(l[j]+1.0-l[i]-EPS)/(j-i));
    }
    if(bl>br) return cout<<"-1",0;
    for(int i=1;;i*=10) if(ceil(bl*i)<=floor(br*i)){
        int b=floor(br*i),res=0;
        F(j,0,n-1) res=max(res,l[j]*i-j*b);
        return cout<<res,0;
    }
    return 0;
}#include<bits/stdc++.h>
#define F(i,a,b) for(int i(a),i##i##end(b);i<=i##i##end;++i)
#define R(i,a,b) for(int i(a),i##i##end(b);i>=i##i##end;--i)
#define ll long long
#define File(a) freopen(#a".in","r",stdin);freopen(#a".out","w",stdout)
using namespace std;
const int MAXN=1e4+2;
const double EPS=1e-9;
char s[MAXN];
int n,d[MAXN],l[MAXN];
set<int>delta;
signed main(){
    ios::sync_with_stdio(0);
    cin.tie(0),cout.tie(0);
    cin>>s;
    n=strlen(s);
    F(i,0,n-1) d[i]=s[i]-'0';
    F(i,0,n-2) delta.insert((d[i+1]-d[i]+10)%10);
    if(delta.size()>2) return cout<<"-1",0;
    if(delta.size()==1||n==1) return cout<<d[0],0;
    int x=*delta.begin(),y=*delta.rbegin();
    if((x-y+10)%10!=1&&(y-x+10)%10!=1) return cout<<"-1",0;
    l[0]=d[0];
    if((y-x+10)%10==9) x=9,y=10;
    F(i,1,n-1){
        if((l[i-1]+x)%10==d[i]) l[i]=l[i-1]+x;
        else l[i]=l[i-1]+y;
    }
    double bl=0,br=1e9;
    F(i,0,n-1) F(j,i+1,n-1){
        bl=max(bl,(l[j]-l[i]-1.0+EPS)/(j-i));
        br=min(br,(l[j]+1.0-l[i]-EPS)/(j-i));
    }
    if(bl>br) return cout<<"-1",0;
    for(int i=1;;i*=10) if(ceil(bl*i)<=floor(br*i)){
        int b=floor(br*i),res=0;
        F(j,0,n-1) res=max(res,l[j]*i-j*b);
        return cout<<res,0;
    }
    return 0;
}

Day 88

考查选手补题能力。

2025/01/27 [ARC166E] Fizz Buzz Difference

题意

所以 f(l,r)=q-p-2-\lfloor\dfrac{(q-p)a-\gcd(a,b)-1}{b}\rfloor\le n。这个式子关于 q-p 单调，可以二分出 q-p 的最大值。下面我们认为 q-p 已经确定且最大。

接下来是最小化 l，即最小化 p。

回到题目本身的限制，我们要有 [l,r] 中出现至少 q-p-1-n 个 b 的倍数。设 k=q-p-1-n，则 r-l-(r\bmod b)\ge b(k-1)，即 (pa+(q-p)a-1)\bmod b\le (q-p)a-2-b(k-1)。

如果 (q-p)a-2-b(k-1)\ge b，无论如何这个不等式都成立，取 l=1 即可。

否则，化简不等式得到 (pa)\bmod b\ge(1-(q-p)a)\bmod b。我们要求满足这个式子的最小的 p。

这是经典问题。设 g(u,v,t,m) 表示满足 u\le(tx)\bmod m\le v 的最小的非负整数 x，其中 u,v,t\in[0,m)。注意此处 t 和上面没有关系。那么此时答案为 g((1-(q-p)a)\bmod b,b-1,a,b)。接下来研究 g 如何求。

如果 [u,v] 中有 t 的倍数，则 g(u,v,t,m)=\lceil\dfrac u t\rceil；如果 u=0，则 g(u,v,t,m)=0，这两种情况好处理的。

否则我们设 s=\lfloor\dfrac{tx}m\rfloor，不等式转化为 u\le tx-sm\le v。

由于 [u,v] 中没有 t 的倍数，不等式两边同时对 t 取模再乘 -1，得到 (-v)\bmod t\le(sm)\bmod t\le(-u)\bmod t。这是一个子问题，取 s=g((-v)\bmod t,(-u)\bmod t,m\bmod t,t) 即可回带求出 x 最小值。

时间复杂度同辗转相除为 O(T\log n)。

代码

#include<bits/stdc++.h>
#define F(i,a,b) for(int i(a),i##i##end(b);i<=i##i##end;++i)
#define R(i,a,b) for(int i(a),i##i##end(b);i>=i##i##end;--i)
#define ll long long
#define File(a) freopen(#a".in","r",stdin);freopen(#a".out","w",stdout)
using namespace std;
ll n,a,b; 
ll g(ll u,ll v,ll t,ll m){
    if(u==0) return 0;
    if(((u-1)/t+1)*t<=v) return (u-1)/t+1;
    ll s=g((t-v%t)%t,(t-u%t)%t,m%t,t);
    return (s*m+u-1)/t+1;
}
signed main(){
    ios::sync_with_stdio(0);
    cin.tie(0),cout.tie(0);
    int T;
    for(cin>>T;T;--T){
        cin>>n>>a>>b;
        ll ql=0,qr=1e12,d=__gcd(a,b);
        while(ql<qr){
            ll mid((ql+qr)>>1);
            if(mid-2-(mid*a-d-1)/b<=n) ql=mid+1;
            else qr=mid;
        }
        ll len=ql-1,p,q;
        if(len*a-2-b*(len-2-n)>=b) p=0;
        else p=g(((1-len*a)%b+b)%b,b-1,a,b);
        q=p+len;
        cout<<p*a+1<<" "<<q*a-1<<"\n";
    }
    return 0;
}

Day 89

春晚太无聊了，来看看很久远的题目，完成五年前的约定。

2025/01/29 [MtOI2019] 幽灵乐团 / 莫比乌斯反演基础练习题

题意

写出所求式 \prod\limits_{i=1}^a\prod\limits_{j=1}^b\prod\limits_{k=1}^c i。

交换求积顺序 \prod\limits_{j=1}^b\prod\limits_{k=1}^c a!。

即有：

预处理阶乘可以做到单次 O(\log n)。

t=0，求 h

写出所求式 \prod\limits_{i=1}^a\prod\limits_{j=1}^b\prod\limits_{k=1}^c\gcd(i,j)。

不妨假设 a\le b。

整理，枚举 \gcd 得到\left(\prod\limits_{k=1}^a\prod\limits_{i=1}^a\prod\limits_{j=1}^b k^{[gcd(i,j)=k]}\right)^c。

化简一下有 \left(\prod\limits_{k=1}^a k^{\sum\limits_{i=1}^a\sum\limits_{j=1}^b[gcd(i,j)=k]}\right)^c。

指数上是这道题，化简过程不再重复。结果是 \left(\prod\limits_{k=1}^a k^{\sum\limits_{d=1}^{\lfloor\frac{a}{k}\rfloor}\mu(d)\lfloor\frac{a}{kd}\rfloor\lfloor\frac{b}{kd}\rfloor}\right)^c。

枚举 k 变为枚举 kd=x，得到最终结果：

括号里的东西可以 O(n\log n) 或狄利克雷前缀和 O(n\log\log n) 预处理，然后就可以整除分块加快速幂 O(\sqrt n\log n) 查询。

t=1，求 g

写出所求式 \prod\limits_{i=1}^a\prod\limits_{j=1}^b\prod\limits_{k=1}^c i^{ijk}。

即 \prod\limits_{i=1}^a i^{i\sum\limits_{j=1}^b\sum\limits_{k=1}^c jk}。

等差数列求和随便化简一下得到 \prod\limits_{i=1}^ai^{i\frac{b(b+1)}{2}\frac{c(c+1)}{2}}

即有：

写出所求式 \prod\limits_{i=1}^a\prod\limits_{j=1}^b\prod\limits_{k=1}^c \gcd(i,j)^{ijk}。

即 \left(\prod\limits_{i=1}^a\prod\limits_{j=1}^b \gcd(i,j)^{ij}\right)^{\frac{c(c+1)}{2}}。

不妨设 a\le b。

然后枚举 \gcd 有 \left(\prod\limits_{k=1}^a k^{\sum\limits_{i=1}^a\sum\limits_{j=1}^b[gcd(i,j)=k]ij}\right)^{\frac{c(c+1)}{2}}。

写出所求式 \prod\limits_{i=1}^a\prod\limits_{j=1}^b\prod\limits_{k=1}^c i^{\gcd(i,j,k)}。

化简有 \prod\limits_{i=1}^a i^{\sum\limits_{j=1}^b\sum\limits_{k=1}^c\gcd(i,j,k)}。

进行欧拉反演，有 \prod\limits_{i=1}^a i^{\sum\limits_{j=1}^b\sum\limits_{k=1}^c\sum\limits_{d|i,d|j,d|k}\varphi(d)}。

枚举因数变为枚举倍数 \prod\limits_{i=1}^a i^{\sum\limits_{d|i}\varphi(d)\lfloor\frac b {d}\rfloor\lfloor\frac c {d}\rfloor}。

再来一次变成 \prod\limits_{d=1}^a\left(\prod\limits_{i=1}^{\lfloor\frac a {d}\rfloor} id\right)^{\varphi(d)\lfloor\frac b {d}\rfloor\lfloor\frac c {d}\rfloor}。

展开得到 \prod\limits_{d=1}^a\left(d^{\lfloor\frac a {d}\rfloor}\lfloor\dfrac a {d}\rfloor!\right)^{\varphi(d)\lfloor\frac b {d}\rfloor\lfloor\frac c {d}\rfloor}。

即有：

可以 O(n\log n) 预处理，单次整除分块 O(\sqrt n\log n)。

t=2，求 h

大的要来了。

写出所求式 \prod\limits_{i=1}^a\prod\limits_{j=1}^b\prod\limits_{k=1}^c \gcd(i,j)^{\gcd(i,j,k)}。

不妨设 a\le b。

枚举 \gcd(i,j)=x，有 \prod\limits_{i=1}^a\prod\limits_{j=1}^b\prod\limits_{k=1}^c \prod\limits_{x=1}^a x^{[\gcd(i,j)=x]\gcd(x,k)}。

化简一下，有 \prod\limits_{x=1}^a x^{\sum\limits_{k=1}^c\gcd(x,k)\sum\limits_{i=1}^a\sum\limits_{j=1}^b [\gcd(i,j)=x]}。

指数的后两个求和是经典问题，上面已经提到过，代入得 \prod\limits_{x=1}^a x^{\sum\limits_{d=1}^{\lfloor\frac{a}{x}\rfloor}\mu(d)\lfloor\frac{a}{xd}\rfloor\lfloor\frac{b}{xd}\rfloor\sum\limits_{k=1}^c\gcd(x,k)}。

对后一个求和号进行欧拉反演，有 \prod\limits_{x=1}^a x^{\left(\sum\limits_{d=1}^{\lfloor\frac{a}{x}\rfloor}\mu(d)\lfloor\frac{a}{xd}\rfloor\lfloor\frac{b}{xd}\rfloor\right)\sum\limits_{k=1}^c\sum\limits_{d|x,d|k}\varphi(d)}。

枚举因数变枚举倍数有 \prod\limits_{x=1}^a x^{\left(\sum\limits_{d=1}^{\lfloor\frac{a}{x}\rfloor}\mu(d)\lfloor\frac{a}{xd}\rfloor\lfloor\frac{b}{xd}\rfloor\right)\sum\limits_{d|x}\varphi(d)\lfloor\frac c d\rfloor}。

化简得到 \prod\limits_{x=1}^a \left(x^{\sum\limits_{d|x}\varphi(d)\lfloor\frac c d\rfloor}\right)^{\sum\limits_{d=1}^{\lfloor\frac{a}{x}\rfloor}\mu(i)\lfloor\frac{a}{xd}\rfloor\lfloor\frac{b}{xd}\rfloor}。

可惜的是，括号里的式子没法预处理，我们必须继续化简。枚举因数变枚举倍数得 \prod\limits_{d=1}^a\prod\limits_{x=1}^{\lfloor\frac a d\rfloor} (xd)^{\varphi(d)\lfloor\frac c d\rfloor\sum\limits_{i=1}^{\lfloor\frac{a}{xd}\rfloor}\mu(i)\lfloor\frac{a}{xdi}\rfloor\lfloor\frac{b}{xdi}\rfloor}。

分成 x,d 分别计算得 \left(\prod\limits_{d=1}^a\prod\limits_{x=1}^{\lfloor\frac a d\rfloor} x^{\varphi(d)\lfloor\frac c d\rfloor\sum\limits_{i=1}^{\lfloor\frac{a}{xd}\rfloor}\mu(i)\lfloor\frac{a}{xdi}\rfloor\lfloor\frac{b}{xdi}\rfloor}\right)\left(\prod\limits_{d=1}^a\prod\limits_{x=1}^{\lfloor\frac a d\rfloor} d^{\varphi(d)\lfloor\frac c d\rfloor\sum\limits_{i=1}^{\lfloor\frac{a}{xd}\rfloor}\mu(i)\lfloor\frac{a}{xdi}\rfloor\lfloor\frac{b}{xdi}\rfloor}\right)。

先处理右边的括号。x 丢到指数上有 \left(\prod\limits_{d=1}^a\prod\limits_{x=1}^{\lfloor\frac a d\rfloor} x^{\varphi(d)\lfloor\frac c d\rfloor\sum\limits_{i=1}^{\lfloor\frac{a}{xd}\rfloor}\mu(i)\lfloor\frac{a}{xdi}\rfloor\lfloor\frac{b}{xdi}\rfloor}\right)\left(\prod\limits_{d=1}^a d^{\varphi(d)\lfloor\frac c d\rfloor\sum\limits_{x=1}^{\lfloor\frac a d\rfloor}\sum\limits_{i=1}^{\lfloor\frac{a}{xd}\rfloor}\mu(i)\lfloor\frac{a}{xdi}\rfloor\lfloor\frac{b}{xdi}\rfloor}\right)。

枚举 xi=k，得 \left(\prod\limits_{d=1}^a\prod\limits_{x=1}^{\lfloor\frac a d\rfloor} x^{\varphi(d)\lfloor\frac c d\rfloor\sum\limits_{i=1}^{\lfloor\frac{a}{xd}\rfloor}\mu(i)\lfloor\frac{a}{xdi}\rfloor\lfloor\frac{b}{xdi}\rfloor}\right)\left(\prod\limits_{d=1}^a d^{\varphi(d)\lfloor\frac c d\rfloor\sum\limits_{k=1}^{\lfloor\frac a d\rfloor}\lfloor\frac{a}{dk}\rfloor\lfloor\frac{b}{dk}\rfloor\sum\limits_{i|k}\mu(i)}\right)。

代入 \sum\limits_{i|k}\mu(i)=[k=1] 得 \left(\prod\limits_{d=1}^a\prod\limits_{x=1}^{\lfloor\frac a d\rfloor} x^{\varphi(d)\lfloor\frac c d\rfloor\sum\limits_{i=1}^{\lfloor\frac{a}{xd}\rfloor}\mu(i)\lfloor\frac{a}{xdi}\rfloor\lfloor\frac{b}{xdi}\rfloor}\right)\left(\prod\limits_{d=1}^a d^{\varphi(d)\lfloor\frac{a}{d}\rfloor\lfloor\frac{b}{d}\rfloor\lfloor\frac c d\rfloor}\right)。

再处理左边的括号。枚举 xi=k，得到最终式子：

左边括号可以做到 O(n\log n) 预处理最内层，整除分块套整除分块单次 O(n^{\frac 3 4}\log n)。右边括号直接整除分块是 O(\sqrt n\log n) 的。

然后你发现 h(a,b,c) 的右边那个括号和 g(a,b,c) 的前半部分消掉了，不需要处理。内层的整除分块就是 t=0 的 h。

总时间复杂度是 O(n\log n+Tn^{\frac 3 4}\log n)。应该可以用杜教筛预处理的技巧做到把指数降到 \dfrac 2 3。

代码

#include<bits/stdc++.h>
#define F(i,a,b) for(int i(a),i##i##end(b);i<=i##i##end;++i)
#define R(i,a,b) for(int i(a),i##i##end(b);i>=i##i##end;--i)
#define ll long long
#define File(a) freopen(#a".in","r",stdin);freopen(#a".out","w",stdout)
using namespace std;
const int MAXN=1e5+1;
int MOD;
inline ll qpow(ll base,ll expo){
    ll res(1);
    while(expo){
        (expo&1)&&(res=res*base%MOD);
        expo>>=1,base=base*base%MOD;
    }
    return res;
}
bool vis[MAXN];
int pr[MAXN],cnt,mu[MAXN],phi[MAXN];
ll fact[MAXN],sum[MAXN],iimul[MAXN],dmu[MAXN],m1[MAXN],m2[MAXN],im1[MAXN],im2[MAXN];
inline void init(){
    mu[1]=phi[1]=1;
    F(i,2,MAXN-1){
        if(!vis[i]) pr[++cnt]=i,mu[i]=-1,phi[i]=i-1;
        F(j,1,cnt){
            ll k=i*pr[j];
            if(k>=MAXN) break;
            vis[k]=1;
            if(i%pr[j]) mu[k]=-mu[i],phi[k]=phi[i]*(pr[j]-1);
            else{
                phi[k]=phi[i]*pr[j];
                break;
            }
        }
    }
    fact[0]=iimul[0]=1;
    F(i,1,MAXN-1) dmu[i]=1,sum[i]=sum[i-1]+phi[i],fact[i]=fact[i-1]*i%MOD,iimul[i]=iimul[i-1]*qpow(i,i)%MOD;
    F(i,1,MAXN-1) if(mu[i]) for(int j=i,k=1;j<MAXN;j+=i,++k) dmu[j]=dmu[j]*(mu[i]==1?k:qpow(k,MOD-2))%MOD;
    m1[0]=m2[0]=1;
    F(i,1,MAXN-1) m1[i]=m1[i-1]*dmu[i]%MOD,m2[i]=m2[i-1]*qpow(dmu[i],i*1ll*i%(MOD-1))%MOD;
    F(i,0,MAXN-1) im1[i]=qpow(m1[i],MOD-2),im2[i]=qpow(m2[i],MOD-2);
    return;
}
#define C2(x) ((x+1ll)*x/2%(MOD-1))
inline ll g0(int a,int b,int c){
    return qpow(fact[a],b*1ll*c%(MOD-1));
}
inline ll h0(int a,int b,int c){
    ll res=1;
    if(a>b) swap(a,b);
    for(int l=1,r;l<=a;l=r+1){
        ll ta=a/l,tb=b/l;
        r=min(a/ta,b/tb);
        res=res*qpow(m1[r]*im1[l-1]%MOD,ta*tb%(MOD-1))%MOD; 
    }
    return qpow(res,c);
}
inline ll g1(int a,int b,int c){
    return qpow(iimul[a],C2(b)*C2(c)%(MOD-1));
}
inline ll h1(int a,int b,int c){
    ll res=1;
    if(a>b) swap(a,b);
    for(int l=1,r;l<=a;l=r+1){
        ll ta=a/l,tb=b/l;
        r=min(a/ta,b/tb);
        res=res*qpow(m2[r]*im2[l-1]%MOD,C2(ta)*C2(tb)%(MOD-1))%MOD; 
    }
    return qpow(res,C2(c)%(MOD-1));
}
inline ll g2(int a,int b,int c){
    ll res=1;
    int lim=min(a,min(b,c));
    for(int l=1,r;l<=lim;l=r+1){
        ll ta=a/l,tb=b/l,tc=c/l;
        r=min(a/ta,min(b/tb,c/tc));
        res=res*qpow(fact[ta],(sum[r]-sum[l-1])%(MOD-1)*tb%(MOD-1)*tc%(MOD-1))%MOD; 
    }
    return res;
}
inline ll h2(int a,int b,int c){
    ll res=1;
    int lim=min(a,min(b,c));
    for(int l=1,r;l<=lim;l=r+1){
        ll ta=a/l,tb=b/l,tc=c/l;
        r=min(a/ta,min(b/tb,c/tc));
        res=res*h0(ta,tb,(sum[r]-sum[l-1])%(MOD-1)*tc%(MOD-1))%MOD; 
    }
    return res;
}
signed main(){
    ios::sync_with_stdio(0);
    cin.tie(0),cout.tie(0);
    int T;
    cin>>T>>MOD;
    for(init();T;--T){
        int A,B,C;
        cin>>A>>B>>C;
        cout<<g0(A,B,C)*g0(B,A,C)%MOD*qpow(h0(A,B,C)*h0(A,C,B)%MOD,MOD-2)%MOD<<" ";
        cout<<g1(A,B,C)*g1(B,A,C)%MOD*qpow(h1(A,B,C)*h1(A,C,B)%MOD,MOD-2)%MOD<<" ";
        cout<<g2(A,B,C)*g2(B,A,C)%MOD*qpow(h2(A,B,C)*h2(A,C,B)%MOD,MOD-2)%MOD<<"\n";
    }
    return 0;
}

Day 90

一坨题目。

2025/01/30~2025/01/31 [ICPC 2018 Qingdao R] Sequence and Sequence

题意

定义序列 P=(1,1,2,2,2,3,3,3,3,\cdots)，即数 i 出现了 i+1 次且所有数不降的数列，下标从 1 开始。定义 Q_1=1，i\ge 2 时 Q_i=Q_{i-1}+Q_{P_i}。T 次给定 n 询问 Q_n。

思路

不妨令 P_0=Q_0=0。

把递推式展开得到 Q_n=\sum\limits\limits_{i=0}^n Q_{P_i}，即 Q_n 可以视为数列 Q_{P_n} 的前缀和。

既然是求前缀和，我们可以使用有限微积分。

先写成有限微积分的形式 Q_n=\sum_{0}^{n+1} Q_{P_x}\Delta x。

类比微积分的分部积分，对于有限微积分，我们有：

【分部求和法】 \sum u\Delta v=uv-\sum\mathrm E v\Delta u，其中 \mathrm E 是位移算子。

对上面的定和式使用分部求和，得到 Q_n=(xQ_{P_x})|^{n+1}_0-\sum_0^{n+1}(x+1)\Delta Q_{P_x}。

展开得到 Q_n=(n+1)Q_{P_{n+1}}-\sum\limits_{i=0}^{n}(i+1)(Q_{P_{i+1}}-Q_{P_i})。

把 i=n 时单独拎出来，再令 i\leftarrow i+1，化简一下，得到：

你发现大部分情况 Q_{P_i}-Q_{P_{i-1}}=0，\neq 0 的位置可以表示成 \dfrac{i(i+1)}2,i\le P_n。于是上面的式子可以写成 Q_n=nQ_{P_n}-\sum\limits_{i=1}^{P_n}(\dfrac{i(i+1)}2-1)(q_j-q_{j-1})。

简单化简得到 Q_n=nQ_{P_n}-\sum\limits_{i=1}^{P_n}(\dfrac{i(i+1)}2-1)q_{p_j}。右边又是一个类似的求和式。

于是设 f_0(n)=1,F_k(n)=\sum\limits_{i=1}^n f_k(i),f_{k+1}(n)=F_k(\dfrac{n(n+1)}{2}-1)，再设 g(n,k)=F_k(n)q_{p_n}-\sum\limits_{i=1}^{p_n}f_{k+1}(i)q_{p_i}。

则有：