蜜月日记 Part11

littlez_meow · 2025-04-13 19:49:59 · 休闲·娱乐

前言

暂且跳过 Day 100，先往后写吧。

Day 101

复健博弈论。

2025/04/14 [SNOI2020] 取石子

题意

有 t 个石子，先手第一次可以拿不超过 k 个，往后每一次可以拿不超过上一个人拿的石子个数的 2 倍，拿走最后一个石子的输。求有多少 t\in[1,n] 满足先手必胜。

思路

比较板的一道博弈论。

首先 n\leftarrow n-1，变成没有石子拿的输。

我们就有如下定理：

【齐肯多夫定理】 任意正整数可以表示为若干不相邻斐波那契数之和。这种表示方法称为齐肯多夫表示法。

然后在没有 k 的限制下是经典的：

【Fibonacci Nim】 有 n 个石子，先手第一次不能一次取完，往后每一次可以拿不超过上一个人拿的石子个数的 2 倍，没有石子的拿的输。则先手必胜当且仅当 n 不是斐波那契数，先手第一次拿的石子个数的最小值为 n 的齐肯多夫表示法下的最小的斐波那契数。

加上 k 的限制，就是最小的斐波那契数要 \le k。同时还去掉了不能一次取完的限制，所以就没有 n 不是斐波那契数的限制。

那么问题变成求 1\sim n 中齐肯多夫表示法下最小的斐波那契数 \le k 的数字个数，容斥变成求 >k 的数字个数，直接数位 dp 就可以做了。具体地，记录是否顶满 n 和上一位是多少就行。

时间复杂度 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 File(a) freopen(#a".in","r",stdin);freopen(#a".out","w",stdout)
#define ll long long
using namespace std;
const int MAXN=90;
ll fib[MAXN],n,k,dp[MAXN][2][2];
int lim,T;
bool nn[MAXN];
ll dfs(int now,bool eq,bool lst){
    if(now<=lim) return 1;
    if(dp[now][eq][lst]!=-1) return dp[now][eq][lst];
    ll ans=0;
    if(nn[now]){
        if(!lst) ans+=dfs(now-1,eq,1);
        ans+=dfs(now-1,0,0);
    }else{
        if(!eq&&!lst) ans+=dfs(now-1,0,1);
        ans+=dfs(now-1,eq,0);
    }
    return dp[now][eq][lst]=ans;
}
signed main(){
    ios::sync_with_stdio(0);
    cin.tie(0),cout.tie(0);
    fib[1]=1;
    F(i,2,89) fib[i]=fib[i-1]+fib[i-2];
    for(cin>>T;T;--T){
        cin>>k>>n;
        --n;
        memset(nn,0,sizeof nn),memset(dp,-1,sizeof dp);
        ll qwq=n;
        R(i,89,1) if(qwq>=fib[i]) nn[i]=1,qwq-=fib[i];
        R(i,89,0) if(fib[i]<=k){lim=i;break;}
        cout<<n-dfs(89,1,0)+1<<"\n";
    }
    return 0;
} 

Day 102

复健生成函数。

2025/04/17 数树（2021 CoE-II E）

题意

定义一棵无标号有根树是 k_1-k_2 叉的，当且仅当每个节点的儿子个数都是 0,k_1,k_2 之一，且儿子之间有顺序。若一棵 k_1-k_2 叉树有 x 个 k_1 叉节点、y 个 k_2 叉节点，则该树的权值为 ak_1+bk_2。求在所有不同构的 n 个点的 k_1-k_2 叉树中随机选择一棵的权值的期望对 998244353 取模的值。

思路

就喜欢这种儿子之间有顺序的树，可以直接上生成函数。

由于我们还要知道 k_1 叉节点个数和 k_2 叉节点个数，我们设 [x^p y^q]F(x) 表示 p 个点 q 个 k_1 叉节点的树的个数。下面会把 F(x) 简记为 F，把 y 视为常数。

直接就有 F=x(1+yF^{k_1}+F^{k_2})。

移项变成 x=\dfrac{F}{1+yF^{k_1}+F^{k_2}}。设 G(x)=\dfrac{x}{1+yx^{k_1}+x^{k_2}}，则 F(G(x))=x。

上拉反：

即 [x^n]F=\dfrac{1}{n}[x^{n-1}](1+yx^{k_1}+x^{k_2})^n。

我们需要对于所有 k\in[0,n] 求出 [x^{n-1}y^k](1+yx^{k_1}+x^{k_2})。然后期望的分子分母就可以直接算了。

枚举 k，设括号里分别有 r,s,t 次乘的是 1,yx^{k_1},x^{k_2}，则有 \left\{\begin{matrix}r+s+t=n \\sk_1+tk_2=n-1 \\s=k\end{matrix}\right.。这是三元一次方程组，可以解出唯一的 r,s,t，然后贡献就是多重组合数 \dbinom{n}{r,s,t}。

预处理阶乘和阶乘逆元，计算时间复杂度 O(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 File(a) freopen(#a".in","r",stdin);freopen(#a".out","w",stdout)
#define ll long long
using namespace std;
const int MAXN=1e7+1,MOD=998244353;
int n,a,b,k1,k2,fact[MAXN],inv[MAXN];
inline ll qpow(ll b,int e){ll res=1;while(e)(e&1)&&(res=res*b%MOD),b=b*b%MOD,e>>=1;return res;}
signed main(){
    ios::sync_with_stdio(0);
    cin.tie(0),cout.tie(0);
    cin>>k1>>k2>>n>>a>>b;
    ll invn=qpow(n,MOD-2);
    fact[0]=1;
    F(i,1,n) fact[i]=fact[i-1]*1ll*i%MOD;
    inv[n]=qpow(fact[n],MOD-2);
    R(i,n,1) inv[i-1]=inv[i]*1ll*i%MOD;
    int ch=0,mo=0;
    F(k,0,n){
        ll r,s=k,t=n-1-s*k1;
        if(t%k2) continue;
        t=t/k2,r=n-s-t;
        if(r<0||s<0||t<0) continue;
        ll qwq=invn*fact[n]%MOD*inv[r]%MOD*inv[s]%MOD*inv[t]%MOD;
        ch=(ch+(s*a+t*b)%MOD*qwq)%MOD;
        (mo+=qwq)>=MOD&&(mo-=MOD);
    }
    cout<<ch*qpow(mo,MOD-2)%MOD;
    return 0;
} 

Day 103

概率不取模，肯定有问题。

2025/04/18 [PA 2025] 考试 / Egzamin

题意

有 n 道题，第 i 道正确率为 p_i，答对、不答、答错分别可以得到 1,0,-1 分。你可以自由选择作答哪些题，求得分 \ge t 的最大概率。

思路

首先 p 从大到小排序，贪心来看每次作答的必然是前缀。

于是设 dp(i,j) 表示当前答前 i 个题获得 j 分的概率。

转移显然为 dp(i,j)=dp(i-1,j-1)p_{i-1}+dp(i-1,j+1)(1-p_{i-1})。可以做到 O(n^2)。

发现有很多 dp(i,j) 都很小，我们不妨只对 dp(i,j)>\epsilon 转移。

为什么这样做时间复杂度是对的？

我们需要如下定理：

【中心极限定理】 对于足够多独立同分布的随机变量，若它们的均值为 \mu、方差为 \sigma^2，它们的和近似于均值为 \mu、方差为 n\sigma^2 的正态分布。

本题中，每道题的得分服从两点分布，方差可以视为常数。根据中心极限定理，总得分应服从正态分布。

如果你不会很多概率论知识，我们可以用人教版数学选择性必修三的内容解释。数学书上提到了 3\sigma 原则——对于均值为 \mu、方差为 \sigma^2 的正态分布，绝大部分样本都落在 [\mu-3\sigma,\mu+3\sigma] 的地方。由于每道题得分的方差视为常数，故总得分服从的正态分布的方差 \sigma^2=O(n)，也就是绝大部分样本落在的区间长度约为 O(\sqrt n)。

如果你学过概率论，你应该听说过切尔诺夫引理。具体内容比较复杂就不展开了，这个引理表明 >\epsilon 的样本落在的区间长度为 O(\sqrt{-n\log\epsilon})。故时间复杂度为 O(n\sqrt{-n\log\epsilon})。如果把精度视为常数，这和上面使用高中课本知识推出来的吻合。

代码

#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 File(a) freopen(#a".in","r",stdin);freopen(#a".out","w",stdout)
#define ll long long
#define fi first
#define se second
using namespace std;
const int MAXN=5e4+1;
const double EPS=1e-12;
int n,t;
double p[MAXN],ans[MAXN]; 
__gnu_pbds::cc_hash_table<int,double>dp[2];
signed main(){
    ios::sync_with_stdio(0);
    cin.tie(0),cout.tie(0);
    cin>>n>>t;
    F(i,1,n) cin>>p[i];
    sort(p+1,p+n+1,greater<double>());
    dp[0][n]=1;
    F(i,0,n-1){
        int now=i&1,nxt=now^1;
        dp[nxt].clear();
        for(auto qwq:dp[now]) if(qwq.se>EPS){
            if(qwq.fi>=n+t) ans[i]+=qwq.se;
            dp[nxt][qwq.fi+1]+=qwq.se*p[i+1];
            dp[nxt][qwq.fi-1]+=qwq.se*(1-p[i+1]);
        }
    }
    for(auto qwq:dp[n&1]) if(qwq.se>EPS) if(qwq.fi>=n+t) ans[n]+=qwq.se;
    double res=0;
    F(i,1,n) res=max(res,ans[i]);
    cout<<fixed<<setprecision(10)<<res;
    return 0;
} 

Day 104

有关这个 3B1B 视频的题目数量还在提升。

2025/04/20 二平方和定理

题意

求 x^2+y^2=n 的所有非负整数解。

思路

我们把问题丢到复平面上，就是求 (x+yi)(x-yi)=n 的所有解。

考虑在高斯整环 \mathbb Z[\sqrt{-1}]（就是实部和虚部都是整数的复数集合）分解质因数。我们接下来先讨论分解质因数的方式，再讨论如何得到答案。

介绍如下定理：

【费马平方和定理】 奇质数 p 能表示为两个平方数之和当且仅当 p\equiv 1\pmod 4，且在不考虑顺序的情况下表示方法唯一。

根据费马平方和定理，对于一个奇质数 p，若 p\equiv 1\pmod 4，那么它可以表示成一对共轭复数之积；否则，p 在高斯整环下仍是质数，即高斯素数。

也就是说，我们只需要把 n 现在整环上质因数分解，再把 n 的模 4 余 1 的质因数分解成一对共轭复数之积即可。

设奇质数 p=z\bar z 满足 p\equiv 1\pmod 4。我们若能找到整数 x 满足 x^2\equiv -1\pmod p，即 p|(x+i)(x-i)，则 z=\gcd(x+i,p)。

根据二次剩余相关知识，x 是存在的。我们随机一个 p 的二次非剩余 y，期望两次随到，那么 x\equiv y^{\frac{p-1}4}\pmod p。

由于高斯整环是辗转相除环，这里的 \gcd 可以用辗转相除法做，a\bmod b 定义为：将 \dfrac a b 的实部和虚部都四舍五入得到 q，则 a\bmod b=a-bq。由于模长每次减半，复杂度仍是 O(\log n)。

那么我们已经可以把 n 在高斯整环上质因数分解了。设 n=z\bar z，我们应该如何分配这些高斯素数到 z,\bar z 中呢？

先处理模 4 余 3 型质数。这些质数本身就是高斯素数，为了让两个因数共轭，必须均分到 z,\bar z 中。如果不能均分，答案就是 0。

然后是模 4 余 1 型质数。设其中一个质数 p=x\bar x，它的指数为 c，那么我们可以分配 t 个 x 和 c-t 个 \bar x 给 z，剩下的给 \bar z。一共 c 种分配方法。对每个这类质数枚举每一种分配方法，时间复杂度的上界是 O(\tau(n))，其中 \tau 是约数个数函数。

最后是 2。你发现虽然 2=(1+i)(1-i)，但无论怎么分配，得到的解都只是旋转了 90^\circ 的区别，我们可以在最后一步全部转到第一象限，所以怎么分都无所谓。

时间复杂度为 O(T(n^{\frac 1 4}+\omega(n)\log n+\tau(n)\log n))，第一项是整环上分解质因数，第二项是在高斯整环上分解（需要 \gcd），第三项是枚举分配方式（需要快速幂）。

代码

#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 File(a) freopen(#a".in","r",stdin);freopen(#a".out","w",stdout)
#define ll long long
#define int ll
#define fi first
#define se second
using namespace std;
__gnu_pbds::gp_hash_table<ll,int>e;
inline ll qpow(__int128 base,ll expo,const ll MOD){
    ll res(1);
    while(expo){
        (expo&1)&&(res=res*base%MOD);
        base=base*base%MOD,expo>>=1;
    }
    return res;
}
namespace Factorize{
    const int MR_Prime[]={2,325,9375,28178,450775,9780504,1795265022};
    inline bool MR(ll n){
        if(n<=2) return n==2;
        if(!(n&1)) return 0;
        ll expo=n-1;
        int r=__builtin_ctzll(expo);
        expo>>=r;
        for(int i:MR_Prime){
            ll v=qpow(i,expo,n);
            if(v<=1||v==n-1) continue;
            F(j,1,r){
                v=__int128(v)*v%n;
                if(v==n-1&&j!=r){
                    v=1;
                    break;
                }
                if(v==1) return 0;
            }
            if(v!=1) return 0;
        }
        return 1;
    }
    __int128 gcd(__int128 x,__int128 y){return y==0?x:gcd(y,x%y);}
    __uint128_t SQUFOF(__uint128_t value) {
        __uint128_t s = sqrtl(1.0l * value) + 1e-6;
        while (s * s > value) s--;
        if (s * s == value) return s;
        __int128 d, po, p, p_prev, q, q_prev, q_, b, r;
        long long l, b_, i;
        int k = 0;
        l = 2 * sqrtl(2.0l * s);
        b_ = 3 * l;
        while (++k) {
            d = k * value;
            po = p_prev = p = sqrtl(1.0l * d);
            q_prev = 1;
            q = d - po * po;
            while (q < 0) {
                po--;
                p_prev--;
                p--;
                q = d - po * po;
            }
            if (!q) return gcd(value, k);
            for (i = 2; i < b_; i++) {
                b = (po + p) / q;
                p = b * q - p;
                q_ = q;
                q = q_prev + b * (p_prev - p);
                r = sqrtl(1.0 * q) + 1e-6;
                if (!(i & 1) && r * r == q) break;
                q_prev = q_;
                p_prev = p;
            }
            if (i >= b_) continue;
            b = (po - p) / r;
            p_prev = p = b * r + p;
            q_prev = r;
            q = (d - p_prev * p_prev) / q_prev;
            i = 0;
            do {
                b = (po + p) / q;
                p_prev = p;
                p = b * q - p;
                q_ = q;
                q = q_prev + b * (p_prev - p);
                q_prev = q_;
                i++;
            } while (p != p_prev && i < b_);
            if (i >= b_) continue;
            r = gcd(value, q_prev);
            if (r != 1) return r;
        }
        return 0;
    }
    inline void factorize(ll n,ll cnt=1){
        if(n==1) return;
        if(MR(n)) return e[n]+=cnt,void();
        ll qaq=SQUFOF(n);
        while(qaq==1||qaq==n) qaq=SQUFOF(n);
        int nc=0;
        while(!(n%qaq)) n/=qaq,++nc;
        factorize(n,cnt),factorize(qaq,nc*cnt);
        return;
    }
} 
struct Cpx{
    __int128 a,b;//a+bi;
    Cpx(const ll&x=0,const ll&y=0):a(x),b(y){}
    Cpx operator+(const Cpx&x)const{return Cpx(a+x.a,b+x.b);}
    Cpx operator-(const Cpx&x)const{return Cpx(a-x.a,b-x.b);}
    Cpx operator*(const Cpx&x)const{return Cpx(a*x.a-b*x.b,a*x.b+b*x.a);}
    Cpx operator/(const Cpx&x)const{
        double qwq=x.a*x.a+x.b*x.b;
        return Cpx(round((a*x.a+b*x.b)/qwq),round((b*x.a-a*x.b)/qwq));
    }
    Cpx operator%(const Cpx&x)const{
        return *this-x*(*this/x);
    }
    Cpx operator^(ll ex)const{
        Cpx res(1),b=*this;
        while(ex) (ex&1)&&(res=res*b,1),b=b*b,ex>>=1;
        return res;
    } 
    bool operator==(const Cpx&x)const{return a==x.a&&b==x.b;}
    __int128 abs()const{return a*a+b*b;}
};
Cpx gcd(const Cpx&x,const Cpx&y){
    if(x.abs()<y.abs()) return gcd(y,x);
    return y==Cpx()?x:gcd(y,x%y);
}
int cnt,ex[100];
ll pr[100];
Cpx coef[100];
mt19937 gen(114514);
vector<pair<ll,ll> >ans;
void dfs(int step,const Cpx&now){
    if(step==cnt+1){
        ll x=now.a,y=now.b;
        if(x>=0&&y>=0) ans.emplace_back(x,y);
        if(x<=0&&y<=0) ans.emplace_back(-x,-y);
        if(y<=0&&x>=0) ans.emplace_back(-y,x);
        if(y>=0&&x<=0) ans.emplace_back(y,-x);
        return;
    }
    if((pr[step]&3)!=1) dfs(step+1,now*coef[step]);
    else{
        Cpx bar(coef[step].a,-coef[step].b),qwq=bar^ex[step];
        F(i,0,ex[step]) dfs(step+1,now*qwq),qwq=qwq/bar*coef[step];
    }
    return;
}
inline void solve(){
    cnt=0;
    for(auto i:e){
        pr[++cnt]=i.fi,ex[cnt]=i.se;
        if(i.fi==2) coef[cnt]=Cpx(1,1)^i.se;
        else if((i.fi&3)==1){
            ll qwq;
            while(1){
                qwq=gen()%(i.fi-1)+1;
                if(qpow(qwq,i.fi>>1,i.fi)==i.fi-1)break;
            }
            qwq=qpow(qwq,i.fi>>2,i.fi);
            coef[cnt]=gcd(Cpx(i.fi),Cpx(qwq,1));
        }else if(i.se&1) return cout<<"0\n",void();
        else coef[cnt]=Cpx(qpow(i.fi,i.se>>1,1e18));
    }
    vector<pair<ll,ll> >().swap(ans);
    dfs(1,Cpx(1,0));
    sort(ans.begin(),ans.end());
    cout<<ans.size()<<"\n";
    for(auto i:ans) cout<<i.fi<<" "<<i.se<<"\n";
    return;
} 
signed main(){
    ios::sync_with_stdio(0);
    cin.tie(0),cout.tie(0);
    ll T,n;
    for(cin>>T;T;--T){
        e.clear();
        cin>>n;
        if(n<=1e10){
            for(int i=2;i*i<=n;++i) if(n%i==0) while(n%i==0) n/=i,++e[i];
            if(n!=1) ++e[n];
        }else Factorize::factorize(n);
        solve();
    }
    return 0;
}

Day 105

感觉自己组合数学有很多知道但是不会推的结论。

2025/04/22 「EZEC-5」「KrOI2021」Chasse Neige 及其加强版

题意

求长为 n 的排列 p 的个数，满足 p_1<p_2,p_{n-1}>p_n，且恰好存在 k 个 i\in[2,n-1] 满足 p_{i-1}<p_i>p_{i+1}。

思路

首先非常显然是排列插入 dp，设 dp(x,y,n,k) 表示 x=[p_1<p_2],y=[p_{n-1}>p_n]、有 k 个峰（不包含首尾）的 1\sim n 的排列个数，然后转移就是往排列里插入 n+1，按照插入的位置分类即可。

你发现 dp(0,0,n,k) 和 dp(1,1,n,k) 很对称，dp(1,0,n,k) 和 dp(0,1,n,k) 很对称。稍加推理可以得到 dp(1,0,n,k)=dp(0,1,n,k),dp(0,0,n,k)=dp(1,1,n,k+1)。化简后这两个的转移方程相当的像。

所以设 f(i,2j)=dp(1,1,i,j),f(i,2j+1)=dp(0,1,i,j)，则有转移方程为 f(i,j)=jf(i-1,j)+(i-j)f(i-1,j-2)+2f(i-1,j-1)。

然后发现数据范围里 n-2k 很小，我们可以先求出对角线 f(i,i-1) 再反过来递推。

而这是经典的 zig-zag 排列计数，我们有 f(i,i-1)=[\dfrac{x^i}{i!}](\tan x+\sec x)=[\dfrac{x^i}{i!}]\dfrac{\sin x+1}{\cos x}，大致思路是分奇偶讨论，列出 dp 方程，发现可以用生成函数求导描述，然后变成微分方程，解得奇偶分别是 \tan x,\sec x。

使用多项式求逆，预处理所有答案，时间复杂度 O(n\log n+n(n-2k)+T)。

可能需要 DIF-DIT 写法的 NTT 来卡常。

代码

#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 File(a) freopen(#a".in","r",stdin);freopen(#a".out","w",stdout)
#define ll long long
using namespace std;
const int MAXN=(1<<23)+5,MOD=998244353,G=3;
inline ll qpow(ll base,ll 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[MAXN];
inline void init(){
    powg[0]=1;
    ll qwq=qpow(G,(MOD-1)>>23);
    F(i,0,21){
        ll qaq=qpow(qwq,1<<(21-i));
        F(j,1<<i,(1<<(i+1))-1) powg[j]=qaq*powg[j-(1<<i)]%MOD;
    }
    return;
}
inline void DIF(int*poly,int len){
    for(int i=len>>1;i>=1;i>>=1) for(int j=0,t=0;j<len;j+=(i<<1),++t) F(k,0,i-1){
        ll x=poly[j|k],y=poly[i|j|k]*1ll*powg[t]%MOD;
        (poly[j|k]=x+y)>=MOD&&(poly[j|k]-=MOD);
        (poly[i|j|k]=x-y)<0&&(poly[i|j|k]+=MOD);
    }
    return;
}
inline void DIT(int*poly,int len){
    for(int i=1;i<len;i<<=1) for(int j=0,t=0;j<len;j+=(i<<1),++t) F(k,0,i-1){
        ll x=poly[j|k],y=poly[i|j|k];
        (poly[j|k]=x+y)>=MOD&&(poly[j|k]-=MOD);
        poly[i|j|k]=(MOD+x-y)*powg[t]%MOD;
    }
    ll invl=qpow(len,MOD-2);
    F(i,0,len-1) poly[i]=poly[i]*invl%MOD;
    reverse(poly+1,poly+len);
    return;
}
inline void inv(int*coef,int*res,int len){
    res[0]=qpow(coef[0],MOD-2);
    for(int l=2;l<=len<<1;l<<=1){
        static int tmp[MAXN];
        int qwq=l<<1;
        memcpy(tmp,coef,sizeof(int)*l);
        memset(tmp+l,0,sizeof(int)*(qwq-l));
        DIF(tmp,qwq),DIF(res,qwq);
        F(i,0,qwq-1) (res[i]=(2-res[i]*1ll*tmp[i])%MOD*res[i]%MOD)<0&&(res[i]+=MOD);
        DIT(res,qwq);
        memset(res+l,0,sizeof(int)*(qwq-l));
    }
    F(i,len+1,len<<1) res[i]=0;
    return;
}
int T,n,k,fact[MAXN],finv[MAXN],mo[MAXN],im[MAXN],ch[MAXN],f[21][MAXN];
signed main(){
    ios::sync_with_stdio(0);
    cin.tie(0),cout.tie(0);
    init();
    cin>>T>>n;
    n+=20;
    fact[0]=1;
    F(i,1,n+1) fact[i]=fact[i-1]*1ll*i%MOD;
    finv[n+1]=qpow(fact[n+1],MOD-2);
    R(i,n+1,1) finv[i-1]=finv[i]*1ll*i%MOD;
    for(int i=0,v=0;i<=n+1;i+=2,v^=1){
        mo[i]=(v?MOD-finv[i]:finv[i]);
        i<=n&&(ch[i+1]=(v?MOD-finv[i+1]:finv[i+1]));
    }
    ch[0]=1;
    inv(mo,im,n+1);
    int len=1<<(__lg(2*n+1)+1);
    DIF(im,len),DIF(ch,len);
    F(i,0,len-1) ch[i]=ch[i]*1ll*im[i]%MOD;
    DIT(ch,len);
    F(i,0,n) f[1][i]=fact[i]*1ll*ch[i]%MOD;
    static ll iv[21]={0,1};
    F(j,2,20) iv[j]=MOD-MOD/j*iv[MOD%j]%MOD;
    F(j,1,19) F(i,1,n) (f[j+1][i-1]=(f[j][i]-(i-j+0ll)*f[j-1][i-1]-2ll*f[j][i-1])%MOD*iv[j]%MOD)<0&&(f[j+1][i-1]+=MOD);
    for(;T;--T){
        cin>>n>>k;
        cout<<f[n-2*k][n]<<"\n";
    }
    return 0;
}

Day 106

是不是好久没放交互了。

2025/04/26 [RMI 2020] 寻线 / Nice Lines

题意

有 n 条斜率不同的直线，第 i 条是 y=a_i x+b_i，满足 |a_i|,|b_i|\in[0,10^4]\cap\mathbb Z。你可以询问不超过 4n+2 次，每次给定一个点 (x,y)，满足 |x|,|y|\in[0,10^{12}]\cap\mathbb R，交互库会返回这个点到所有直线的距离和。你需要找出这 n 条直线。

思路

计算几何值域却开这么小，肯定有说法。下面记 V 为值域。

如果我们能找到一个 x_0>>V（这里 >> 是远大于）并求出 y_0=a_i x_0+b_i，就可以把 y_0 对 x_0 取模得到 b_i，然后就可以算出 a_i。

那么我们就只询问直线 l:x=x_0>>V 上的点。设我们询问 (x_0,t) 返回的结果是 f(t)。对于每条直线而言，到 l 上某点的距离是一个绝对值函数，所以 f(t) 是 n 个绝对值函数加在一起，是一个下凸壳，拐点就是某条直线于 l 相交了。问题变为找到所有拐点。

考虑分治。设当前分治到值域区间 [l,r]，如果 \dfrac{f(l)+f(r)}{2}=f(\dfrac{l+r}{2}) 则区间内没有斜率变化，直接返回，否则从中点递归下去。这样的询问次数是 O(n\log V) 的。

一个非常自然的想法是把分治从值域上搬到拐点序列上。我们希望如果只有区间内只有一个拐点就找到它并且不再递归。首先，我们在足够大和足够小的坐标上找到区间内最靠左和最靠右的组成凸壳的直线，找到它们的交点 (p,q)。如果这个交点在凸壳上（即 q=f(p)），则代表这个区间内只有这一个交点，可以求出来对应的直线然后返回。否则在 p 处找到组成凸壳的直线并递归下去。

这里说的在某处找到组成凸壳的直线可以使用两次询问直接求斜率。当 x_0>>V 时，一个小区间内（具体地，当 x\ge 3V 时，长度为 2V 的区间中至多存在一处斜率变化，证明考虑极端情况就行）是不存在斜率变化的，所以这样做是对的。

分析这样的询问次数，我们初始要 4 次询问出最左最右，递归树的每个节点（共 2n-1 个）都需要 1 次判断是否递归，每个非叶子节点（共 n-1 个）都还需要 2 次找交点处的直线。总的次数是 4+(2n-1)+2(n-1)=4n+1。这是可以接受的。

代码

#include<bits/stdc++.h>
using namespace std;
long double query(long double x, long double y);
void the_lines_are(std::vector<int> a, std::vector<int> b);
namespace MySol{
#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)
#define double long double
#define Line pair<double,double>
#define k first
#define b second
const int x0=1e5,V=2e4;
const double EPS=1e-5;
int n;
vector<int>ansa,ansb;
inline Line find(double t){
    ll qwq=floor(t/x0);
    double x=qwq*x0+V,y=query(x0,x),k=(y-query(x0,(qwq+1)*x0-V))/(2*V-x0);
    return {k,y-k*x};

}
void solve(Line l,Line r){
    double mid=(r.b-l.b)/(l.k-r.k),fmid=query(x0,mid);
    if(fabsl(l.k*mid+l.b-fmid)<EPS&&fabsl(r.k*mid+r.b-fmid)<EPS){
        ll k=round(mid/x0);
        ansa.push_back(k),ansb.push_back(round(mid-x0*k));
        return;
    }
    Line m=find(mid);
    solve(l,m),solve(m,r);
    return;
}
void main(){
    solve(find(-2e9),find(2e9));
    the_lines_are(ansa,ansb);
    return;
}
#undef double
} 
void solve(int subtask_id, int N){MySol::n=N;MySol::main();}

Day 107

之前被卡常了，拿新多项式板子写写。

2025/05/05 [集训队互测 2012] calc 及其加强版

题意

给定 k，对于所有 n\in[1,m]\cap\mathbb Z，求元素两两不相同的值域在 [1,k]\cap\mathbb Z 长度为 n 的所有序列的积的和对 998244353 取模的值。

思路

求出无序的答案再乘 n! 即可。下面研究无序的答案。

显然可以写出答案的生成函数为：

使用指数对数技巧化简为：

接下来处理 \sum\limits_{i=1}^k\ln(1+ix)，最后多项式 \exp 就行。

使用处理对数的经典技巧，先求导再积分，我们有 \left(\ln(1+ix)\right)'=\dfrac{i}{1+ix}，可以得到：

对积分里面的做幂级数展开，交换积分与求和号：

把积分积出来：

补上一个 -1，整理：

交换求和号：

只要能快速求一列自然数幂和就能求出这个多项式。考虑到伯努利数的形式，我们分析一列自然数幂和的 EGF：\sum\limits_{t\ge 0}\dfrac{x^t}{t!}\sum\limits_{i=1}^k i^t。

交换求和号得到 \sum\limits_{i=1}^k\sum\limits_{t\ge 0}\dfrac{(ix)^t}{t!}。

发现求和号后是 e^{ix}，使用等比数列求和得到最终的生成函数 \dfrac{e^{(k+1)x}-e^x}{e^x-1}。

上多项式全家桶即可 O(n\log n)/O(n\log^2 n) 取决于你的 \exp 写的什么。

代码

#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 File(a) freopen(#a".in","r",stdin);freopen(#a".out","w",stdout)
#define ll long long
using namespace std;
const int MAXN=(1<<22)+5,MOD=998244353,G=3;
inline ll qpow(ll base,ll 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[MAXN];
inline void init(){
    powg[0]=1;
    ll qwq=qpow(G,(MOD-1)>>23);
    F(i,0,21){
        ll qaq=qpow(qwq,1<<(21-i));
        F(j,1<<i,(1<<(i+1))-1) powg[j]=qaq*powg[j-(1<<i)]%MOD;
    }
    return;
}
inline void DIF(int*poly,int len){
    for(int i=len>>1;i>=1;i>>=1) for(int j=0,t=0;j<len;j+=(i<<1),++t) F(k,0,i-1){
        ll x=poly[j|k],y=poly[i|j|k]*1ll*powg[t]%MOD;
        (poly[j|k]=x+y)>=MOD&&(poly[j|k]-=MOD);
        (poly[i|j|k]=x-y)<0&&(poly[i|j|k]+=MOD);
    }
    return;
}
inline void DIT(int*poly,int len){
    for(int i=1;i<len;i<<=1) for(int j=0,t=0;j<len;j+=(i<<1),++t) F(k,0,i-1){
        ll x=poly[j|k],y=poly[i|j|k];
        (poly[j|k]=x+y)>=MOD&&(poly[j|k]-=MOD);
        poly[i|j|k]=(MOD+x-y)*powg[t]%MOD;
    }
    ll invl=qpow(len,MOD-2);
    F(i,0,len-1) poly[i]=poly[i]*invl%MOD;
    reverse(poly+1,poly+len);
    return;
}
inline void Inv(int*coef,int*res,int len){
    res[0]=qpow(coef[0],MOD-2);
    for(int l=2;l<=len<<1;l<<=1){
        static int tmp[MAXN];
        int qwq=l<<1;
        memcpy(tmp,coef,sizeof(int)*l);
        memset(tmp+l,0,sizeof(int)*(qwq-l));
        DIF(tmp,qwq),DIF(res,qwq);
        F(i,0,qwq-1) (res[i]=(2-res[i]*1ll*tmp[i])%MOD*res[i]%MOD)<0&&(res[i]+=MOD);
        DIT(res,qwq);
        memset(res+l,0,sizeof(int)*(qwq-l));
    }
    F(i,len+1,len<<1) res[i]=0;
    return;
}
int fact[MAXN],finv[MAXN],inv[MAXN],ans[MAXN],res[MAXN];
inline void Exp(int l,int r){
    if(l+1==r){
        if(l) ans[l]=ans[l]*1ll*inv[l]%MOD;
        else ans[l]=1;
        return;
    }
    int mid=(l+r)>>1;
    Exp(l,mid);
    static int f[MAXN],g[MAXN],len;
    len=1<<(__lg(r-l)+1);
    F(i,0,mid-l-1) f[i]=ans[i+l];
    F(i,0,r-l-2) g[i]=res[i];
    F(i,mid-l,len-1) f[i]=0;
    F(i,r-l-1,len-1) g[i]=0;
    DIF(f,len),DIF(g,len);
    F(i,0,len-1) f[i]=f[i]*1ll*g[i]%MOD;
    DIT(f,len);
    F(i,mid,r-1) (ans[i]+=f[i-l-1])>=MOD&&(ans[i]-=MOD);
    Exp(mid,r);
}
int k,m,mo[MAXN],ch[MAXN];
signed main(){
    ios::sync_with_stdio(0);
    cin.tie(0),cout.tie(0);
    init();
    cin>>k>>m;
    m+=2;
    fact[0]=fact[1]=finv[0]=inv[1]=1;
    F(i,2,m) fact[i]=fact[i-1]*1ll*i%MOD,inv[i]=MOD-MOD/i*1ll*inv[MOD%i]%MOD;
    F(i,1,m) finv[i]=finv[i-1]*1ll*inv[i]%MOD;
    int l=1<<(__lg(2*m)+1);
    F(i,1,m) ch[i-1]=finv[i];
    Inv(ch,mo,m);
    ll pw=k+1;
    F(i,1,m) ch[i-1]=finv[i]*(pw-1)%MOD,pw=pw*(k+1)%MOD;
    DIF(ch,l),DIF(mo,l);
    F(i,0,l-1) ch[i]=ch[i]*1ll*mo[i]%MOD;
    DIT(ch,l);
    F(j,1,m) res[j]=((j&1)?1ll:MOD-1ll)*inv[j]%MOD*ch[j]%MOD*fact[j]%MOD;
    F(j,1,m) res[j-1]=res[j]*1ll*j%MOD;
    res[m]=0;
    Exp(0,m);
    F(i,1,m-2) cout<<ans[i]*1ll*fact[i]%MOD<<"\n";
    return 0;
}

Day 108

感觉很典的逻辑题。

2025/05/18 [PA 2022] Mędrcy

题意

有 n 个人和 m 条咒语，每条咒语恰好被 n-2 个人知道，第 i 条咒语不被 a_i,b_i 知道。对于自己知道的咒语，每个人都知道其他哪些人知道这条咒语。如果一个人知道自己存在不知道的咒语就会在当天晚上走掉。第一天早上大家都知道至少有一个人知道一条咒语，求前 k 天第一次有人走是在第几天以及走了哪些人。

思路

既然是 n-2 个，那把不知道的两个人 a_i,b_i 连边。

一个人如果什么都不知道，他第 1 天就会走。在图上对应着删掉这个点连着的边就没有边了。

如果 x 看到 y 应该在第一天走而 y 没有，那么 x,y 就会一起在第 2 天走。在图上对应着，删掉 y 连着的边，再删掉 x 连着的边，就没有边了。

类似的，如果删掉 x_1,x_2,\cdots,x_t 连着的所有边就没有边了，那么 x_1,\cdots,x_t 就会在第 t 天走。

由于是求最早走的，相当于是求这个图的所有最小点覆盖。如果最小点覆盖 >k 就在前 k 天不会走。

这显然不存在多项式做法。由于 k 很小，我们来研究如何搞出一个和 k 有关的搜索。

设现在最小点覆盖还有 t 个点的空位，时间复杂度为 T(t)，初始 t=k。

枚举度数最大的点在或不在最小点覆盖中，设其度数为 D，如果在就递归到规模 -1 的子问题，如果不在这 D 个相邻点必选递归到规模 -D 的子问题。当 D=2 时只有环和链，环所有点都有可能在最小点覆盖，长度为奇数的链只有从头到尾偶数位置的在，长度为偶数的链都有可能在。特判 D\le 2，则有 T(t)=T(t-1)+T(t-3)。

总复杂度 O(T(k)(n+m))，可以通过。

代码

#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=1005;
int n,m,k,deg[MAXN]; 
vector<int>g[MAXN];
bitset<MAXN>ans,isin;
inline void mdf(int now,int v){for(int i:g[now]) deg[i]+=v;}
bool reach[MAXN],del[MAXN];
vector<int>p[MAXN];
void find(int now,int bg){
    reach[now]=1;p[bg].push_back(now);
    for(int i:g[now]) if(!reach[i]&&!del[i]) find(i,bg);
}
int mn,now;
void dfs(){
    if(now>mn) return;
    int qwq=0;
    F(i,1,n) if(!del[i]&&deg[i]>deg[qwq]) qwq=i;
    if(!qwq){
        if(now<mn) ans.reset(),mn=now;
        if(now==mn) ans|=isin;
        return;
    }else if(deg[qwq]>2){
        ++now,isin[qwq]=del[qwq]=1,mdf(qwq,-1);
        dfs();
        isin[qwq]=0,--now;
        vector<int>ne;
        for(int i:g[qwq]) if(!del[i]) ne.push_back(i),isin[i]=del[i]=1,mdf(i,-1),++now;
        dfs();
        for(int i:ne) isin[i]=del[i]=0,mdf(i,1),--now;
        del[qwq]=0,mdf(qwq,1);
    }else{
        memset(reach,0,n+1);
        vector<int>chain,cycle;
        int delta=0;
        F(i,1,n) if(!reach[i]&&!del[i]&&deg[i]<=1){
            vector<int>().swap(p[i]);
            find(i,i);
            chain.push_back(i);
            delta+=(p[i].size()>>1);
        }
        F(i,1,n) if(!reach[i]&&!del[i]){
            vector<int>().swap(p[i]);
            find(i,i);
            cycle.push_back(i);
            delta+=(p[i].size()-1)/2+1;
        }
        now+=delta;
        if(now<mn) ans.reset(),mn=now;
        if(now==mn){
            ans|=isin;
            for(int i:chain){
                if(p[i].size()&1){
                    bool flag=0;
                    for(int j:p[i]){
                        if(flag) ans[j]=1;
                        flag^=1;
                    }
                }else for(int j:p[i]) ans[j]=1;
            }
            for(int i:cycle) for(int j:p[i]) ans[j]=1;
        }
        now-=delta;
    }
}
signed main(){
    ios::sync_with_stdio(0);
    cin.tie(0),cout.tie(0);
    int T;
    for(cin>>T;T;--T){
        F(i,1,n) vector<int>().swap(g[i]);
        cin>>n>>m>>k,mn=k+1;
        unordered_set<int>vis;
        F(i,1,m){
            int u,v;
            cin>>u>>v;
            if(u>v) swap(u,v);
            if(vis.find(u*(n+1)+v)!=vis.end()) continue;
            else vis.insert(u*(n+1)+v);
            g[u].push_back(v),g[v].push_back(u);
        }
        F(i,1,n) deg[i]=g[i].size();
        ans.reset(); 
        dfs();
        if(mn==k+1) cout<<"-1\n";
        else{
            cout<<mn<<" "<<ans.count()<<"\n";
            F(i,1,n) if(ans[i]) cout<<i<<" ";
            cout<<"\n";
        }
    }
    return 0;
}

Day 109

每次看到新的模拟费用流题，都会觉得之前蜜月日记里总结的模拟费用流非常不典型。感觉之前蜜月日记里的全都只能是反悔贪心，没有用到费用流建模出的图，最多只用到了流量。所以放一道从建模的图来分析的模拟费用流。

2025/05/21 [Ynoi Easy Round 2018] 星野瑠美衣

题意

你有一个 n 个点的环，第 i 个点有坐标 (x_i,y_i)。你还有 m 个点，第 i 个坐标为 (x'_i,y'_i) 且有权值 v_i。你可以往每条边上插入至多 1 个这 m 个点之一。对于 i=1,2,\cdots,n，求恰好插入 i 个点时插入点的权值和加上环上相邻点的曼哈顿距离和的最大值。

思路

下面认为 x_{n+1}=x_1,y_{n+1}=y_1 以及 n,m 同阶。

先来研究 y 全部相等的情况。

数轴上我们把权值为 v 的 x 插入到 x_1,x_2 之间，则总权值有如下几种变化量：

1. 当 x<\min\{x_1,x_2\}，变化量为 -|x_1-x_2|+x_1+x_2-2x+v；

2. 当 \min\{x_1,x_2\}\le x\le \max\{x_1,x_2\}，变化量为 v；

3. 当 x>\max\{x_1,x_2\}，变化量为 -|x_1-x_2|-x_1-x_2+2x+v。

画图发现，我们可以不管前提条件，直接对这三个式子取 \max 就可以得到最大值。

再加上 y 的贡献，一共就是 9 种情况。

那么就可以进行费用流建模了。

我们一共有五类点：源点 S、汇点 T、代表 m 个点的 A_1,\cdots,A_m、代表可插入的 n 个空位的 B_1,\cdots,B_n、计算权值的虚点 C_1,\cdots,C_9。

有如下四类容量均为 1 的边：

点数边数都是 O(n)，没有正环。直接跑最大费用最大流即可，时间复杂度 O(n^3)。

直接费用流跑不过去，考虑模拟费用流，我们每次找一条增广路。

这张图是二分图，左部点有 S,T,C_1,\cdots,C_9 一共 11 个，剩下的是右部点。增广时我们会找一条从 S 到 T 的路径，由于是二分图，路径必然是左部点右部点交错，因此长度是 O(1) 的。

那么直接维护任意两个左部点之间经过一个右部点的最长路，以这个为边权在左部点内部跑 SPFA，找出来增广路后暴力修改影响到的路径。我们维护任意两个左部点之间经过一个右部点的所有路径，每次增广是把 O(1) 条边反转了，相当于是删除 O(1) 条路径再加入 O(1) 条，使用可删堆维护即可。

时间复杂度 O(n\log 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 fi first
#define se second
using namespace std;
const int MAXN=2e5+1;
const ll INF=0x3f3f3f3f3f3f3f3f;
int n,m,px[MAXN],py[MAXN],x[MAXN],y[MAXN],v[MAXN];
//S->0,T->10
priority_queue<pair<ll,int> >q[11][11],del[11][11];
ll ans,val[11][MAXN],dis[11][11],dep[11];
int dir[11][MAXN],fr[11];
inline void era(int u,int mid,int v){
    F(i,0,10) if(dir[i][mid]*dir[u][mid]==-1){
        if(dir[u][mid]==1) del[u][i].emplace(val[u][mid]+val[i][mid],mid);
        else del[i][u].emplace(val[u][mid]+val[i][mid],mid);
    }
    F(i,0,10) if(i!=u&&dir[i][mid]*dir[v][mid]==-1){
        if(dir[v][mid]==-1) del[i][v].emplace(val[v][mid]+val[i][mid],mid);
        else del[v][i].emplace(val[v][mid]+val[i][mid],mid);
    }
}
inline void add(int u,int mid,int v){
    F(i,0,10) if(dir[i][mid]*dir[v][mid]==-1){
        if(dir[v][mid]==1) q[v][i].emplace(val[v][mid]+val[i][mid],mid);
        else q[i][v].emplace(val[v][mid]+val[i][mid],mid);
    }
    F(i,0,10) if(i!=v&&dir[i][mid]*dir[u][mid]==-1){
        if(dir[u][mid]==-1) q[i][u].emplace(val[u][mid]+val[i][mid],mid);
        else q[u][i].emplace(val[u][mid]+val[i][mid],mid);
    }
}
signed main(){
    ios::sync_with_stdio(0);
    cin.tie(0),cout.tie(0);
    cin>>n>>m;
    F(i,1,n) cin>>px[i]>>py[i];
    px[n+1]=px[1],py[n+1]=py[1];
    F(i,1,m) cin>>x[i]>>y[i]>>v[i];
    F(i,1,m){
        val[0][i]=v[i],dir[0][i]=1;
        val[10][i]=-INF;
        F(j,1,9){
            dir[j][i]=-1;
            switch((j-1)/3){
                case 0:{val[j][i]-=2*x[i];break;}
                case 2:{val[j][i]+=2*x[i];break;}
            }
            switch(j%3){
                case 1:{val[j][i]-=2*y[i];break;}
                case 0:{val[j][i]+=2*y[i];break;}
            }
        }
    }
    F(i,1,n){
        val[0][i+m]=-INF;
        val[10][i+m]=-abs(px[i]-px[i+1])-abs(py[i]-py[i+1]),dir[10][i+m]=-1;
        F(j,1,9){
            dir[j][i+m]=1;
            switch((j-1)/3){
                case 0:{val[j][i+m]+=px[i]+px[i+1];break;}
                case 1:{val[j][i+m]+=abs(px[i]-px[i+1]);break;}
                case 2:{val[j][i+m]+=-px[i]-px[i+1];break;}
            }
            switch(j%3){
                case 1:{val[j][i+m]+=py[i]+py[i+1];break;}
                case 2:{val[j][i+m]+=abs(py[i]-py[i+1]);break;}
                case 0:{val[j][i+m]+=-py[i]-py[i+1];break;}
            }
        }
    }
    F(i,1,m) F(j,1,9) q[0][j].emplace(val[0][i]+val[j][i],i);
    F(i,1,n) F(j,1,9) q[j][10].emplace(val[10][i+m]+val[j][i+m],i+m);
    F(i,1,n) ans+=abs(px[i]-px[i+1])+abs(py[i]-py[i+1]);
    F(k,1,n){
        F(i,0,10) F(j,0,10){
            while(!q[i][j].empty()&&!del[i][j].empty()&&q[i][j].top()==del[i][j].top()) q[i][j].pop(),del[i][j].pop();
            if(q[i][j].empty()) dis[i][j]=-INF;
            else dis[i][j]=q[i][j].top().fi;
        }
        memset(dep,-0x3f,sizeof(dep));
        static bool isin[11];
        queue<int>qu;
        qu.push(0);
        dep[0]=0,isin[0]=1;
        while(!qu.empty()){
            int now=qu.front();
            qu.pop();isin[now]=0;
            F(i,0,10) if(dep[i]<dep[now]+dis[now][i]){
                dep[i]=dep[now]+dis[now][i],fr[i]=now;
                if(!isin[i]) isin[i]=1,qu.push(i);
            }
        }
        ans+=dep[10];
        cout<<ans<<" ";
        int now=10,tp=0;
        static int u[11],v[11],mid[11];
        while(now){
            int i=fr[now],t=q[i][now].top().se;
            era(i,t,now);
            ++tp,u[tp]=i,mid[tp]=t,v[tp]=now;
            now=i;
        }
        R(i,tp,1) dir[u[i]][mid[i]]*=-1,dir[v[i]][mid[i]]*=-1,val[u[i]][mid[i]]*=-1,val[v[i]][mid[i]]*=-1;
        while(tp) add(u[tp],mid[tp],v[tp]),--tp;
    }
    return 0;
}

Day 110

还有几道难炸了的集训队互测，先放一道没有炸的。

2025/05/27 [集训队互测 2024] 生命的循环

题意

给定一张 n 个点 m 条边的有向图，边有边权。每个点可以是激发或不激发。若有一条边从 u 指向 v 边权为 w 且 u 在时刻 t 激发，则 v 在时刻 t+w 激发。若没有其他点在下一时刻激发自己，则一个激发的点在下一时刻会变成不激发。0 时刻时节点 1 被激发，求充分长的时间后节点 n 激发情况的周期对 10^9+9 取模的值。保证所有简单环的边权和不超过 B。

思路

对于强连通的情况，周期显然为所有环环长的 \gcd。这是经典问题，做法是拉出一棵 dfs 树，求出每个点的深度 d_i，对于所有非树边 (u,v,w) 求 w+d_u-d_v 的 \gcd 即可。

对于其余情况，首先先对强连通分量缩点，把第 i 个 SCC 的周期 p_i 求出来，从一个 SCC 走到另一个 SCC 的边的新边权也可以用两个点在两个强连通分量的深度求出来。

我们预处理 e'_{i,j,k} 表示从 1 所在 SCC 到第 i 个 SCC，路径长度模 j 余 k 的路径是否存在，转移枚举边加入即可。

然后再预处理 e_{i,j,k} 表示从 1 所在 SCC 到第 i 个 SCC，路径的周期为 j，路径长度模 j 余 k 的路径是否存在。发现一条路径的周期就是经过所有的 SCC 的周期 \gcd，所以也是好转移的。

全部预处理出来之后图论的部分就结束了。把 e_{*,i} 视为一个字符串，你发现单个 i 的答案为 e_{*,i} 全部并起来无限复制的最小周期。如果所有 i 两两互质最终答案就是每个 i 答案的 \operatorname{lcm}。

如果不互质就可以用这个题的 Trick。我们枚举长度模 M=2520 的余数，则原本周期为 i 的新周期会变成 i/\gcd(i,M)，在 M=2520 时全部为质数的幂，要么两两互质，要么互为倍数。互质的求 \operatorname{lcm}，倍数的把短的复制几遍直到一样长。

最后就是怎么求并起来之后的周期的问题了。你发现并不好求，所以取反变成交，相当于是一个和上面说的那道题类似的同余方程组，用相同的方法做就行。

代码

#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)
#define fi first
#define se second
using namespace std;
const int MAXN=5001,MAXM=252001,M=2520;
int n,m,type,scc[MAXN],nn,per[MAXN],deg[MAXN];
vector<pair<int,int> >gs[MAXN];
bool exi[MAXN][101][101],s[MAXN][MAXM];
int bor[MAXM];
namespace SCC{
    vector<pair<int,int> >g[MAXN],gt[MAXN];
    int stk[MAXN],tp,dis[MAXN];
    bool vis[MAXN];
    void dfs1(int now){vis[now]=1;for(auto i:g[now])if(!vis[i.fi])dfs1(i.fi);stk[++tp]=now;}
    void dfs2(int now,int id){scc[now]=id;for(auto i:gt[now])if(!scc[i.fi])dfs2(i.fi,id);}
    void dfs3(int now){vis[now]=1;for(auto i:g[now])if(!vis[i.fi]&&scc[i.fi]==scc[now])dis[i.fi]=dis[now]+i.se,dfs3(i.fi);};
    inline void main(){
        F(i,1,n) if(!vis[i]) dfs1(i);
        memset(vis,0,n+1);
        while(tp){
            if(!scc[stk[tp]]) dfs2(stk[tp],++nn);
            --tp;
        }
        F(i,1,n) if(!vis[i]) dfs3(i);
        F(i,1,n) for(auto j:g[i]){
            if(scc[i]==scc[j.fi]) per[scc[i]]=__gcd(per[scc[i]],abs(j.se-dis[j.fi]+dis[i]));    
            else gs[scc[i]].emplace_back(scc[j.fi],j.se-dis[j.fi]+dis[i]),++deg[scc[j.fi]];
        }
    }
}
inline bool isp(int x){for(int i=2;i*i<=x;++i)if(x%i==0)return 0;return 1;}
int id[101],cnt,len[26],e[101];
inline void solve(int rem){ 
    bool avai[26][101]={};
    F(i,2,100){
        int qwq=__gcd(i,M);
        if(qwq==i){
            if(exi[scc[n]][i][rem%i]) return;
            continue;
        }
        F(j,0,len[id[i/qwq]]-1) avai[id[i/qwq]][j]|=exi[scc[n]][i][(j*M+rem)%i];
    }
    F(i,1,25){
        bool flag=0;
        F(j,0,len[i]-1) if(!avai[i][j]) flag=1;
        if(!flag) return;
    }
    F(i,1,25) F(j,0,len[i]-1) if(!avai[i][j]) s[i][j*M+rem]=1;
}
signed main(){
    ios::sync_with_stdio(0);
    cin.tie(0),cout.tie(0);
    cin>>n>>m>>type;
    F(i,1,m){
        int u,v,w;
        cin>>u>>v>>w;
        SCC::g[u].emplace_back(v,w);SCC::gt[v].emplace_back(u,w);
    }
    SCC::main();
    queue<int>q;
    F(i,1,100) exi[scc[1]][i][0]=1;
    q.push(scc[1]);
    while(!q.empty()){
        int now=q.front();q.pop();
        for(auto i:gs[now]){
            F(j,1,100) F(r,0,j-1) if(exi[now][j][r]) exi[i.fi][j][((r+i.se)%j+j)%j]=1;
            !(--deg[i.fi])&&(q.push(i.fi),1);
        }
    }
    F(i,1,nn) for(auto j:gs[i]) ++deg[j.fi];
    F(i,1,nn) F(j,1,100) if(j!=per[i]) F(k,0,j-1) exi[i][j][k]=0;
    q.push(scc[1]);
    while(!q.empty()){
        int now=q.front();q.pop();
        for(auto i:gs[now]){
            F(j,1,100) F(r,0,j-1) if(exi[now][j][r]){
                int qwq=__gcd(j,per[i.fi]);
                exi[i.fi][qwq][((r+i.se)%qwq+qwq)%qwq]=1;
            }
            !(--deg[i.fi])&&(q.push(i.fi),1);
        }
    }
    F(i,2,100) if(isp(i)){
        int qwq=1;
        ++cnt;
        while(qwq*i<=100) qwq*=i,id[qwq]=cnt;
        len[id[qwq]]=qwq;
    }
    F(rem,0,M-1) solve(rem);
    F(t,1,25){
        static bool ch[MAXM];
        int qwq=len[t]*M,qaq=qwq;
        F(j,1,qwq) ch[j]=s[t][j-1];
        for(int i(2),j(0);i<=qwq;++i){
            while(j&&ch[i]!=ch[j+1]) j=bor[j];
            if(ch[i]==ch[j+1]) ++j;
            bor[i]=j;
        }
        if(qwq%(qwq-bor[qwq])==0) qaq=qwq-bor[qwq];
        for(int j=2;j*j<=qaq;++j) if(qaq%j==0){
            int ee=0;
            while(qaq%j==0) qaq/=j,++ee;
            e[j]=max(e[j],ee);
        }
        if(qaq!=1) e[qaq]=max(e[qaq],1);
    }
    ll ans=1;
    F(i,2,100) while(e[i]) ans=ans*i%1000000009,--e[i];
    cout<<ans;
    return 0;
}

后记

最近有很多好题啊，慢慢往蜜月日记里面放吧。