Thank you for participation and we hope you enjoy this round :)
D2A Submission is All You Need
Hint 1
Tutorial
solution
#include <bits/stdc++.h>
using namespace std;
void solve(){
int n, tmp, ans = 0;
cin >> n;
for(int i=0; i<n; i++){cin >> tmp; ans += tmp + (tmp == 0);}
cout << ans << '\n';
}
int main(){
ios::sync_with_stdio(false);
cin.tie(NULL);
int t;
cin >> t;
while(t--) solve();
}
D2B Pathless
Hint1
Hint2
Hint3
Tutorial
solution
#include <bits/stdc++.h>
using namespace std;
void solve(){
int n, k, tmp;
cin >> n >> k;
vector<int> a(3);
for(int i=0; i<n; i++){cin >> tmp; a[tmp]++;}
if(k == a[1] + 2*a[2] || k >= a[1] + 2*a[2] + 2){cout << -1 << '\n'; return;}
for(int i=0; i<a[0]; i++) cout << 0 << ' ';
for(int i=0; i<a[2]; i++) cout << 2 << ' ';
for(int i=0; i<a[1]; i++) cout << 1 << ' ';
cout << '\n';
}
int main(){
ios::sync_with_stdio(false);
cin.tie(NULL);
int t;
cin >> t;
while(t--) solve();
}
D1A Double Perspective
Hint1
Hint2
Tutorial
solution
#include <vector>
#include <list>
#include <map>
#include <set>
#include <deque>
#include <vector>
#include <list>
#include <map>
#include <set>
#include <deque>
#include <queue>
#include <stack>
#include <bitset>
#include <algorithm>
#include <functional>
#include <numeric>
#include <utility>
#include <sstream>
#include <iostream>
#include <iomanip>
#include <cstdio>
#include <cmath>
#include <cstdlib>
#include <cctype>
#include <string>
#include <cstring>
#include <ctime>
#include <random>
#include <chrono>
using namespace std;
#define _int64 long long
#define mo 998244353
pair<int,int> ma[610000];
int main()
{
int i,j,k,n,l,t,m,x,y,o,b1;
vector<int> ans;
scanf("%d",&t);
for (l=0;l<t;l++)
{
scanf("%d",&n);
for (i=0;i<n+n;i++)
{
ma[i]=make_pair(-1,-1);
}
for (i=0;i<n;i++)
{
scanf("%d%d",&x,&y);
x--;y--;
ma[x]=max(ma[x],make_pair(y,i));
}
ans.clear();
for (i=0;i<n+n;i++)
{
if (ma[i]==make_pair(-1,-1)) continue;
ans.push_back(ma[i].second);
}
printf("%d\n",(int)ans.size());
for (i=0;i<ans.size();i++)
printf("%d ",ans[i]+1);
printf("\n");
}
}
D1B Stay or Mirror
Hint1
Hint2
Tutorial
solution
#include <bits/stdc++.h>
using namespace std;
void solve(){
int n, ans = 0;
cin >> n;
vector<int> p(n);
for(int i=0; i<n; i++) cin >> p[i];
for(int i=0; i<n; i++){
int left = 0, right = 0;
for(int j=0; j<i; j++) if(p[j] > p[i]) left++;
for(int j=i+1; j<n; j++) if(p[j] > p[i]) right++;
ans += min(left, right);
}
cout << ans << '\n';
}
int main(){
ios::sync_with_stdio(false);
cin.tie(NULL);
int t;
cin >> t;
while(t--) solve();
}
D1C1 Interactive RBS (Easy Version), D1C2 Interactive RBS (Medium Version) and D1C3 Interactive RBS (Hard Version)
Hint1
Hint2
Hint3
Hint4
Hint5
Tutorial
solution (easy version)
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef long long ll;
typedef long double ld;
typedef pair<ll, ll> pll;
typedef pair<int, int> pii;
typedef vector<int> vi;
typedef vector<ll> vll;
const long long kk = 1000;
const long long ml = kk * kk;
const long long mod = ml * kk + 7;
const long long inf = ml * ml * ml + 7;
#define T(x, i) get<i>(x)
#define F first
#define S second
#define all(v) v.begin(), v.end()
#define forn(i, n) for (int i = 0; i < n; i++)
template<class T, class U> inline istream& operator>>(istream& str, pair<T, U> &p) { return str >> p.F >> p.S; }
template<class T, class U> inline ostream& operator<<(ostream& str, const pair<T, U> &p) { str << p.F << ' ' << p.S << ' '; return str; }
template<class T> inline istream& operator>>(istream& str, vector<T> &a) { for (auto &i : a) str >> i; return str; }
template<class T> inline ostream& operator<<(ostream& str, const vector<T> &a) { for (auto &i : a) str << i << ' '; return str; }
void flush() { cout << flush; }
void flushln() { cout << endl; }
template<class T> void read(T &x) { cin >> x; }
template<class T, class ...U> void read(T &x, U&... u) { read(x); read(u...); }
template<class T> void print(const T &x) { cout << x; }
template<class T, class ...U> void print(const T &x, const U&... u) { print(x); print(u...); }
template<class ...T> void println(const T&... u) { print(u..., '\n'); }
template<class T> void eprint(const T &x) { cerr << x; }
template<class T, class ...U> void eprint(const T &x, const U&... u) { eprint(x); eprint(u...); }
template<class ...T> void eprintln(const T&... u) { eprint(u..., '\n'); }
#ifdef DEBUG
mt19937 rng(1033);
#define var(x) cerr << #x << ": " << x << '\n';
#define range(a, b) cerr << #a << ", " << #b << ": "; for (auto _it = a; _it != b; ++_it) cerr << *_it << ' '; cerr << '\n';
#else
mt19937 rng(chrono::steady_clock::now().time_since_epoch().count());
#define cerr if (false) cerr
#define var(x)
#define range(a, b)
#define eprintln if(false)eprintln
#define eprint if(false)eprint
#endif
int rnd(int mod) { return uniform_int_distribution<int>(0, mod - 1)(rng); }
int n;
map<vector<int>, int> mem;
int ask(vector<int> ind) {
if (!mem.count(ind)) {
cout << "? ";
cout << ind.size() << ' ';
for (auto i : ind)
cout << i + 1 << ' ';
cout << endl;
int res;
cin >> res;
mem[ind] = res;
}
return mem[ind];
}
bool same(vector<int> ind) {
if (ask(ind))
return false;
ind.push_back(ind[0]);
ind.erase(ind.begin());
if (ask(ind))
return false;
return true;
}
pii search(vector<int> ind) {
if (ind.size() == 2) {
if (ask(ind))
return {ind[0], ind[1]};
return {ind[1], ind[0]};
}
vector<int> inda, indb;
int sz = (ind.size() + 1) / 2;
for (int i = 0; i < sz; i++) {
inda.push_back(ind[i]);
indb.push_back(ind[ind.size() - 1 - i]);
}
if (!same(inda))
return search(inda);
if (!same(indb))
return search(indb);
inda.push_back(indb[0]);
return search(inda);
}
pii find_reg2() {
vector<int> ind(n);
iota(all(ind), 0);
shuffle(all(ind), rng);
return search(ind);
}
void solve() {
mem.clear();
auto [po, pc] = find_reg2();
vector<int> to_ask(n);
iota(all(to_ask), 0);
if (to_ask.size() & 1)
to_ask.push_back(0);
string ans = string(n, '.');
ans[po] = '(';
ans[pc] = ')';
for (int i = 0; i < to_ask.size(); i += 2) {
int a = to_ask[i], b = to_ask[i + 1];
vector<int> q = {po, pc, po, a, b, pc};
int res = ask(q);
// ()(ab) - x ab
// ()(()) - 4 ()
// ()((() - 2 ((
// ()())) - 3 ))
// ()()() - 6 )(
map<int, pair<char, char>> results = {
{4, {'(', ')'}},
{2, {'(', '('}},
{3, {')', ')'}},
{6, {')', '('}},
};
assert(results.count(res));
auto [resa, resb] = results[res];
ans[a] = resa;
ans[b] = resb;
}
cout << "! " << ans << endl;
}
signed main() {
#ifdef DEBUG
// freopen("input.txt", "r", stdin);
// freopen("output.txt", "w", stdout);
#endif
ios_base::sync_with_stdio(0); cin.tie(0);
cout << fixed << setprecision(20);
int t = 1;
cin >> t;
for (int i = 0; i < t; i++) {
cin >> n;
solve();
}
#ifdef DEBUG
cerr << "Runtime is: " << clock() * 1.0 / CLOCKS_PER_SEC << endl;
#endif
}
/*
()(ab)
()(()) - 4 ()
()((() - 2 ((
()())) - 3 ))
()()() - 6 )(
2
3
)((
2
()
*/
solution (medium version)
#include <vector>
#include <list>
#include <map>
#include <set>
#include <deque>
#include <vector>
#include <list>
#include <map>
#include <set>
#include <deque>
#include <queue>
#include <stack>
#include <bitset>
#include <algorithm>
#include <functional>
#include <numeric>
#include <utility>
#include <sstream>
#include <iostream>
#include <iomanip>
#include <cstdio>
#include <cmath>
#include <cstdlib>
#include <cctype>
#include <string>
#include <cstring>
#include <ctime>
#include <random>
#include <chrono>
using namespace std;
#define _int64 long long
#define mo 998244353
char tar[110000];
int cnt;
int doask(vector<int>& a)
{
int i,ret;
/*
{
int v,tmp,j;
string ss;
cnt++;
ss="";
for (i=0;i<a.size();i++)
{
ss.push_back(tar[a[i]]);
}
tmp=0;
for (i=0;i<ss.length();i++)
{
v=0;
for (j=i;j<ss.length();j++)
{
if (ss[j]=='(') v++;
else v--;
if (v<0) break;
if (v==0) tmp++;
}
}
return tmp;
}
*/
printf("? %d",(int)a.size());
for (i=0;i<a.size();i++)
printf(" %d",a[i]+1);
printf("\n");
fflush(stdout);
scanf("%d",&ret);
return ret;
}
char ans[1100];
void doanswer()
{
printf("! %s\n",ans);
cerr<<"T "<<tar<<endl;
cerr<<"cnt:"<<cnt<<endl;
if (string(ans)==string(tar))
{
cerr<<"good"<<endl;
}
else cerr<<"bad"<<endl;
fflush(stdout);
}
int main()
{
int i,j,k,n,l,t,m,x,y,o,b1,indl,indr,ll,rr,mid,tmp,v,st;
vector<int> a;
scanf("%d",&t);
for (l=0;l<t;l++)
{
scanf("%d",&n);
for (i=0;i<n;i++)
{
x=rand()%2;
if (x==0) tar[i]='(';
else tar[i]=')';
}
tar[n]='\0';
cnt=0;
st=0;ans[n]='\0';
indl=-1;
indr=-1;
// while (st<n)
// {
ll=st;rr=n;
while (rr-ll>1)
{
mid=(ll+rr)/2;
a.clear();
for (i=st;i<=mid;i++)
a.push_back(i);
tmp=doask(a);
if (tmp==0) ll=mid;
else rr=mid;
}
if (rr==n)
{
indr=st;
/*
if ((indl==-1)&&(indr==-1))
{
indr=st;
}
ll=st-1;rr=n;
while (rr-ll>1)
{
mid=(ll+rr)/2;
a.clear();
a.push_back(mid);
a.push_back(indr);
tmp=doask(a);
if (tmp==0) ll=mid;
else rr=mid;
}
for (i=st;i<=ll;i++)
ans[i]=')';
for (i=rr;i<n;i++)
ans[i]='(';
st=n;
break;
*/
}
else
{
indr=rr;
/*
v=rr;
ans[v]=')';
ll=st-1;rr=v-1;
while (rr-ll>1)
{
mid=(ll+rr)/2;
a.clear();
a.push_back(mid);
a.push_back(v);
tmp=doask(a);
if (tmp==0) ll=mid;
else rr=mid;
}
for (i=st;i<=ll;i++)
ans[i]=')';
for (i=rr;i<v;i++)
ans[i]='(';
st=v+1;
indl=v-1;
indr=v;
*/
}
// }
//cerr<<"indr:"<<indr<<endl;
for (i=0;i<n;i+=8)
{
a.clear();
for (j=0;j<8;j++)
{
if (i+j>=n) break;
for (k=0;k<(1<<j);k++)
{
a.push_back(i+j);
a.push_back(indr);
a.push_back(indr);
}
}
tmp=doask(a);
for (j=0;j<8;j++)
{
if (i+j>=n) break;
if (((1<<j)&tmp)!=0) ans[i+j]='(';
else ans[i+j]=')';
}
}
doanswer();
}
}
solution (hard version)
#include <vector>
#include <list>
#include <map>
#include <set>
#include <deque>
#include <vector>
#include <list>
#include <map>
#include <set>
#include <deque>
#include <queue>
#include <stack>
#include <bitset>
#include <algorithm>
#include <functional>
#include <numeric>
#include <utility>
#include <sstream>
#include <iostream>
#include <iomanip>
#include <cstdio>
#include <cmath>
#include <cstdlib>
#include <cctype>
#include <string>
#include <cstring>
#include <ctime>
#include <random>
#include <chrono>
using namespace std;
#define _int64 long long
#define mo 998244353
char tar[110000];
int cnt;
int doask(vector<int>& a)
{
int i,ret;
/*
{
int v,tmp,j;
string ss;
cnt++;
ss="";
for (i=0;i<a.size();i++)
{
ss.push_back(tar[a[i]]);
}
tmp=0;
for (i=0;i<ss.length();i++)
{
v=0;
for (j=i;j<ss.length();j++)
{
if (ss[j]=='(') v++;
else v--;
if (v<0) break;
if (v==0) tmp++;
}
}
return tmp;
}
*/
printf("? %d",(int)a.size());
for (i=0;i<a.size();i++)
printf(" %d",a[i]+1);
printf("\n");
fflush(stdout);
scanf("%d",&ret);
return ret;
}
char ans[1100];
void doanswer()
{
printf("! %s\n",ans);
cerr<<"T "<<tar<<endl;
cerr<<"cnt:"<<cnt<<endl;
if (string(ans)==string(tar))
{
cerr<<"good"<<endl;
}
else cerr<<"bad"<<endl;
fflush(stdout);
}
int main()
{
int i,j,k,n,l,t,m,x,y,o,b1,indl,indr,ll,rr,mid,tmp,v,st,sum,len;
vector<int> a;
vector<int> p;
sum=0;len=0;
for (i=1;i<=130;i++)
{
tmp=i*(i+1)/2;
if (tmp>sum)
{
p.push_back(i);
//cerr<<"i:"<<i<<endl;
sum+=tmp;
len+=i+i;
}
}
//cerr<<"p.sz:"<<p.size()<<" "<<len<<endl;
scanf("%d",&t);
for (l=0;l<t;l++)
{
scanf("%d",&n);
for (i=0;i<n;i++)
{
x=rand()%2;
if (x==0) tar[i]='(';
else tar[i]=')';
}
tar[n]='\0';
cnt=0;
st=0;ans[n]='\0';
indl=-1;
indr=-1;
// while (st<n)
// {
ll=st;rr=n;
while (rr-ll>1)
{
mid=(ll+rr)/2;
a.clear();
for (i=st;i<=mid;i++)
a.push_back(i);
tmp=doask(a);
if (tmp==0) ll=mid;
else rr=mid;
}
if (rr==n)
{
indr=st;
/*
if ((indl==-1)&&(indr==-1))
{
indr=st;
}
ll=st-1;rr=n;
while (rr-ll>1)
{
mid=(ll+rr)/2;
a.clear();
a.push_back(mid);
a.push_back(indr);
tmp=doask(a);
if (tmp==0) ll=mid;
else rr=mid;
}
for (i=st;i<=ll;i++)
ans[i]=')';
for (i=rr;i<n;i++)
ans[i]='(';
st=n;
break;
*/
}
else
{
indr=rr;
/*
v=rr;
ans[v]=')';
ll=st-1;rr=v-1;
while (rr-ll>1)
{
mid=(ll+rr)/2;
a.clear();
a.push_back(mid);
a.push_back(v);
tmp=doask(a);
if (tmp==0) ll=mid;
else rr=mid;
}
for (i=st;i<=ll;i++)
ans[i]=')';
for (i=rr;i<v;i++)
ans[i]='(';
st=v+1;
indl=v-1;
indr=v;
*/
}
// }
//cerr<<"indr:"<<indr<<endl;
for (i=0;i<n;i+=p.size())
{
a.clear();
for (j=0;j<p.size();j++)
{
if (i+j>=n) break;
for (k=0;k<p[j];k++)
{
a.push_back(i+j);
a.push_back(indr);
}
a.push_back(indr);
}
tmp=doask(a);
for (j=p.size()-1;j>=0;j--)
{
if (i+j>=n) continue;
v=p[j]*(p[j]+1)/2;
if (tmp>=v)
{
tmp-=v;ans[i+j]='(';
}
else ans[i+j]=')';
}
}
doanswer();
}
}
D1D Permutation Blackhole
Hint1
Tutorial
solution
#include <bits/stdc++.h>
using namespace std;
constexpr int N = 5e5 + 10, V = 30, inf = INT_MAX, L = 18;
int n, q, S, a[N], b[N];
struct Rmq
{
int st[L][N];
void init()
{
for (int i = 1; i <= n; i++)
st[0][i] = a[i];
for (int i = 1; i < L; i++)
{
for (int j = 1; j <= n - (1 << i) + 1; j++)
{
st[i][j] = max(st[i - 1][j], st[i - 1][j + (1 << i - 1)]);
}
}
}
int qry(int l, int r)
{
int k = __lg(r - l + 1);
return max(st[k][l], st[k][r - (1 << k) + 1]);
}
} ds;
struct Info
{
int pre[V], suf[V];
int ans, l, r;
Info()
{
memset(pre, 0, sizeof(pre));
memset(suf, 0, sizeof(suf));
ans = inf, l = r = 0;
}
Info(const Info &x, const Info &y)
{
for (int i = 0; i < V; i++)
{
pre[i] = x.pre[i] ? x.pre[i] : y.pre[i];
suf[i] = y.suf[i] ? y.suf[i] : x.suf[i];
}
ans = inf, l = r = 0;
}
friend Info operator+(const Info &x, const Info &y)
{
Info z(x, y);
z.ans = min(x.ans, y.ans), z.l = x.l, z.r = y.r;
int pl = x.r, pr = y.l;
if (!pl)
return z;
for (int i = V - 1; i >= 0; i--)
{
int u = x.suf[i], v = y.pre[i];
if (u)
u = min(u, pl);
if (v)
v = max(v, pr);
int lans = u ? ds.qry(u, pr) : inf;
int rans = v ? ds.qry(pl, v) : inf;
if (lans < rans)
{
if (S >> i & 1)
{
if (lans < z.ans)
pl = u;
else
break;
}
else
z.ans = min(z.ans, lans);
}
else
{
if (S >> i & 1)
{
if (rans < z.ans)
pr = v;
else
break;
}
else
z.ans = min(z.ans, rans);
}
}
return z;
}
void clear()
{
*this = Info();
}
void upd(int p)
{
l = r = p, ans = b[p] > S ? a[p] : inf;
memset(pre, 0, sizeof(pre));
memset(suf, 0, sizeof(suf));
for (int j = 0; j < V; j++)
{
if (b[p] >> j & 1)
{
suf[j] = p;
pre[j] = p;
}
}
}
} ans;
void reply(const Info &cur)
{
ans = ans + cur;
}
struct Node
{
Node *ls, *rs;
int l, r, mid;
Info info;
void up()
{
info = ls->info + rs->info;
}
void update(int p)
{
if (l == r)
return info.upd(l);
p <= mid ? ls->update(p) : rs->update(p);
up();
}
void query(int ql, int qr)
{
if (l >= ql && r <= qr)
return reply(info);
if (ql <= mid)
ls->query(ql, qr);
if (qr > mid)
rs->query(ql, qr);
}
} mempool[N + 5 << 1], *cnt = mempool, *rt;
Node *build(int l, int r)
{
Node *u = cnt++;
int mid = l + r >> 1;
u->l = l, u->r = r, u->mid = mid;
if (l == r)
{
u->info.upd(l);
return u;
}
u->ls = build(l, mid);
u->rs = build(mid + 1, r);
u->up();
return u;
}
void solve()
{
cin >> n >> S;
S--;
for (int i = 1; i <= n; i++)
cin >> a[i];
ds.init();
for (int i = 1; i <= n; i++)
cin >> b[i];
rt = build(1, n);
cin >> q;
for (int i = 0; i < q; i++)
{
int opt;
cin >> opt;
if (opt == 1)
{
int p, x;
cin >> p >> x;
b[p] = x;
rt->update(p);
}
else
{
int l, r;
cin >> l >> r;
ans.clear();
rt->query(l, r);
cout << (ans.ans == inf ? -1 : ans.ans) << ' ';
}
}
cout << "\n";
}
int main()
{
#ifndef ONLINE_JUDGE
freopen("input.txt", "r", stdin);
freopen("output.txt", "w", stdout);
#endif
ios_base::sync_with_stdio(0);
cin.tie(0);
cout.tie(0);
ios::sync_with_stdio(false), cin.tie(0);
int t = 0;
cin >> t;
while (t--)
{
solve();
}
// cerr << 1.0 * clock() / CLOCKS_PER_SEC << '\n';
return 0;
}
D1E — Yet Yet Another Permutation Problem
Idea : wuhudsm
Hint1
We found it difficult to directly calculate all valid permutations. How about calculating all invalid permutations?
Hint2
Can you come up with an $$$O(n^2)$$$ DP solution?
Hint3
How to optimize the $$$O(n^2)$$$ DP solution?
Tutorial
We found it difficult to directly calculate all valid permutations. Consider calculating all invalid permutations and subtract it from $$$n!$$$.
Let's first make some notes. We note $$$P_i=max(p_1,...,p_i),Q_i=max(q_1,...,q_i)$$$,and $$$last_i=max(j)(P_j \neq P_i)$$$.And note the values of all extreme points is $$$val[1,...,m]$$$ and positions is $$$pos[1,...,m]$$$.
For example,$$$p=[3,1,4,2,6,5]$$$,we get $$$m=3$$$, $$$val=[3,4,6]$$$ and $$$pos=[1,3,5]$$$.
Let’s call a permutation $$$q$$$ “i-invalid” if there’s an index $$$j$$$ satisfying $$$Q_j=P_j=val_i$$$.
Note $$$S_i$$$ as the set of all “i-invalid” permutations.According to the Inclusion-Exclusion Principle,the answer is $$$n!-(|S_1|+|S_2|+...)+(|S_1∩S_2|+|S_1∩S_3|+...)-...$$$
Consider put $$$val_i$$$ in $$$q$$$ to make $$$q$$$ “i-invalid”.There’re $$$2$$$ kinds of ways:
$$$q_j=val_i(1 \le j \lt pos_i)$$$ and $$$Q_{pos_i}=val_i$$$;
$$$q_j=val_i(pos_i \le j \lt pos_{i+1})$$$ and $$$Q_{pos_i}=val_i$$$.
In both cases,we call the first $$$max(j,pos_i)$$$ numbers are “determined”. In other words,they’re some numbers in $$$[1,val_i]$$$.
Then we can find an $$$O(n^2)$$$ DP. Note $$$dp_{i,j}$$$ — $$$i$$$ stands for the first $$$i$$$ numbers is “determined”,and $$$j$$$ stands for $$$j$$$ extreme values are “invalid”. In fact we only care about the parity of $$$j$$$,so $$$j=0,1$$$. Note for the first $$$i$$$ numbers are “determined” automatically means $$$Q_i=P_i$$$. In addition,$$$q_i=P_i(P_i=P_{i-1})$$$ or $$$q_i \lt P_i(P_i \neq P_{i-1})$$$.
We get a DP formula $$$dp_{i,j}= \Sigma_{k=0}^{last_i} dp_{k,j \oplus 1} \cdot A^{P_i-k-1}_{i-k-1} \cdot (P_i==P_{i+1}?1:i-k)$$$
Here is a small trick to reduce the constant. If we note $$$f_i=dp_{i,0}-dp_{i,1}$$$, we get a more concise formula $$$f_{i}= -\Sigma_{j=0}^{last_i} f_{j} \cdot A^{P_i-j-1}_{i-j-1} \cdot (P_i==P_{i+1}?1:i-j)$$$
We perform the following transformation on the formula.
① For $$$i$$$ satisfying $$$P_i=P_{i-1}$$$, we get
$$$f_{i}= -\Sigma_{j=0}^{last_i} f_{j} \cdot A^{P_i-j-1}_{i-j-1} = -[(P_i-i)!]^{-1} \Sigma_{j=0}^{last_i} f_{j} \cdot (P_i-j-1)!$$$
② For $$$i$$$ satisfying $$$P_i \neq P_{i-1}$$$, we get
$$$f_{i}= -\Sigma_{j=0}^{last_i} f_{j} \cdot A^{P_i-j-1}_{i-j-1} \cdot (i-j) = -i[(P_i-i)!]^{-1} \Sigma_{j=0}^{last_i} f_{j} \cdot (P_i-j-1)! + [(P_i-i)!]^{-1} \Sigma_{j=0}^{last_i} jf_{j} \cdot (P_i-j-1)!$$$
After that, we have transformed the formula into a classical form $$$f(i)= h(i) \Sigma_{j=0}^{i-1} f(j) \cdot g(i-j)$$$, which can be calculate in $$$O(n log^2 n)$$$ by D&C+FFT.
solution
#include<bits/stdc++.h>
using namespace std;
#define all(a) a.begin(),a.end()
#define pb push_back
#define sz(a) ((int)a.size())
using ll=long long;
using u32=unsigned int;
using u64=unsigned long long;
using i128=__int128;
using u128=unsigned __int128;
using f128=__float128;
using pii=pair<int,int>;
using pll=pair<ll,ll>;
template<typename T> using vc=vector<T>;
template<typename T> using vvc=vc<vc<T>>;
template<typename T> using vvvc=vc<vvc<T>>;
using vi=vc<int>;
using vll=vc<ll>;
using vvi=vc<vi>;
using vvll=vc<vll>;
#define vv(type,name,n,...) \
vector<vector<type>> name(n,vector<type>(__VA_ARGS__))
#define vvv(type,name,n,m,...) \
vector<vector<vector<type>>> name(n,vector<vector<type>>(m,vector<type>(__VA_ARGS__)))
template<typename T> using min_heap=priority_queue<T,vector<T>,greater<T>>;
template<typename T> using max_heap=priority_queue<T>;
// https://trap.jp/post/1224/
#define rep1(n) for(ll i=0; i<(ll)(n); ++i)
#define rep2(i,n) for(ll i=0; i<(ll)(n); ++i)
#define rep3(i,a,b) for(ll i=(ll)(a); i<(ll)(b); ++i)
#define rep4(i,a,b,c) for(ll i=(ll)(a); i<(ll)(b); i+=(c))
#define cut4(a,b,c,d,e,...) e
#define rep(...) cut4(__VA_ARGS__,rep4,rep3,rep2,rep1)(__VA_ARGS__)
#define per1(n) for(ll i=((ll)n)-1; i>=0; --i)
#define per2(i,n) for(ll i=((ll)n)-1; i>=0; --i)
#define per3(i,a,b) for(ll i=((ll)a)-1; i>=(ll)(b); --i)
#define per4(i,a,b,c) for(ll i=((ll)a)-1; i>=(ll)(b); i-=(c))
#define per(...) cut4(__VA_ARGS__,per4,per3,per2,per1)(__VA_ARGS__)
#define rep_subset(i,s) for(ll i=(s); i>=0; i=(i==0?-1:(i-1)&(s)))
template<typename T, typename S> constexpr T ifloor(const T a, const S b){return a/b-(a%b&&(a^b)<0);}
template<typename T, typename S> constexpr T iceil(const T a, const S b){return ifloor(a+b-1,b);}
template<typename T>
void sort_unique(vector<T> &vec){
sort(vec.begin(),vec.end());
vec.resize(unique(vec.begin(),vec.end())-vec.begin());
}
template<typename T, typename S> constexpr bool chmin(T &a, const S b){if(a>b) return a=b,true; return false;}
template<typename T, typename S> constexpr bool chmax(T &a, const S b){if(a<b) return a=b,true; return false;}
template<typename T, typename S> istream& operator >> (istream& i, pair<T,S> &p){return i >> p.first >> p.second;}
template<typename T, typename S> ostream& operator << (ostream& o, const pair<T,S> &p){return o << p.first << ' ' << p.second;}
#ifdef i_am_noob
#define bug(...) cerr << "#" << __LINE__ << ' ' << #__VA_ARGS__ << "- ", _do(__VA_ARGS__)
template<typename T> void _do(vector<T> x){for(auto i: x) cerr << i << ' ';cerr << "\n";}
template<typename T> void _do(set<T> x){for(auto i: x) cerr << i << ' ';cerr << "\n";}
template<typename T> void _do(unordered_set<T> x){for(auto i: x) cerr << i << ' ';cerr << "\n";}
template<typename T> void _do(T && x) {cerr << x << endl;}
template<typename T, typename ...S> void _do(T && x, S&&...y) {cerr << x << ", "; _do(y...);}
#else
#define bug(...) 777771449
#endif
template<typename T> void print(vector<T> x){for(auto i: x) cout << i << ' ';cout << "\n";}
template<typename T> void print(set<T> x){for(auto i: x) cout << i << ' ';cout << "\n";}
template<typename T> void print(unordered_set<T> x){for(auto i: x) cout << i << ' ';cout << "\n";}
template<typename T> void print(T && x) {cout << x << "\n";}
template<typename T, typename... S> void print(T && x, S&&... y) {cout << x << ' ';print(y...);}
template<typename T> istream& operator >> (istream& i, vector<T> &vec){for(auto &x: vec) i >> x; return i;}
vvi read_graph(int n, int m, int base=1){
vvi adj(n);
for(int i=0,u,v; i<m; ++i){
cin >> u >> v,u-=base,v-=base;
adj[u].pb(v),adj[v].pb(u);
}
return adj;
}
vvi read_tree(int n, int base=1){return read_graph(n,n-1,base);}
template<typename T, typename S> pair<T,S> operator + (const pair<T,S> &a, const pair<T,S> &b){return {a.first+b.first,a.second+b.second};}
template<typename T> constexpr T inf=0;
template<> constexpr int inf<int> = 0x3f3f3f3f;
template<> constexpr ll inf<ll> = 0x3f3f3f3f3f3f3f3f;
template<typename T> vector<T> operator += (vector<T> &a, int val){for(auto &i: a) i+=val; return a;}
template<typename T> T isqrt(const T &x){T y=sqrt(x+2); while(y*y>x) y--; return y;}
#define ykh mt19937 rng(chrono::steady_clock::now().time_since_epoch().count())
#include <utility>
namespace atcoder {
namespace internal {
// @param m `1 <= m`
// @return x mod m
constexpr long long safe_mod(long long x, long long m) {
x %= m;
if (x < 0) x += m;
return x;
}
// Fast modular multiplication by barrett reduction
// Reference: https://en.wikipedia.org/wiki/Barrett_reduction
// NOTE: reconsider after Ice Lake
struct barrett {
unsigned int _m;
unsigned long long im;
// @param m `1 <= m < 2^31`
barrett(unsigned int m) : _m(m), im((unsigned long long)(-1) / m + 1) {}
// @return m
unsigned int umod() const { return _m; }
// @param a `0 <= a < m`
// @param b `0 <= b < m`
// @return `a * b % m`
unsigned int mul(unsigned int a, unsigned int b) const {
// [1] m = 1
// a = b = im = 0, so okay
// [2] m >= 2
// im = ceil(2^64 / m)
// -> im * m = 2^64 + r (0 <= r < m)
// let z = a*b = c*m + d (0 <= c, d < m)
// a*b * im = (c*m + d) * im = c*(im*m) + d*im = c*2^64 + c*r + d*im
// c*r + d*im < m * m + m * im < m * m + 2^64 + m <= 2^64 + m * (m + 1) < 2^64 * 2
// ((ab * im) >> 64) == c or c + 1
unsigned long long z = a;
z *= b;
#ifdef _MSC_VER
unsigned long long x;
_umul128(z, im, &x);
#else
unsigned long long x =
(unsigned long long)(((unsigned __int128)(z)*im) >> 64);
#endif
unsigned int v = (unsigned int)(z - x * _m);
if (_m <= v) v += _m;
return v;
}
};
// @param n `0 <= n`
// @param m `1 <= m`
// @return `(x ** n) % m`
constexpr long long pow_mod_constexpr(long long x, long long n, int m) {
if (m == 1) return 0;
unsigned int _m = (unsigned int)(m);
unsigned long long r = 1;
unsigned long long y = safe_mod(x, m);
while (n) {
if (n & 1) r = (r * y) % _m;
y = (y * y) % _m;
n >>= 1;
}
return r;
}
// Reference:
// M. Forisek and J. Jancina,
// Fast Primality Testing for Integers That Fit into a Machine Word
// @param n `0 <= n`
constexpr bool is_prime_constexpr(int n) {
if (n <= 1) return false;
if (n == 2 || n == 7 || n == 61) return true;
if (n % 2 == 0) return false;
long long d = n - 1;
while (d % 2 == 0) d /= 2;
constexpr long long bases[3] = {2, 7, 61};
for (long long a : bases) {
long long t = d;
long long y = pow_mod_constexpr(a, t, n);
while (t != n - 1 && y != 1 && y != n - 1) {
y = y * y % n;
t <<= 1;
}
if (y != n - 1 && t % 2 == 0) {
return false;
}
}
return true;
}
template <int n> constexpr bool is_prime = is_prime_constexpr(n);
// @param b `1 <= b`
// @return pair(g, x) s.t. g = gcd(a, b), xa = g (mod b), 0 <= x < b/g
constexpr std::pair<long long, long long> inv_gcd(long long a, long long b) {
a = safe_mod(a, b);
if (a == 0) return {b, 0};
// Contracts:
// [1] s - m0 * a = 0 (mod b)
// [2] t - m1 * a = 0 (mod b)
// [3] s * |m1| + t * |m0| <= b
long long s = b, t = a;
long long m0 = 0, m1 = 1;
while (t) {
long long u = s / t;
s -= t * u;
m0 -= m1 * u; // |m1 * u| <= |m1| * s <= b
// [3]:
// (s - t * u) * |m1| + t * |m0 - m1 * u|
// <= s * |m1| - t * u * |m1| + t * (|m0| + |m1| * u)
// = s * |m1| + t * |m0| <= b
auto tmp = s;
s = t;
t = tmp;
tmp = m0;
m0 = m1;
m1 = tmp;
}
// by [3]: |m0| <= b/g
// by g != b: |m0| < b/g
if (m0 < 0) m0 += b / s;
return {s, m0};
}
// Compile time primitive root
// @param m must be prime
// @return primitive root (and minimum in now)
constexpr int primitive_root_constexpr(int m) {
if (m == 2) return 1;
if (m == 167772161) return 3;
if (m == 469762049) return 3;
if (m == 754974721) return 11;
if (m == 998244353) return 3;
int divs[20] = {};
divs[0] = 2;
int cnt = 1;
int x = (m - 1) / 2;
while (x % 2 == 0) x /= 2;
for (int i = 3; (long long)(i)*i <= x; i += 2) {
if (x % i == 0) {
divs[cnt++] = i;
while (x % i == 0) {
x /= i;
}
}
}
if (x > 1) {
divs[cnt++] = x;
}
for (int g = 2;; g++) {
bool ok = true;
for (int i = 0; i < cnt; i++) {
if (pow_mod_constexpr(g, (m - 1) / divs[i], m) == 1) {
ok = false;
break;
}
}
if (ok) return g;
}
}
template <int m> constexpr int primitive_root = primitive_root_constexpr(m);
} // namespace internal
} // namespace atcoder
#include <cassert>
#include <numeric>
#include <type_traits>
namespace atcoder {
namespace internal {
#ifndef _MSC_VER
template <class T>
using is_signed_int128 =
typename std::conditional<std::is_same<T, __int128_t>::value ||
std::is_same<T, __int128>::value,
std::true_type,
std::false_type>::type;
template <class T>
using is_unsigned_int128 =
typename std::conditional<std::is_same<T, __uint128_t>::value ||
std::is_same<T, unsigned __int128>::value,
std::true_type,
std::false_type>::type;
template <class T>
using make_unsigned_int128 =
typename std::conditional<std::is_same<T, __int128_t>::value,
__uint128_t,
unsigned __int128>;
template <class T>
using is_integral = typename std::conditional<std::is_integral<T>::value ||
is_signed_int128<T>::value ||
is_unsigned_int128<T>::value,
std::true_type,
std::false_type>::type;
template <class T>
using is_signed_int = typename std::conditional<(is_integral<T>::value &&
std::is_signed<T>::value) ||
is_signed_int128<T>::value,
std::true_type,
std::false_type>::type;
template <class T>
using is_unsigned_int =
typename std::conditional<(is_integral<T>::value &&
std::is_unsigned<T>::value) ||
is_unsigned_int128<T>::value,
std::true_type,
std::false_type>::type;
template <class T>
using to_unsigned = typename std::conditional<
is_signed_int128<T>::value,
make_unsigned_int128<T>,
typename std::conditional<std::is_signed<T>::value,
std::make_unsigned<T>,
std::common_type<T>>::type>::type;
#else
template <class T> using is_integral = typename std::is_integral<T>;
template <class T>
using is_signed_int =
typename std::conditional<is_integral<T>::value && std::is_signed<T>::value,
std::true_type,
std::false_type>::type;
template <class T>
using is_unsigned_int =
typename std::conditional<is_integral<T>::value &&
std::is_unsigned<T>::value,
std::true_type,
std::false_type>::type;
template <class T>
using to_unsigned = typename std::conditional<is_signed_int<T>::value,
std::make_unsigned<T>,
std::common_type<T>>::type;
#endif
template <class T>
using is_signed_int_t = std::enable_if_t<is_signed_int<T>::value>;
template <class T>
using is_unsigned_int_t = std::enable_if_t<is_unsigned_int<T>::value>;
template <class T> using to_unsigned_t = typename to_unsigned<T>::type;
} // namespace internal
} // namespace atcoder
#include <cassert>
#include <numeric>
#include <type_traits>
#ifdef _MSC_VER
#include <intrin.h>
#endif
namespace atcoder {
namespace internal {
struct modint_base {};
struct static_modint_base : modint_base {};
template <class T> using is_modint = std::is_base_of<modint_base, T>;
template <class T> using is_modint_t = std::enable_if_t<is_modint<T>::value>;
} // namespace internal
template <int m, std::enable_if_t<(1 <= m)>* = nullptr>
struct static_modint : internal::static_modint_base {
using mint = static_modint;
public:
static constexpr int mod() { return m; }
static mint raw(int v) {
mint x;
x._v = v;
return x;
}
static_modint() : _v(0) {}
template <class T, internal::is_signed_int_t<T>* = nullptr>
static_modint(T v) {
long long x = (long long)(v % (long long)(umod()));
if (x < 0) x += umod();
_v = (unsigned int)(x);
}
template <class T, internal::is_unsigned_int_t<T>* = nullptr>
static_modint(T v) {
_v = (unsigned int)(v % umod());
}
static_modint(bool v) { _v = ((unsigned int)(v) % umod()); }
unsigned int val() const { return _v; }
mint& operator++() {
_v++;
if (_v == umod()) _v = 0;
return *this;
}
mint& operator--() {
if (_v == 0) _v = umod();
_v--;
return *this;
}
mint operator++(int) {
mint result = *this;
++*this;
return result;
}
mint operator--(int) {
mint result = *this;
--*this;
return result;
}
mint& operator+=(const mint& rhs) {
_v += rhs._v;
if (_v >= umod()) _v -= umod();
return *this;
}
mint& operator-=(const mint& rhs) {
_v -= rhs._v;
if (_v >= umod()) _v += umod();
return *this;
}
mint& operator*=(const mint& rhs) {
unsigned long long z = _v;
z *= rhs._v;
_v = (unsigned int)(z % umod());
return *this;
}
mint& operator/=(const mint& rhs) { return *this = *this * rhs.inv(); }
mint operator+() const { return *this; }
mint operator-() const { return mint() - *this; }
mint pow(long long n) const {
assert(0 <= n);
mint x = *this, r = 1;
while (n) {
if (n & 1) r *= x;
x *= x;
n >>= 1;
}
return r;
}
mint inv() const {
if (prime) {
assert(_v);
return pow(umod() - 2);
} else {
auto eg = internal::inv_gcd(_v, m);
assert(eg.first == 1);
return eg.second;
}
}
friend mint operator+(const mint& lhs, const mint& rhs) {
return mint(lhs) += rhs;
}
friend mint operator-(const mint& lhs, const mint& rhs) {
return mint(lhs) -= rhs;
}
friend mint operator*(const mint& lhs, const mint& rhs) {
return mint(lhs) *= rhs;
}
friend mint operator/(const mint& lhs, const mint& rhs) {
return mint(lhs) /= rhs;
}
friend bool operator==(const mint& lhs, const mint& rhs) {
return lhs._v == rhs._v;
}
friend bool operator!=(const mint& lhs, const mint& rhs) {
return lhs._v != rhs._v;
}
private:
unsigned int _v;
static constexpr unsigned int umod() { return m; }
static constexpr bool prime = internal::is_prime<m>;
};
template <int id> struct dynamic_modint : internal::modint_base {
using mint = dynamic_modint;
public:
static int mod() { return (int)(bt.umod()); }
static void set_mod(int m) {
assert(1 <= m);
bt = internal::barrett(m);
}
static mint raw(int v) {
mint x;
x._v = v;
return x;
}
dynamic_modint() : _v(0) {}
template <class T, internal::is_signed_int_t<T>* = nullptr>
dynamic_modint(T v) {
long long x = (long long)(v % (long long)(mod()));
if (x < 0) x += mod();
_v = (unsigned int)(x);
}
template <class T, internal::is_unsigned_int_t<T>* = nullptr>
dynamic_modint(T v) {
_v = (unsigned int)(v % mod());
}
dynamic_modint(bool v) { _v = ((unsigned int)(v) % mod()); }
unsigned int val() const { return _v; }
mint& operator++() {
_v++;
if (_v == umod()) _v = 0;
return *this;
}
mint& operator--() {
if (_v == 0) _v = umod();
_v--;
return *this;
}
mint operator++(int) {
mint result = *this;
++*this;
return result;
}
mint operator--(int) {
mint result = *this;
--*this;
return result;
}
mint& operator+=(const mint& rhs) {
_v += rhs._v;
if (_v >= umod()) _v -= umod();
return *this;
}
mint& operator-=(const mint& rhs) {
_v += mod() - rhs._v;
if (_v >= umod()) _v -= umod();
return *this;
}
mint& operator*=(const mint& rhs) {
_v = bt.mul(_v, rhs._v);
return *this;
}
mint& operator/=(const mint& rhs) { return *this = *this * rhs.inv(); }
mint operator+() const { return *this; }
mint operator-() const { return mint() - *this; }
mint pow(long long n) const {
assert(0 <= n);
mint x = *this, r = 1;
while (n) {
if (n & 1) r *= x;
x *= x;
n >>= 1;
}
return r;
}
mint inv() const {
auto eg = internal::inv_gcd(_v, mod());
assert(eg.first == 1);
return eg.second;
}
friend mint operator+(const mint& lhs, const mint& rhs) {
return mint(lhs) += rhs;
}
friend mint operator-(const mint& lhs, const mint& rhs) {
return mint(lhs) -= rhs;
}
friend mint operator*(const mint& lhs, const mint& rhs) {
return mint(lhs) *= rhs;
}
friend mint operator/(const mint& lhs, const mint& rhs) {
return mint(lhs) /= rhs;
}
friend bool operator==(const mint& lhs, const mint& rhs) {
return lhs._v == rhs._v;
}
friend bool operator!=(const mint& lhs, const mint& rhs) {
return lhs._v != rhs._v;
}
private:
unsigned int _v;
static internal::barrett bt;
static unsigned int umod() { return bt.umod(); }
};
template <int id> internal::barrett dynamic_modint<id>::bt = 998244353;
using modint998244353 = static_modint<998244353>;
using modint1000000007 = static_modint<1000000007>;
using modint = dynamic_modint<-1>;
namespace internal {
template <class T>
using is_static_modint = std::is_base_of<internal::static_modint_base, T>;
template <class T>
using is_static_modint_t = std::enable_if_t<is_static_modint<T>::value>;
template <class> struct is_dynamic_modint : public std::false_type {};
template <int id>
struct is_dynamic_modint<dynamic_modint<id>> : public std::true_type {};
template <class T>
using is_dynamic_modint_t = std::enable_if_t<is_dynamic_modint<T>::value>;
} // namespace internal
} // namespace atcoder
#include <algorithm>
#include <array>
#ifdef _MSC_VER
#include <intrin.h>
#endif
namespace atcoder {
namespace internal {
// @param n `0 <= n`
// @return minimum non-negative `x` s.t. `n <= 2**x`
int ceil_pow2(int n) {
int x = 0;
while ((1U << x) < (unsigned int)(n)) x++;
return x;
}
// @param n `1 <= n`
// @return minimum non-negative `x` s.t. `(n & (1 << x)) != 0`
int bsf(unsigned int n) {
#ifdef _MSC_VER
unsigned long index;
_BitScanForward(&index, n);
return index;
#else
return __builtin_ctz(n);
#endif
}
} // namespace internal
} // namespace atcoder
#include <cassert>
#include <type_traits>
#include <vector>
namespace atcoder {
namespace internal {
template <class mint, internal::is_static_modint_t<mint>* = nullptr>
void butterfly(std::vector<mint>& a) {
static constexpr int g = internal::primitive_root<mint::mod()>;
int n = int(a.size());
int h = internal::ceil_pow2(n);
static bool first = true;
static mint sum_e[30]; // sum_e[i] = ies[0] * ... * ies[i - 1] * es[i]
if (first) {
first = false;
mint es[30], ies[30]; // es[i]^(2^(2+i)) == 1
int cnt2 = bsf(mint::mod() - 1);
mint e = mint(g).pow((mint::mod() - 1) >> cnt2), ie = e.inv();
for (int i = cnt2; i >= 2; i--) {
// e^(2^i) == 1
es[i - 2] = e;
ies[i - 2] = ie;
e *= e;
ie *= ie;
}
mint now = 1;
for (int i = 0; i <= cnt2 - 2; i++) {
sum_e[i] = es[i] * now;
now *= ies[i];
}
}
for (int ph = 1; ph <= h; ph++) {
int w = 1 << (ph - 1), p = 1 << (h - ph);
mint now = 1;
for (int s = 0; s < w; s++) {
int offset = s << (h - ph + 1);
for (int i = 0; i < p; i++) {
auto l = a[i + offset];
auto r = a[i + offset + p] * now;
a[i + offset] = l + r;
a[i + offset + p] = l - r;
}
now *= sum_e[bsf(~(unsigned int)(s))];
}
}
}
template <class mint, internal::is_static_modint_t<mint>* = nullptr>
void butterfly_inv(std::vector<mint>& a) {
static constexpr int g = internal::primitive_root<mint::mod()>;
int n = int(a.size());
int h = internal::ceil_pow2(n);
static bool first = true;
static mint sum_ie[30]; // sum_ie[i] = es[0] * ... * es[i - 1] * ies[i]
if (first) {
first = false;
mint es[30], ies[30]; // es[i]^(2^(2+i)) == 1
int cnt2 = bsf(mint::mod() - 1);
mint e = mint(g).pow((mint::mod() - 1) >> cnt2), ie = e.inv();
for (int i = cnt2; i >= 2; i--) {
// e^(2^i) == 1
es[i - 2] = e;
ies[i - 2] = ie;
e *= e;
ie *= ie;
}
mint now = 1;
for (int i = 0; i <= cnt2 - 2; i++) {
sum_ie[i] = ies[i] * now;
now *= es[i];
}
}
for (int ph = h; ph >= 1; ph--) {
int w = 1 << (ph - 1), p = 1 << (h - ph);
mint inow = 1;
for (int s = 0; s < w; s++) {
int offset = s << (h - ph + 1);
for (int i = 0; i < p; i++) {
auto l = a[i + offset];
auto r = a[i + offset + p];
a[i + offset] = l + r;
a[i + offset + p] =
(unsigned long long)(mint::mod() + l.val() - r.val()) *
inow.val();
}
inow *= sum_ie[bsf(~(unsigned int)(s))];
}
}
}
} // namespace internal
template <class mint, internal::is_static_modint_t<mint>* = nullptr>
std::vector<mint> convolution(std::vector<mint> a, std::vector<mint> b) {
int n = int(a.size()), m = int(b.size());
if (!n || !m) return {};
if (std::min(n, m) <= 60) {
if (n < m) {
std::swap(n, m);
std::swap(a, b);
}
std::vector<mint> ans(n + m - 1);
for (int i = 0; i < n; i++) {
for (int j = 0; j < m; j++) {
ans[i + j] += a[i] * b[j];
}
}
return ans;
}
int z = 1 << internal::ceil_pow2(n + m - 1);
a.resize(z);
internal::butterfly(a);
b.resize(z);
internal::butterfly(b);
for (int i = 0; i < z; i++) {
a[i] *= b[i];
}
internal::butterfly_inv(a);
a.resize(n + m - 1);
mint iz = mint(z).inv();
for (int i = 0; i < n + m - 1; i++) a[i] *= iz;
return a;
}
template <unsigned int mod = 998244353,
class T,
std::enable_if_t<internal::is_integral<T>::value>* = nullptr>
std::vector<T> convolution(const std::vector<T>& a, const std::vector<T>& b) {
int n = int(a.size()), m = int(b.size());
if (!n || !m) return {};
using mint = static_modint<mod>;
std::vector<mint> a2(n), b2(m);
for (int i = 0; i < n; i++) {
a2[i] = mint(a[i]);
}
for (int i = 0; i < m; i++) {
b2[i] = mint(b[i]);
}
auto c2 = convolution(move(a2), move(b2));
std::vector<T> c(n + m - 1);
for (int i = 0; i < n + m - 1; i++) {
c[i] = c2[i].val();
}
return c;
}
std::vector<long long> convolution_ll(const std::vector<long long>& a,
const std::vector<long long>& b) {
int n = int(a.size()), m = int(b.size());
if (!n || !m) return {};
static constexpr unsigned long long MOD1 = 754974721; // 2^24
static constexpr unsigned long long MOD2 = 167772161; // 2^25
static constexpr unsigned long long MOD3 = 469762049; // 2^26
static constexpr unsigned long long M2M3 = MOD2 * MOD3;
static constexpr unsigned long long M1M3 = MOD1 * MOD3;
static constexpr unsigned long long M1M2 = MOD1 * MOD2;
static constexpr unsigned long long M1M2M3 = MOD1 * MOD2 * MOD3;
static constexpr unsigned long long i1 =
internal::inv_gcd(MOD2 * MOD3, MOD1).second;
static constexpr unsigned long long i2 =
internal::inv_gcd(MOD1 * MOD3, MOD2).second;
static constexpr unsigned long long i3 =
internal::inv_gcd(MOD1 * MOD2, MOD3).second;
auto c1 = convolution<MOD1>(a, b);
auto c2 = convolution<MOD2>(a, b);
auto c3 = convolution<MOD3>(a, b);
std::vector<long long> c(n + m - 1);
for (int i = 0; i < n + m - 1; i++) {
unsigned long long x = 0;
x += (c1[i] * i1) % MOD1 * M2M3;
x += (c2[i] * i2) % MOD2 * M1M3;
x += (c3[i] * i3) % MOD3 * M1M2;
// B = 2^63, -B <= x, r(real value) < B
// (x, x - M, x - 2M, or x - 3M) = r (mod 2B)
// r = c1[i] (mod MOD1)
// focus on MOD1
// r = x, x - M', x - 2M', x - 3M' (M' = M % 2^64) (mod 2B)
// r = x,
// x - M' + (0 or 2B),
// x - 2M' + (0, 2B or 4B),
// x - 3M' + (0, 2B, 4B or 6B) (without mod!)
// (r - x) = 0, (0)
// - M' + (0 or 2B), (1)
// -2M' + (0 or 2B or 4B), (2)
// -3M' + (0 or 2B or 4B or 6B) (3) (mod MOD1)
// we checked that
// ((1) mod MOD1) mod 5 = 2
// ((2) mod MOD1) mod 5 = 3
// ((3) mod MOD1) mod 5 = 4
long long diff =
c1[i] - internal::safe_mod((long long)(x), (long long)(MOD1));
if (diff < 0) diff += MOD1;
static constexpr unsigned long long offset[5] = {
0, 0, M1M2M3, 2 * M1M2M3, 3 * M1M2M3};
x -= offset[diff % 5];
c[i] = x;
}
return c;
}
} // namespace atcoder
using namespace atcoder;
using mint=modint998244353;
//using mint=modint1000000007;
template<int mod>
struct nCr{
vector<static_modint<mod>> fac,inv,ifac;
void calc(int n){
fac.resize(n+1),inv.resize(n+1),ifac.resize(n+1);
fac[0]=inv[1]=ifac[0]=1;
for(int i=1; i<=n; ++i) fac[i]=fac[i-1]*static_modint<mod>::raw(i);
for(int i=2; i<=n; ++i) inv[i]=inv[mod%i]*static_modint<mod>::raw(mod-mod/i);
for(int i=1; i<=n; ++i) ifac[i]=ifac[i-1]*inv[i];
}
static_modint<mod> C(int n, int m){
if(m<0||m>n) return 0;
return fac[n]*ifac[m]*ifac[n-m];
}
};
nCr<998244353> de;
void ahcorz(){
// (a_i - j)! / (a_i - i)! (a_i = a_j)
// (a_i - j - 1)! / (a_i - i)! * (i - j) (a_i != a_j)
int n; cin >> n;
de.calc(n+1);
vi a; a.pb(0);
rep(n){
int x; cin >> x;
a.pb(max(a.back(),x));
}
vc<mint> dp(n+1),aux(n+1);
auto solve=[&](auto &self, int l, int r) -> void{
if(l+1==r){
if(l==0) dp[l]=-1;
else dp[l]+=aux[l]*l;
dp[l]*=de.ifac[a[l]-l];
return;
}
int mid=l+r>>1;
self(self,l,mid);
if(a[mid]==a[mid-1]){
int x=mid-1,y=mid;
while(x-1>=l&&a[x-1]==a[mid]) x--;
while(y+1<r&&a[y+1]==a[mid]) y++;
mint tot,tot1,tot2;
rep(i,x,mid) tot+=dp[i]*de.fac[a[mid]-i],tot1+=dp[i]*de.fac[a[mid]-i-1],tot2+=dp[i]*de.fac[a[mid]-i-1]*i;
rep(i,mid,y+1) dp[i]-=tot-(tot1*i-tot2);
}
vc<mint> vec;
rep(i,a[mid]-mid,a[r-1]-(l+1)+1) vec.pb(de.fac[i]);
vc<mint> vec1,vec2;
rep(i,l,mid) vec1.pb(dp[i]),vec2.pb(dp[i]*i);
auto res1=convolution(vec,vec1),res2=convolution(vec,vec2);
int base=a[mid]-mid+l+1;
rep(i,mid,r) if(a[i]-base>=0&&a[i]-base<sz(res1)) aux[i]-=res1[a[i]-base],dp[i]+=res2[a[i]-base];
self(self,mid,r);
};
solve(solve,0,n+1);
print(dp[n].val());
}
signed main(){
ios_base::sync_with_stdio(0),cin.tie(0);
cout << fixed << setprecision(20);
int t=1;
cin >> t;
while(t--) ahcorz();
}
D1F — Grand Finale: Circles
Idea : chromate00
UPD: The TL was a bit too loose, and unintended solutions with $$$\mathcal{O}(n \log^2 (1/\epsilon))$$$ complexity passed during the round. For the sake of the problem itself, I suggest that you try to solve assuming a $$$1$$$ second time limit.
Tutorial$$$f(x,y)=\min_{1 \le i \le n} {\left({r_i-\sqrt{(x_i-x)^2+(y_i-y)^2}}\right)}$$$ $$$\begin{split} \max & \ r \\ \text{s.t.} & \ \sqrt{(x_i-x)^2+(y_i-y)^2} \le r_i-r \ \ \forall \ 1 \le i \le n \end{split}$$$ $$$(x_A-x)^2+(y_A-y)^2 = (r_A-r)^2 \\ (x_B-x)^2+(y_B-y)^2 = (r_B-r)^2 \\ (x_C-x)^2+(y_C-y)^2 = (r_C-r)^2$$$
For a given center coordinate $$$(x,y)$$$, we can model the objective function $$$r=f(x,y)$$$ to maximize as follows.
Formally, this can be modeled as follows:
This is a second order cone program. Generally, this is not as simple to solve as linear or quadratic programs. There are multiple solvers out there, but all of them are either commercial or simply unfit for use on Codeforces. And even if we use general approaches to second order cone programming, their time complexity is at least $$$\mathcal{O}(n^{3.5})$$$, so they cannot work. Almost all general approaches to convex programming do not help very much for this task also, because of the following reasons.
Gradient Descent cannot work in this task. The objective function is not smooth.
Coordinate Descent also cannot work in this task for the same reason, and it is not hard to find a countercase.
Subgradient methods can solve this task, but they take $$$\mathcal{O}(1/\epsilon^2)$$$ iterations to converge. This is a very big dependency on the precision, so they take too much time.
We do not expect anyone to write Adaptive Coordinate Descent during contest time, but we still tested the PRAXIS algorithm and checked that it does not pass.
Using nested ternary searches does not help very much. The time complexity will be $$$\mathcal{O}(n \log^2 (1/\epsilon))$$$, but the constants are large due to multiple function calls and floating point operations.
Heuristics such as the Nelder-Mead method are not proven to converge. In many cases we found it not converging successfully.
In the editorial, we will explain a solution of expected $$$\mathcal{O}(n)$$$ time complexity, with constants depending on the number of dimensions (which is $$$2$$$ in our case).
First, we observe that any answer can be described as a "basis" of at most $$$3$$$ selected circles. For $$$k=0,1,2,3$$$, the basis is as follows.
- Case 1: $$$k=0$$$
This is the case of no circles, which can be regarded as an "identity element" of bases. Some very huge circle that will enclose any given circle, for example the circle with $$$10^{18}$$$ as radius and $$$(0,0)$$$ as center, can represent this case.
- Case 2: $$$k=1$$$
This is the case of a single circle. The circle that represents this case is the one selected circle.
- Case 3: $$$k=2$$$
This is the case of two circles, for which the basis is the largest circle in the intersection of two circles. We can find this by choosing the center as the midpoint of the two points $$$A$$$ and $$$B$$$, where $$$A$$$ and $$$B$$$ are the points where the line segment between circle centers and the circle itself intersect.
- Case 4: $$$k=3$$$
This is the case of three circles, and the most tricky of all four cases. It can be found as one solution to the following system of quadratic equations:
This can be solved algebraically. If you are used to solving geometry problems in Olympiads, you may know the ''Problem of Apollonius''. This is one special case of the problem, which is not very often mentioned because often the three circles do not intersect with each other (thus most of the time one solution for $$$r$$$ must be negative). Or if one does not know the problem, they can find that the locus of the circle center that meets with two circle edges (as in Case 3) is a conic section, precisely one side of a pair of hyperbola. By finding the intersections of two hyperbola, one can find the center of this circle. Thus, this circle can be found in $$$\mathcal{O}(1)$$$ time, with constants based on number of dimensions.
By the system of quadratic equations stated above, we see that there are at most $$$2$$$ solutions. It is sufficient to take the one that happens to be inside the intersection. (If there are two such solutions, take the larger one. For the current version of the data this check is not necessary, but it does not hurt to try.) Thus, for some $$$n \gt 3$$$, the bases can be represented with some $$$k \le 3$$$ circles.
Do note that, already, we have an $$$\mathcal{O}(n^3)$$$ solution based on this. We enumerate all possible bases, and find the smallest one out of them (since that one will be the one that satisfies the conditions).
Now, for some basis circle $$$(x,y,r)$$$, we can check whether some given circle $$$(x_i,y_i,r_i)$$$ is "violated" by this basis (i.e. this basis is not inside this circle) in $$$\mathcal{O}(1)$$$ by comparing the distance. We take an incremental approach. We add circles to the basis one by one, up to $$$3$$$ circles. We choose to change the basis if only the next added circle is violated by the current basis. Of course, still this takes $$$\mathcal{O}(n^3)$$$ time, and sorting the circles based on any argument will not help anyways.
However, adding the circles in random order will make it only take expected $$$\mathcal{O}(n)$$$ time. This can be analyzed by the probability that a new basis will be added. The precise analysis for the time complexity is left as a practice for the reader.
If you have seen the minimum enclosing circle problem, the Welzl's algorithm to solve that precise problem will seem very similar to this. In fact, the suggested algorithm itself is not very different from Welzl's algorithm. If you want, though, you may take an iterative approach instead of modifying Welzl's algorithm. Using three nested loops will do.
There is only two issues left. The first is the issue of numerical instability. Sometimes, due to numerical instability, the algorithm will return a NaN value. It is not that hard to solve, we can simply shuffle again and rerun the algorithm. The given time limit will be sufficient for doing this unless the constant is too big. The second is more tricky, and it is about hacks.
If you use a fixed seed for randomization (or time-based seeds to some extent), someone may hack you by adding an adversarial test case which makes the time complexity explode back to $$$\mathcal{O}(n^3)$$$. This is very bad. How do we solve this? To intervent this, we will make the algorithm halt when the number of iterations exceed a certain limit. Formally, let us choose a constant $$$c$$$, and the algorithm will halt and rerun when the number of iterations exceeds $$$c \cdot n$$$. Then, if the probability of the algorithm successfully running in this number of iterations is $$$p$$$, the expected number of reruns needed is $$$\mathcal{O}(1/p)$$$, and hacking is much harder. The algorithm will successfully find the solution in $$$\mathcal{O}((c/p) \cdot n)$$$ time. Of course, the probability $$$p$$$ depends on $$$c$$$, and you must tune the value as needed. Empirically $$$c=30$$$ runs very well, almost always taking no longer than $$$0.4$$$ seconds under proper optimizations, and $$$1.5$$$ seconds without.
solution
#include<bits/stdc++.h>
using namespace std;
using ll=long long;
using lf=long double;
using pt=pair<lf,lf>;
lf& real(pt& p){return p.first;}
lf& imag(pt& p){return p.second;}
pt midp(pt a,pt b){return pt{(real(a)+real(b))/2,(imag(a)+imag(b))/2};}
pt addi(pt a,pt b){return pt{(real(a)+real(b)),(imag(a)+imag(b))};}
pt subt(pt a,pt b){return pt{(real(a)-real(b)),(imag(a)-imag(b))};}
pt mult(pt a,lf b){return pt{real(a)*b,imag(a)*b};}
lf abs(pt a){return sqrtl(powl(real(a),2)+powl(imag(a),2));}
lf dist(pt a,pt b){return sqrtl(powl(real(a)-real(b),2)+powl(imag(a)-imag(b),2));}
struct circ{pt p;lf r;};
const lf inf=1e18;
circ basis0(){return {pt{0,0},inf};}
circ basis1(circ a){return a;}
circ basis2(circ a,circ b)
{
pt aa=a.p,bb=b.p;
pt ab=subt(bb,aa),ba=subt(aa,bb);
lf aab=abs(ab),aba=abs(ba);
real(ab)/=aab;
imag(ab)/=aab;
real(ba)/=aba;
imag(ba)/=aba;
pt ar=addi(aa,mult(ab,a.r)),br=addi(bb,mult(ba,b.r));
return {midp(ar,br),dist(ar,br)/2.0L};
}
circ basis3(circ a,circ b,circ c)
{
lf x1=real(a.p),y1=imag(a.p),r1=a.r;
lf x2=real(b.p),y2=imag(b.p),r2=b.r;
lf x3=real(c.p),y3=imag(c.p),r3=c.r;
lf a2=x1-x2,a3=x1-x3,b2=y1-y2,b3=y1-y3,c2=r2-r1,c3=r3-r1;
lf d1=x1*x1+y1*y1-r1*r1,d2=d1-x2*x2-y2*y2+r2*r2,d3=d1-x3*x3-y3*y3+r3*r3;
lf ab=a3*b2-a2*b3;
lf xa=(b2*d3-b3*d2)/(ab*2)-x1;
lf xb=(b3*c2-b2*c3)/ab;
lf ya=(a3*d2-a2*d3)/(ab*2)-y1;
lf yb=(a2*c3-a3*c2)/ab;
lf A=xb*xb+yb*yb-1;
lf B=2*(r1+xa*xb+ya*yb);
lf C=xa*xa+ya*ya-r1*r1;
lf r=-(A?(B-sqrtl(B*B-4*A*C))/(2*A):C/B);
return {pt{x1+xa+xb*r,y1+ya+yb*r},r};
}
bool viol(circ p,circ a)
{
return a.r<p.r+dist(a.p,p.p);
}
circ solve(vector<circ>&v)
{
lf _nan=nan("aaa");
mt19937_64 mt(1999999);
//shuffle(begin(v),end(v),mt);
vector<circ>basis;
auto trivial=[&]()
{
if(size(basis)==0)return basis0();
if(size(basis)==1)return basis[0];
if(size(basis)==2)return basis2(basis[0],basis[1]);
return basis3(basis[0],basis[1],basis[2]);
};
int counter=0;
auto rec=[&](auto rec,int n)->circ
{
if(n==0||size(basis)==3)return trivial();
counter++;
if(counter>30*size(v))return {pt{0,0},_nan};
auto c=rec(rec,n-1);
auto p=v[n-1];
if(!viol(c,p)||isnan(c.r))return c;
basis.push_back(p);
c=rec(rec,n-1);
basis.pop_back();
return c;
};
auto c=rec(rec,size(v));
while(isnan(c.r))
{
counter=0;
shuffle(begin(v),end(v),mt);
c=rec(rec,size(v));
}
return c;
}
int main()
{
cin.tie(0)->sync_with_stdio(0);
int n;cin>>n;
vector<circ>v(n);
for(auto&[p,r]:v)
{
ll x,y,rr;cin>>x>>y>>rr;
p={x,y};r=rr;
}
auto ans=solve(v);
cout<<setprecision(16)<<fixed<<real(ans.p)<<" "<<imag(ans.p)<<" "<<ans.r<<"\n";
}
