結果
| 問題 |
No.3182 recurrence relation’s intersection sum
|
| コンテスト | |
| ユーザー |
 Today03
|
| 提出日時 | 2025-11-07 08:28:21 |
| 言語 | C++23 (gcc 13.3.0 + boost 1.87.0) |
| 結果 |
AC
|
| 実行時間 | 8 ms / 2,000 ms |
| コード長 | 18,374 bytes |
| コンパイル時間 | 4,299 ms |
| コンパイル使用メモリ | 311,252 KB |
| 実行使用メモリ | 7,724 KB |
| 最終ジャッジ日時 | 2025-11-07 08:28:27 |
| 合計ジャッジ時間 | 6,122 ms |
|
ジャッジサーバーID （参考情報） |
judge3 / judge2 |
（要ログイン）
| ファイルパターン | 結果 |
|---|---|
| sample | AC * 3 |
| other | AC * 40 |
ソースコード
#include<bits/stdc++.h>
using namespace std;
#define rep(i,n) for(int i=0; i<n; i++)
#define all(x) (x).begin(), (x).end()
bool chmax(auto& a, auto b) { return a<b ? a=b, true : false; }
bool chmin(auto& a, auto b) { return a>b ? a=b, true : false; }
using ll=long long; const int INF=1e9+10; const ll INFL=4e18;
#ifdef DEBUG
#include "./debug.hpp"
#else
#define debug(...)
#define print_line
#endif
/// @brief ModInt
template<ll MOD>
struct ModInt {
ModInt(ll x=0){ value=(x>=0?x%MOD:MOD-(-x)%MOD); }
ModInt operator-() const { return ModInt(-value); }
ModInt operator+() const { return ModInt(*this); }
ModInt& operator+=(const ModInt& other) {
value+=other.value;
if(value>=MOD) value-=MOD;
return*this;
}
ModInt& operator-=(const ModInt& other) {
value+=MOD-other.value;
if(value>=MOD) value-=MOD;
return*this;
}
ModInt& operator*=(const ModInt other) {
value=value*other.value%MOD;
return*this;
}
ModInt& operator/=(ModInt other) {
(*this)*=other.inv();
return*this;
}
ModInt operator+(const ModInt& other) const { return ModInt(*this)+=other; }
ModInt operator-(const ModInt& other) const { return ModInt(*this)-=other; }
ModInt operator*(const ModInt& other) const { return ModInt(*this)*=other; }
ModInt operator/(const ModInt& other) const { return ModInt(*this)/=other; }
bool operator==(const ModInt& other) const { return value==other.value; }
bool operator!=(const ModInt& other) const { return value!=other.value; }
friend ostream& operator<<(ostream& os, const ModInt& x) { return os<<x.value; }
friend istream& operator>>(istream& is, ModInt& x) {
ll v;
is>>v;
x=ModInt<MOD>(v);
return is;
}
ModInt pow(ll x) const {
ModInt ret(1),mul(value);
while(x) {
if(x&1) ret*=mul;
mul*=mul;
x>>=1;
}
return ret;
}
ModInt inv() const { return pow(MOD-2); }
ll val() {return value; }
static constexpr ll mod() { return MOD; }
private:
ll value;
};
using Mod998=ModInt<998244353>;
using Mod107=ModInt<1000000007>;
/// @brief 数論変換
/// @note O(N log(N))
/// @details f(x) = Σ a[i]x^i, w^N = 1 とすると、F(t) = Σ f(w^i)t^i の各係数 mod 998244353 に変換する
void NTT998(vector<Mod998>& a, bool inv=false) {
static bool flag=false;
static int divide_max;
static vector<Mod998> roots, inv_roots, tmp;
if(!flag) {
flag=true;
divide_max=1;
ll n=998244353-1;
while(n%2==0) n>>=1,divide_max++;
roots=vector<Mod998>(divide_max+1);
inv_roots=vector<Mod998>(divide_max+1);
roots[0]=inv_roots[0]=1;
for(int i=1; i<=divide_max; i++) {
roots[i]=Mod998(3).pow((998244353-1)/(1<<i));
inv_roots[i]=roots[i].inv();
}
}
tmp=vector<Mod998>(a.size());
int n=a.size(), mask=n-1, p=0;
for(int i=n>>1; i>=1; i>>=1) {
auto& cur=((p&1) ? tmp : a);
auto& nxt=((p&1) ? a : tmp);
Mod998 e=(inv ? roots[p+1] : inv_roots[p+1]);
Mod998 w=1;
for(int j=0; j<n; j+=i) {
rep(k,i) nxt[j+k]=cur[((j<<1)&mask)+k]+w*cur[(((j<<1)+i)&mask)+k];
w*=e;
}
p++;
}
if(p&1) swap(a,tmp);
if(inv) {
Mod998 inv_sz=Mod998(n).inv();
for(int i=0; i<n; i++) a[i]*=inv_sz;
}
}
/// @brief AとBの畳み込み C[i] = Σ A[j]B[i-j] mod 998244353 を返す
/// @note O(N log(N))
/// @attention |a|+|b| <= 2^23 が必要
vector<Mod998> Convolve998(vector<Mod998> a,vector<Mod998> b) {
int n=1;
while(n+1<a.size()+b.size()) n<<=1;
vector<Mod998> fa(n), fb(n);
for(int i=0; i<a.size(); i++) fa[i]=a[i];
for(int i=0; i<b.size(); i++) fb[i]=b[i];
NTT998(fa); NTT998(fb);
for(int i=0; i<n; i++) fa[i]*=fb[i];
NTT998(fa,true);
while(fa.size()+1>a.size()+b.size()) fa.pop_back();
return fa;
}
using Fps=vector<Mod998>; //< 多項式
using FpsSparse=vector<pair<int,Mod998>>; //< 疎な多項式（(次数, 係数)の配列）
/// @brief 形式的冪級数
/// @ref https://potato167.github.io/po167_library
namespace FPS {
/// @brief 多項式 f, g の和を返す
Fps Add(const Fps& a, const Fps& b) {
Fps res(max(a.size(),b.size()));
for(int i=0; i<res.size(); i++) {
if(i<a.size()) res[i]+=a[i];
if(i<b.size()) res[i]+=b[i];
}
return res;
}
/// @brief 多項式 f, g の差を返す
Fps Sub(const Fps& f, const Fps& g, int len=-1) {
if(len==-1) len=max(f.size(),g.size());
Fps res(len);
for(int i=0; i<len; i++) {
if(i<f.size()) res[i]+=f[i];
if(i<g.size()) res[i]-=g[i];
}
return res;
}
/// @brief 多項式 f, g の積を返す
Fps Mul(const Fps& f, const Fps& g, int len=-1) {
auto fg=Convolve998(f,g);
if(len!=-1) fg.resize(len+1);
return fg;
}
/// @brief 多項式 f に対し、f*g = 1 なる g を返す
/// @note O(N log(N))
/// @ref https://judge.yosupo.jp/problem/inv_of_formal_power_series
Fps Inv(Fps f, int len=-1) {
/**
* 方程式 h(g) = 1/g - f = 0 の解を求める問題に変換する
* ニュートン法より
* g_{n+1} = g_n - h(g_n)/h'(g_n)
* h'(g) = - 1/g^2 より
* g_{n+1} = g_n - (1/g_n - f) / (-1/g_n^2)
* = g_n + g_n - f*g_n^2
* = g_n * (2 - f*g_n)
* ここで、f*g_n を上位次数の項と下位次数の項 fg_low, fg_high に分けると、
* fg_low + fg_high = f*g_n, fg_low = 1 (mod x^d)
* よって、
* g_{n+1} = g_n * (2 - fg_low - fg_high)
* = g_n * (1 - fg_high)
* = g_n - g_n*fg_high
*/
if(len==-1) len=f.size();
assert(f[0]!=0);
Fps g={f[0].inv()};
int s=1;
while(s<len) {
Fps nxtg(s*2,0), res(s*2,0);
g.resize(s*2);
for(int i=0; i<s*2; i++) {
if(int(f.size())>i) res[i]=f[i];
nxtg[i]=g[i];
}
//fg_high を計算
NTT998(g); NTT998(res);
for(int i=0; i<s*2; i++) res[i]*=g[i];
NTT998(res,true); for(int i=0; i<s; i++) res[i]=0; //fg_high
//fg_high * g_n を計算
NTT998(res); for(int i=0; i<s*2; i++) res[i]*=g[i]; NTT998(res,true);
for(int i=s; i<s*2; i++) nxtg[i]-=res[i];
swap(nxtg,g);
s*=2;
}
g.resize(len);
return g;
}
/// @brief 多項式 f の積分を返す
Fps Integral(Fps f) {
if(f.empty()) return f;
Fps num_inv((int)f.size()+1);
num_inv[0]=1; num_inv[1]=1;
auto m=Mod998::mod();
for(int i=2; i<=(int)f.size(); i++) num_inv[i]=(-num_inv[m%i])*(Mod998)(m/i);
f.reserve((int)f.size()+1);
f.push_back(0);
for(int i=(int)f.size()-1; i>0; i--) f[i]=f[i-1]*num_inv[i];
f[0]=0;
return f;
}
/// @brief 多項式 f の微分を返す
Fps Differential(Fps f) {
if(f.empty()) return f;
for(int i=0; i<(int)f.size()-1; i++) f[i]=f[i+1]*(Mod998)(i+1);
f.pop_back();
return f;
}
/// @brief 多項式 f, g について、`f = gq + r` なる q, r を返す
pair<Fps,Fps> Div(Fps f, Fps g) {
int n=f.size(),m=g.size();
if(n<m) return{{},f};
Fps r=f;
reverse(all(f)); reverse(all(g));
f.resize(n-m+1);
Fps q=Mul(f,Inv(g,n-m+1));
q.resize(n-m+1);
reverse(all(q)); reverse(all(g));
r=Sub(r,Mul(q,g));
while(!q.empty()&&q.back()==0) q.pop_back();
while(!r.empty()&&r.back()==0) r.pop_back();
return {q,r};
}
/// @brief 多項式 f, g の積を返す（ただし、g は疎な多項式として与える）
Fps MulSparse(Fps f, FpsSparse g) {
auto itr=find_if(g.begin(),g.end(),[&](auto p) { return p.first==0; });
Mod998 x0=0;
if(itr!=g.end()) {
x0=itr->second;
g.erase(itr);
}
for(int i=(int)f.size()-1; i>=0; i--) {
for(auto& [d,c]: g) {
if(i+d>=f.size()) continue;
f[i+d]+=f[i]*c;
}
f[i]*=x0;
}
return f;
}
/// @brief 多項式 f, g に対し、 f / g を返す（ただし、g は疎な多項式として与える）
Fps DivSparse(Fps f, FpsSparse g) {
auto itr=find_if(g.begin(),g.end(),[&](auto p) { return p.first==0; });
assert(itr!=g.end());
Mod998 x0_inv=itr->second.inv();
g.erase(itr);
for(int i=0; i<f.size(); i++) {
f[i]*=x0_inv;
for(auto& [d,c]: g) {
if(i+d>=f.size()) continue;
f[i+d]-=f[i]*c;
}
}
return f;
}
namespace Internal {
const int PRIMITIVE_ROOT=3;
Fps CyclicConvolution(Fps f, Fps g) {
NTT998(f); NTT998(g);
for(int i=0; i<(int)f.size(); i++) f[i]*=g[i];
NTT998(f,true);
return f;
}
//in :DFT(v)(len(v)=z)
//out:DFT(v)(len(v)=2*z)
void Extend(Fps& v) {
int z=v.size();
Mod998 e=Mod998(PRIMITIVE_ROOT).pow(Mod998::mod()/(2*z));
auto cp=v;
NTT998(cp,true);
Mod998 tmp=1;
for(int i=0; i<z; i++) {
cp[i]*=tmp;
tmp*=e;
}
NTT998(cp);
for(int i=0; i<z; i++) v.push_back(cp[i]);
}
//s.t|v|=2^s(no assert)
void PickEvenOdd(Fps& v, int odd) {
int z=v.size()/2;
Mod998 half=(Mod998)(1)/(Mod998)(2);
if(odd==0) {
for(int i=0; i<z; i++) v[i]=(v[i*2]+v[i*2+1])*half;
v.resize(z);
}else{
Mod998 e=Mod998(PRIMITIVE_ROOT).pow(Mod998::mod()/(2*z));
Mod998 ie=Mod998(1)/e;
Fps es={half};
while((int)es.size()!=z) {
Fps n_es((int)es.size()*2);
for(int i=0; i<(int)es.size(); i++) {
n_es[i*2]=(es[i]);
n_es[i*2+1]=(es[i]*ie);
}
ie*=ie;
swap(n_es,es);
}
for(int i=0; i<z; i++) v[i]=(v[i*2]-v[i*2+1])*es[i];
v.resize(z);
}
}
}
/// @brief 多項式 f について、e^f = Σ[k=0~len-1](f(x)^k/k!) を返す
Fps Exp(Fps f, int len=-1) {
/**
* 方程式 h(g) = log(g) - f = 0 の解を求める問題に変換する。
* ニュートン法より
* g_{n+1} = g_n - h(g_n)/h'(g_n)
* h(g) = log(g) - f, h'(g) = 1/g より
* g_{n+1} = g_n - (log(g_n) - f) / (1/g_n)
* = g_n - g_n(log(g_n) - f)
* = g_n * (1 - log(g_n) + f)
* ここで、log(g_n) を上位項と下位項 lg_low, lg_high に分けると、
* log(g_n) = lg_low + lg_high, lg_low = f_low
* よって、
* g_{n+1} = g_n(1 - f_low - lg_high + f)
* = g_n(1 + (f - log(g))_high)
* = g_n + g_n * (f - log(g))_high
*/
if(len==-1) len=f.size();
if(len==0) return{};
if(len==1) return{Mod998(1)};
assert(!f.empty() && f[0]==0);
int s=1;
Fps g={Mod998(1)};
while(s<len) {
Fps A=g,B=g;
A=Differential(A); //A = g'
B=Inv(B,2*s); //B = 1/g
A.resize(2*s);
A=Internal::CyclicConvolution(A,B); //A = g'/g
A.pop_back();
A=Integral(A); //A = ∫(g'/g)dx = log(g)
//A = (f-log(g_n))_high
for(int i=0; i<s; i++) A[i]=0;
for(int i=s; i<s*2; i++) A[i]=(i<(int)f.size() ? f[i] : 0)-A[i];
//g_{n+1} = g_n + g_n * (f-log(g))_high
g.resize(2*s);
B=Internal::CyclicConvolution(A,g);
for(int i=s; i<s*2; i++) g[i]=B[i];
s*=2;
}
g.resize(len);
return g;
}
/// @brief 多項式 f について、log(f) を返す
Fps Log(Fps f, int len=-1) {
if(len==-1) len=f.size();
if(len==0) return{};
if(len==1) return{Mod998(0)};
assert(!f.empty()&&f[0]==1);
Fps res=Convolve998(Differential(f),Inv(f,len));
res.resize(len-1);
return Integral(res);
}
/// @brief 多項式 f^M を返す
Fps Pow(Fps f, ll M, int len=-1) {
if(len==-1) len=f.size();
Fps res(len,0);
if(M==0) {
res[0]=1;
return res;
}
for(int i=0; i<(int)f.size(); i++) {
if(f[i]==0) continue;
if(i>(len-1)/M) break;
Fps g((int)f.size()-i);
Mod998 v=(Mod998)(1)/(Mod998)(f[i]);
for(int j=i; j<(int)f.size(); j++) g[j-i]=f[j]*v;
ll zero=i*M;
if(i) len-=i*M;
g=Log(g,len);
for(Mod998& x:g) x*=M;
g=Exp(g,len);
v=(Mod998)(1)/v;
Mod998 c=1;
while(M) {
if(M&1) c=c*v;
v=v*v;
M>>=1;
}
for(int i=0; i<len; i++) res[i+zero]=g[i]*c;
return res;
}
return res;
}
/// @brief `[x^k](P/Q)` を返す
Mod998 BostanMori(ll k, Fps P, Fps Q) {
assert(!Q.empty()&&Q[0]!=0);
int z=1;
while(z<(int)max(P.size(),Q.size())) z*=2;
P.resize(z*2,0); Q.resize(z*2,0);
NTT998(P); NTT998(Q);
//fast
while(k) {
//Q(-x)
Fps Q_n(z*2);
for(int i=0; i<z; i++) {
Q_n[i*2]=Q[i*2+1];
Q_n[i*2+1]=Q[i*2];
}
for(int i=0; i<z*2; i++) {
P[i]*=Q_n[i];
Q[i]*=Q_n[i];
}
Internal::PickEvenOdd(P,k&1);
Internal::PickEvenOdd(Q,0);
k/=2;
if(k==0) break;
Internal::Extend(P);
Internal::Extend(Q);
}
Mod998 SP=0,SQ=0;
for(int i=0; i<z; i++) SP+=P[i],SQ+=Q[i];
return SP/SQ;
}
//0=a[i]*c[0]+a[i-1]*c[1]+a[i-2]*c[2]+...+a[i-d]*c[d]
//a.size()+1==c.size()
//c[0]=-1?
//return a[k]
Mod998 KthLinear(ll k, Fps a, Fps c) {
int d=a.size();
assert(d+1==int(c.size()));
Fps P=Convolve998(a,c);
P.resize(d);
return BostanMori(k,P,c);
}
}
vector<Mod998> BerlekampMassey(const vector<Mod998>& s) {
int n=s.size();
vector<Mod998> c{1}, b{1};
int l=0, m=1;
Mod998 bb=1;
for(int i=0; i<n; i++) {
Mod998 d=0;
for(int j=0; j<=l; j++) d+=c[j]*s[i-j];
if(d==Mod998(0)) m++;
else if(2*l<=i) {
auto t=c; Mod998 coef=d*bb.inv();
if(c.size()<b.size()+m) c.resize(b.size()+m,0);
for(int j=0; j<b.size(); j++) c[j+m]-=coef*b[j];
l=i+1-l; b=t; bb=d; m=1;
} else {
Mod998 coef=d*bb.inv();
if(c.size()<b.size()+m) c.resize(b.size()+m,0);
for(int j=0; j<b.size(); j++) c[j+m]-=coef*b[j];
m++;
}
}
if(!c.empty()) c.erase(c.begin());
for(auto& x: c) x=-x;
return c;
// const int N = (int)f.size(); vector<Mod998> b, c;
// b.reserve(N + 1); c.reserve(N + 1); b.push_back(Mod998(1)); c.push_back(Mod998(1));
// Mod998 y = Mod998(1);
// for (int ed = 1; ed <= N; ed++) {
// int l = int(c.size()), m = int(b.size()); Mod998 x = 0;
// for (int i = 0; i < l; i++) x += c[i] * f[ed - l + i];
// b.emplace_back(Mod998(0)); m++;
// if (x == Mod998(0)) continue;
// Mod998 freq = x / y;
// if (l < m) {
// auto tmp = c; c.insert(begin(c), m - l, Mod998(0));
// for (int i = 0; i < m; i++) c[m - 1 - i] -= freq * b[m - 1 - i];
// b = tmp; y = x;
// } else {
// for (int i = 0; i < m; i++) c[l - 1 - i] -= freq * b[m - 1 - i];
// }
// }
// reverse(begin(c), end(c));
// return c;
}
Mod998 bostanMori(std::vector<Mod998> num, std::vector<Mod998> den, long long n) {
while (n > 0) {
// den_neg = den with odd coefficients negated
std::vector<Mod998> den_neg = den;
for (size_t i = 0; i < den_neg.size(); ++i) {
if (i % 2 == 1) den_neg[i] = Mod998(0) - den_neg[i];
}
// f2 = num * den_neg
// g2 = den * den_neg
std::vector<Mod998> f2 = Convolve998(num, den_neg);
std::vector<Mod998> g2 = Convolve998(den, den_neg);
std::vector<Mod998> num2, den2;
// pick even or odd coefficients depending on parity of n
if (n % 2 == 0) {
for (size_t i = 0; i < f2.size(); i += 2) num2.push_back(f2[i]);
} else {
for (size_t i = 1; i < f2.size(); i += 2) num2.push_back(f2[i]);
}
for (size_t i = 0; i < g2.size(); i += 2) den2.push_back(g2[i]);
num = std::move(num2);
den = std::move(den2);
n /= 2;
}
return num[0] * den[0].inv();
}
Mod998 LinearRecurrence(vector<Mod998> a, vector<Mod998> c, ll k) {
int n=c.size();
if(n==0) return 0;
vector<Mod998> dnm(n+1); dnm[0]=1;
rep(i,n) dnm[i+1]=-c[i];
a.resize(n,0);
auto num=Convolve998(dnm,a); if(num.size()>n) num.resize(n);
// return FPS::BostanMori(k,num,dnm);
return bostanMori(num,dnm,k);
}
Mod998 BMBM(const vector<Mod998> &s, ll n) {
auto bm=BerlekampMassey(s);
return LinearRecurrence(s,bm,n);
}
//----------------------------------------------------------
void solve() {
ll K,L,R; cin>>K>>L>>R;
ll N=500;
vector<Mod998> A(N+1);
A[0]=1; rep(i,N) A[i+1]=A[i]*K+Mod998(i).pow(K)+Mod998(K).pow(i);
vector<Mod998> B(N+2);
rep(i,N+1) B[i+1]=B[i]+A[i];
cout<<BMBM(B,R+1)-BMBM(B,L)<<endl;
}
int main() {
ios::sync_with_stdio(false); cin.tie(nullptr);
//cout<<fixed<<setprecision(15);
int T=1; //cin>>T;
while(T--) solve();
}