【The 13th Chinese Northeast Collegiate Programming Contest】I. Temperature Survey
题目描述
[题目链接](https://codeforces.com/gym/102220/problem/I)
给定长度为 $n$ 的 $a$ 序列，保证 $a_n \le n$，求有多少个长度为 $n$ 的 $b$ 序列，满足 $b_{i-1} \le b_i \le a_i$，答案模 $998244353$
数据范围：$1 \le n \le 2 \times 10^5$
题解
显然 $a$ 序列也是满足 $a_{i-1} \le a_i$，否则在 $a_{i-1} > a_i$ 的时候把 $a_{i-1}$ 变为 $a_i$
考虑模拟题意，设 $f_{i,j}$ 表示考虑完 $b_1,b_2,\cdots,b_i$，其中 $b_i=j$ 的方案数
显然有 $f_{i,j}=\sum_{k=0}^{j} f_{i-1,k}=f_{i-1,j}+f_{i,j-1}$
换句话说，相当于一开始有 $b_i=[i=1]$，依次做 $n$ 次前缀和，每次把 $b$ 的前 $a_i$ 个元素做一次前缀和（然而这个转化是做不了的似乎）
从另一个角度看，相当于从 $(n,j)$ 走到 $(1,1)$ 的方案数，其中 $1 \le j \le a_n$
然而这个还是太麻烦，不如再加上 $a_{n+1}=n+1$，这样就是 $(n+1,n+1)$ 走到 $(1,1)$ 的方案数了
为了方便起见，把这个 $a$ 序列翻转一下，也就是从 $(1,n+1)$ 走到 $(n+1,1)$ 的方案数
考虑分治，对于区间 $[l,r]$，就提取出来一个高度为 $a_{mid}$ 的极大矩形，然后递归做即可
考虑实现函数 $rect(l,r,bot,top)$，表示解决矩形 $[l,r] \times [bot,top]$ 的答案（即解决右边界和下边界）
设 $n=height,m=width$，且上边界的答案为 $a$，左边界的答案为 $b$，那么一共就如下几种情况：
如果是上边界到下边界，设生成函数为 $P(x)$，则：
$$
\begin{aligned}
P(x)=\sum_{i=0}^{m-1}\sum_{j=0}^{i} a_j {i-j+n-1 \choose n-1} x^i
\end{aligned}
$$
如果是左边界到右边界，设生成函数为 $Q(x)$，则：
$$
\begin{aligned}
Q(x)=\sum_{i=0}^{n-1}\sum_{j=0}^{i} b_j {i-j+m-1 \choose m-1} x^i
\end{aligned}
$$
如果是左边界到下边界，考虑到 $i-j$ 可能会小于 $0$，不妨先把它平移一下，设原先的生成函数为 $B(x)$，则平移后为 $x^{n-1}B(x)$，即：
$$
\begin{aligned}
&B(x)=\sum_{i=0}^{m-1}\sum_{j=0}^{n-1}b_j {i+n-1-j \choose i} x^i \\
\Rightarrow
&x^{n-1}B(x)=\sum_{i=n-1}^{n+m-2}\sum_{j=0}^{n-1} b_j {i-j \choose i-n+1} x^i \\
\Rightarrow
&x^{n-1}B(x)=\sum_{i=n-1}^{n+m-2}\sum_{j=0}^{n-1} b_j {i-j \choose i-n+1} x^i \\
\Rightarrow
&x^{x-1}B(x)=\sum_{i=n-1}^{n+m-2}\sum_{j=0}^{n-1}b_j \frac{(i-j)!}{(i-n+1)!(n-1-j)!} x^i \\
\Rightarrow
&x^{x-1}B(x)=\sum_{i=n-1}^{n+m-2} \frac{x^i}{(i-n+1)!} \sum_{j=0}^{i} \frac{b_j}{(n-1-j)!} (i-j)!
\end{aligned}
$$
如果是上边界到右边界，依然考虑 $i-j$ 可能会小于 $0$，因此把它平移一下，设原先的生成函数为 $A(x)$，平移后就是 $x^{m-1}A(x)$，即：
$$
\begin{aligned}
&A(x)=\sum_{i=0}^{n-1}\sum_{j=0}^{m-1} a_j {i+m-1-j \choose i} x^i \\
\Rightarrow
&x^{m-1}A(x)=\sum_{i=m-1}^{n+m-2}\sum_{j=0}^{m-1} a_j {i-j \choose i-m+1} x^i \\
\Rightarrow
&x^{m-1}A(x)=\sum_{i=m-1}^{n+m-2}\sum_{j=0}^{m-1}a_j \frac{(i-j)!}{(i-m+1)!(m-1-j)!}x^i \\
\Rightarrow
&x^{m-1}A(x)=\sum_{i=m-1}^{n+m-2} \frac{x^i}{(i-m+1)!} \sum_{j=0}^{m-1} \frac{a_j}{(m-1-j)!} (i-j)!
\end{aligned}
$$
这些部分都可以多项式乘法做完，于是这一部分的时间复杂度是 $O((n+m) \log (n+m))$
考虑到每次分治过程中，同一层的矩形互不相交，于是 $n+m$ 之和是 $a$ 的长度级别的，因此同一层的时间复杂度为 $O(|a| \log |a|)$，由于一共有 $O(\log |a|)$ 层，因此总的时间复杂度为 $O(|a| \log^2 |a|)$
代码
1 #include <bits/stdc++.h> 2 using namespace std; 3 typedef long long ll; 4 5 #pragma GCC target ("avx2") 6 7 #define __AVX__ 1 8 #define __AVX2__ 1 9 #define __SSE__ 1 10 #define __SSE2__ 1 11 #define __SSE2_MATH__ 1 12 #define __SSE3__ 1 13 #define __SSE4_1__ 1 14 #define __SSE4_2__ 1 15 #define __SSE_MATH__ 1 16 #define __SSSE3__ 1 17 18 19 #include <immintrin.h> 20 21 #include "bits/stdc++.h" 22 23 using namespace std; 24 25 typedef unsigned char U8; 26 typedef int I32; 27 typedef unsigned int U32; 28 typedef long long I64; 29 typedef unsigned long long U64; 30 31 const int P=998244353; 32 const int inv2=(P+1)>>1; 33 const int G=3; 34 const int GInv=332748118; 35 const int MaxExp=25; 36 const int MAXN=1048576<<2; 37 const int MAXK=500005; 38 const int BmM=288737297; 39 const int BmW=29; 40 41 #define maxl (524300 * 2) 42 typedef long long ll; 43 using namespace std; 44 45 46 int r[maxl],len,l; 47 U32 A[maxl],B[maxl]; 48 const ll p=998244353,_g[2]={3,332748118}; 49 ll pw(ll a,ll b){ 50 // printf("a=%lld b=%lld\n",a,b); 51 ll ret=1; 52 while(b){ 53 if(b&1) 54 ret=ret*a%p; 55 a=a*a%p,b>>=1; 56 } 57 return ret; 58 } 59 void modify(int x){ 60 for(len=1,l=0;len<=x;len<<=1,l++); 61 for(int i=0;i<len;i++) 62 r[i]=(r[i>>1]>>1)|((i&1)<<(l-1)); 63 } 64 65 namespace NTT{ 66 67 68 U32 GcdEx(U32 A, U32 B, I32& x, I32& y) { 69 if (!B) { 70 x = 1; 71 y = 0; 72 return A; 73 } 74 U32 d = GcdEx(B, A % B, y, x); 75 y -= x * (I32) (A / B); 76 return d; 77 } 78 79 inline U32 MNorm(I32 V) { 80 V = V % P; 81 return (U32) (V < 0 ? V + P : V); 82 } 83 84 inline U32 MAdd(U32 A, U32 B) { 85 U32 res = A + B; 86 return res < P ? res : res - P; 87 } 88 89 inline U32 MSub(U32 A, U32 B) { 90 U32 res = A - B; 91 return A < B ? res + P : res; 92 } 93 94 inline U32 MMul(U32 A, U32 B) { 95 return (U32) ((U64) A * B % P); 96 } 97 98 inline U32 MPow(U32 A, U32 B) { 99 U32 res = 1; 100 while (B) { 101 if (B & 1) 102 res = MMul(res, A); 103 A = MMul(A, A); 104 B >>= 1; 105 } 106 return res; 107 } 108 109 inline U32 MInv(U32 N) { 110 I32 x, y; 111 GcdEx(N, P, x, y); 112 x %= P; 113 return (U32) (x < 0 ? x + P : x); 114 } 115 116 inline __m256i VLod(const U32* __restrict__ A) { 117 return _mm256_load_si256((const __m256i*) A); 118 } 119 120 inline void VSto(U32* __restrict__ A, __m256i v) { 121 _mm256_store_si256((__m256i*) A, v); 122 } 123 124 inline __m256i VEx0(__m256i v) { 125 const __m256i vm0 = _mm256_set_epi64x( 126 0x111111111b1a1918, 0x1111111113121110, 127 0x111111110b0a0908, 0x1111111103020100 128 ); 129 return _mm256_shuffle_epi8(v, vm0); 130 } 131 132 inline __m256i VEx1(__m256i v) { 133 const __m256i vm1 = _mm256_set_epi64x( 134 0x111111111f1e1d1c, 0x1111111117161514, 135 0x111111110f0e0d0c, 0x1111111107060504 136 ); 137 return _mm256_shuffle_epi8(v, vm1); 138 } 139 140 inline __m256i VIntlv(__m256i v0, __m256i v1) { 141 return _mm256_blend_epi32(v0, _mm256_shuffle_epi32(v1, 0xb1), 0xaa); 142 } 143 144 inline __m256i VAdd(__m256i va, __m256i vb) { 145 const __m256i vm32 = _mm256_set1_epi32(P); 146 __m256i vra = _mm256_add_epi32(va, vb); 147 __m256i vrb = _mm256_sub_epi32(vra, vm32); 148 return _mm256_min_epu32(vra, vrb); 149 } 150 151 inline __m256i VSub(__m256i va, __m256i vb) { 152 const __m256i vm32 = _mm256_set1_epi32(P); 153 __m256i vra = _mm256_sub_epi32(va, vb); 154 __m256i vrb = _mm256_add_epi32(vra, vm32); 155 return _mm256_min_epu32(vra, vrb); 156 } 157 158 inline __m256i VMul(__m256i va0, __m256i va1, __m256i vb0, __m256i vb1) { 159 const __m256i vm32 = _mm256_set1_epi32(P); 160 const __m256i vm64 = _mm256_set1_epi64x(P); 161 const __m256i vbmm = _mm256_set1_epi64x(BmM); 162 __m256i vmul0 = _mm256_mul_epi32(va0, vb0); 163 __m256i vmul1 = _mm256_mul_epi32(va1, vb1); 164 __m256i vlow = VIntlv(vmul0, vmul1); 165 __m256i vquo0 = _mm256_srli_epi64(_mm256_mul_epi32(_mm256_srli_epi64(vmul0, 29), vbmm), BmW); 166 __m256i vquo1 = _mm256_srli_epi64(_mm256_mul_epi32(_mm256_srli_epi64(vmul1, 29), vbmm), BmW); 167 __m256i vval0 = _mm256_mul_epi32(vquo0, vm64); 168 __m256i vval1 = _mm256_mul_epi32(vquo1, vm64); 169 __m256i vval = VIntlv(vval0, vval1); 170 __m256i vra = _mm256_sub_epi32(vlow, vval); 171 __m256i vrb = _mm256_add_epi32(vra, vm32); 172 __m256i vrc = _mm256_sub_epi32(vra, vm32); 173 __m256i vmin = _mm256_min_epu32(vra, vrb); 174 return _mm256_min_epu32(vmin, vrc); 175 } 176 177 inline __m256i VMul(__m256i va, __m256i vb0, __m256i vb1) { 178 return VMul(VEx0(va), VEx1(va), vb0, vb1); 179 } 180 181 inline __m256i VMul(__m256i va, __m256i vb) { 182 return VMul(va, VEx0(vb), VEx1(vb)); 183 } 184 185 inline void VMul(U32* __restrict__ A, U32 Len, U32 W) { 186 if (Len < 8) { 187 for (U32 i = 0; i < Len; ++i) 188 A[i] = MMul(A[i], W); 189 return; 190 } 191 __m256i vw = _mm256_set1_epi64x(W); 192 for (U32 i = 0; i < Len; i += 8) 193 VSto(A + i, VMul(VLod(A + i), vw, vw)); 194 } 195 196 inline void VMul(U32* __restrict__ A, const U32* __restrict__ B, U32 Len) { 197 if (Len < 8) { 198 for (U32 i = 0; i < Len; ++i) 199 A[i] = MMul(A[i], B[i]); 200 return; 201 } 202 for (U32 i = 0; i < Len; i += 8) 203 VSto(A + i, VMul(VLod(A + i), VLod(B + i))); 204 } 205 206 inline void VSqr(U32* __restrict__ A, U32 Len) { 207 if (Len < 8) { 208 for (U32 i = 0; i < Len; ++i) 209 A[i] = MMul(A[i], A[i]); 210 return; 211 } 212 for (U32 i = 0; i < Len; i += 8) { 213 __m256i va = VLod(A + i); 214 __m256i v0 = VEx0(va); 215 __m256i v1 = VEx1(va); 216 VSto(A + i, VMul(v0, v1, v0, v1)); 217 } 218 } 219 220 U32 WbFwd[MaxExp + 1]; 221 U32 WbInv[MaxExp + 1]; 222 U32 LenInv[MaxExp + 1]; 223 224 inline void NttInitAll(int Max) { 225 for (int Exp = 0; Exp <= Max; ++Exp) { 226 WbFwd[Exp] = MPow(G, (P - 1) >> Exp); 227 WbInv[Exp] = MPow(GInv, (P - 1) >> Exp); 228 LenInv[Exp] = MInv(1u << Exp); 229 } 230 } 231 232 inline void NttImpl1(U32* __restrict__ A, U32 Len) { 233 for (U32 j = 0; j < Len; j += 2) { 234 U32 a0 = MAdd(A[j + 0], A[j + 1]); 235 U32 b0 = MSub(A[j + 0], A[j + 1]); 236 A[j + 0] = a0; 237 A[j + 1] = b0; 238 } 239 } 240 241 inline void NttFwd2(U32* __restrict__ A, U32 Len, U32 Wn) { 242 for (U32 j = 0; j < Len; j += 4) { 243 U32 a0 = MAdd(A[j + 0], A[j + 2]); 244 U32 a1 = MAdd(A[j + 1], A[j + 3]); 245 U32 b0 = MSub(A[j + 0], A[j + 2]); 246 U32 b1 = MSub(A[j + 1], A[j + 3]); 247 A[j + 0] = a0; 248 A[j + 1] = a1; 249 A[j + 2] = b0; 250 A[j + 3] = MMul(b1, Wn); 251 } 252 } 253 254 inline void NttFwd3(U32* __restrict__ A, U32 Len, U32 Wn) { 255 U32 W2 = MMul(Wn, Wn); 256 U32 W3 = MMul(W2, Wn); 257 const __m128i vm32 = _mm_set1_epi32(P); 258 for (U32 j = 0; j < Len; j += 8) { 259 __m128i va = _mm_load_si128((const __m128i*) (A + j)); 260 __m128i vb = _mm_load_si128((const __m128i*) (A + j + 4)); 261 __m128i vc = _mm_add_epi32(va, vb); 262 __m128i vd = _mm_sub_epi32(va, vb); 263 __m128i ve = _mm_sub_epi32(vc, _mm_andnot_si128(_mm_cmpgt_epi32(vm32, vc), vm32)); 264 __m128i vf = _mm_add_epi32(vd, _mm_and_si128(_mm_cmpgt_epi32(vb, va), vm32)); 265 _mm_store_si128((__m128i*) (A + j), ve); 266 _mm_store_si128((__m128i*) (A + j + 4), vf); 267 A[j + 5] = MMul(Wn, A[j + 5]); 268 A[j + 6] = MMul(W2, A[j + 6]); 269 A[j + 7] = MMul(W3, A[j + 7]); 270 } 271 } 272 273 inline void NttFwd(U32* __restrict__ A, int Exp) { 274 U32 Len = 1u << Exp; 275 U32 Wn = WbFwd[Exp]; 276 for (int i = Exp - 1; i >= 3; --i) { 277 U32 ChkSiz = 1u << i; 278 U32 tw2 = MMul(Wn, Wn); 279 U32 tw3 = MMul(tw2, Wn); 280 U32 tw4 = MMul(tw3, Wn); 281 U32 tw5 = MMul(tw4, Wn); 282 U32 tw6 = MMul(tw5, Wn); 283 U32 tw7 = MMul(tw6, Wn); 284 U32 twn = MMul(tw7, Wn); 285 __m256i vw32 = _mm256_set_epi32(tw7, tw6, tw5, tw4, tw3, tw2, Wn, 1); 286 __m256i vwn = _mm256_set1_epi64x(twn); 287 for (U32 j = 0; j < Len; j += 2u << i) { 288 U32* A_ = A + j; 289 U32* B_ = A_ + ChkSiz; 290 __m256i vw = vw32; 291 for (U32 k = 0; k < ChkSiz; k += 8) { 292 __m256i va = VLod(A_ + k); 293 __m256i vb = VLod(B_ + k); 294 __m256i vw0 = VEx0(vw); 295 __m256i vw1 = VEx1(vw); 296 __m256i vc = VAdd(va, vb); 297 __m256i vd = VSub(va, vb); 298 VSto(A_ + k, vc); 299 VSto(B_ + k, VMul(vd, vw0, vw1)); 300 vw = VMul(vw0, vw1, vwn, vwn); 301 } 302 } 303 Wn = MMul(Wn, Wn); 304 } 305 if (Exp >= 3) { 306 NttFwd3(A, Len, Wn); 307 Wn = MMul(Wn, Wn); 308 } 309 if (Exp >= 2) 310 NttFwd2(A, Len, Wn); 311 if (Exp) 312 NttImpl1(A, Len); 313 } 314 315 inline void NttInv2(U32* __restrict__ A, U32 Len, U32 Wn) { 316 for (U32 j = 0; j < Len; j += 4) { 317 U32 a0 = A[j + 0]; 318 U32 a1 = A[j + 1]; 319 U32 b0 = A[j + 2]; 320 U32 b1 = MMul(A[j + 3], Wn); 321 A[j + 0] = MAdd(a0, b0); 322 A[j + 1] = MAdd(a1, b1); 323 A[j + 2] = MSub(a0, b0); 324 A[j + 3] = MSub(a1, b1); 325 } 326 } 327 328 inline void NttInv3(U32* __restrict__ A, U32 Len, U32 Wn) { 329 U32 W2 = MMul(Wn, Wn); 330 U32 W3 = MMul(W2, Wn); 331 const __m128i vm32 = _mm_set1_epi32(P); 332 for (U32 j = 0; j < Len; j += 8) { 333 A[j + 5] = MMul(Wn, A[j + 5]); 334 A[j + 6] = MMul(W2, A[j + 6]); 335 A[j + 7] = MMul(W3, A[j + 7]); 336 __m128i va = _mm_load_si128((const __m128i*) (A + j)); 337 __m128i vb = _mm_load_si128((const __m128i*) (A + j + 4)); 338 __m128i vc = _mm_add_epi32(va, vb); 339 __m128i vd = _mm_sub_epi32(va, vb); 340 __m128i ve = _mm_sub_epi32(vc, _mm_andnot_si128(_mm_cmpgt_epi32(vm32, vc), vm32)); 341 __m128i vf = _mm_add_epi32(vd, _mm_and_si128(_mm_cmpgt_epi32(vb, va),vm32)); 342 _mm_store_si128((__m128i*) (A + j), ve); 343 _mm_store_si128((__m128i*) (A + j + 4), vf); 344 } 345 } 346 347 inline void NttInv(U32* __restrict__ A, int Exp) { 348 if (!Exp) 349 return; 350 U32 Len = 1u << Exp; 351 NttImpl1(A, Len); 352 if (Exp == 1) { 353 VMul(A, Len, LenInv[1]); 354 return; 355 } 356 U32 Ws[MaxExp]; 357 Ws[0] = WbInv[Exp]; 358 for (int i = 1; i < Exp; ++i) 359 Ws[i] = MMul(Ws[i - 1], Ws[i - 1]); 360 NttInv2(A, Len, Ws[Exp - 2]); 361 if (Exp == 2) { 362 VMul(A, Len, LenInv[2]); 363 return; 364 } 365 NttInv3(A, Len, Ws[Exp - 3]); 366 if (Exp == 3) { 367 VMul(A, Len, LenInv[3]); 368 return; 369 } 370 for (int i = 3; i < Exp; ++i) { 371 U32 ChkSiz = 1u << i; 372 U32 Wn = Ws[Exp - 1 - i]; 373 U32 tw2 = MMul(Wn, Wn); 374 U32 tw3 = MMul(tw2, Wn); 375 U32 tw4 = MMul(tw3, Wn); 376 U32 tw5 = MMul(tw4, Wn); 377 U32 tw6 = MMul(tw5, Wn); 378 U32 tw7 = MMul(tw6, Wn); 379 U32 twn = MMul(tw7, Wn); 380 __m256i vw32 = _mm256_set_epi32(tw7, tw6, tw5, tw4, tw3, tw2, Wn, 1); 381 __m256i vwn = _mm256_set1_epi64x(twn); 382 for (U32 j = 0; j < Len; j += 2u << i) { 383 U32* A_ = A + j; 384 U32* B_ = A_ + ChkSiz; 385 __m256i vw = vw32; 386 for (U32 k = 0; k < ChkSiz; k += 8) { 387 __m256i vw0 = VEx0(vw); 388 __m256i vw1 = VEx1(vw); 389 __m256i vb = VMul(VLod(B_ + k), vw0, vw1); 390 vw = VMul(vw0, vw1, vwn, vwn); 391 __m256i va = VLod(A_ + k); 392 __m256i vc = VAdd(va, vb); 393 __m256i vd = VSub(va, vb); 394 VSto(A_ + k, vc); 395 VSto(B_ + k, vd); 396 } 397 } 398 } 399 VMul(A, Len, LenInv[Exp]); 400 } 401 402 inline int Log2Ceil(U32 N) { 403 static const U8 Table[32] = { 404 0, 9, 1, 10, 13, 21, 2, 29, 405 11, 14, 16, 18, 22, 25, 3, 30, 406 8, 12, 20, 28, 15, 17, 24, 7, 407 19, 27, 23, 6, 26, 5, 4, 31, 408 }; 409 N = (N << 1) - 1; 410 N |= N >> 1; 411 N |= N >> 2; 412 N |= N >> 4; 413 N |= N >> 8; 414 N |= N >> 16; 415 return (int) Table[(N * 0x07c4acddu) >> 27]; 416 } 417 } 418 419 using namespace NTT; 420 421 const int mod = 998244353, N = 5e6 + 10; 422 423 int n, a[N]; 424 425 ll getinv(ll n) { 426 return pw(n, mod - 2); 427 } 428 429 ll fac[N], invfac[N], inv[N]; 430 void init(int n) { 431 fac[0] = invfac[0] = inv[0] = 1; 432 for(int i = 1 ; i <= n ; ++ i) { 433 fac[i] = fac[i - 1] * i % mod; 434 inv[i] = getinv(i); 435 } 436 invfac[n] = pw(fac[n], mod - 2); 437 for(int i = n - 1 ; i ; -- i) { 438 invfac[i] = invfac[i + 1] * (i + 1) % mod; 439 } 440 } 441 442 void upd(int &x, int y) { 443 x = ((ll) x + y) % mod; 444 } 445 446 void mul(vector<int> &a, vector<int> &b) { 447 // 计算这两个序列的卷积 448 vector<int> c(a.size() + b.size() - 1); 449 // for(int i = 0 ; i < a.size() ; ++ i) { 450 // for(int j = 0 ; j < b.size() ; ++ j) { 451 // upd(c[i + j], (ll) a[i] * b[j] % mod); 452 // } 453 // } 454 // a = c; return ; 455 456 int len = 1; while(len <= c.size() * 2) len <<= 1; 457 int lg2 = Log2Ceil(len); 458 459 // printf("%d, %d %d, %d\n", len, a.size(), b.size(), c.size()); 460 461 for(int i = 0 ; i < len ; ++ i) { 462 A[i] = B[i] = 0; 463 } 464 for(int i = 0 ; i < a.size() ; ++ i) { 465 A[i] = a[i]; 466 } 467 for(int i = 0 ; i < b.size() ; ++ i) { 468 B[i] = b[i]; 469 } 470 NttFwd(A, lg2), NttFwd(B, lg2); 471 for(int i = 0 ; i < len ; ++ i) { 472 A[i] = (ll) A[i] * B[i] % mod; 473 } 474 NttInv(A, lg2); 475 for(int i = 0 ; i < c.size() ; ++ i) { 476 c[i] = A[i]; 477 } 478 a = c; 479 } 480 ll C(int n, int m) { 481 return n < m || m < 0 ? 0 : fac[n] * invfac[m] % mod * invfac[n - m] % mod; 482 } 483 484 map<int, int> dp[N]; 485 486 void sol_rect(int l, int r, int bot, int top) { 487 // 左右边界，下边界，上边界 488 489 // printf("处理： [%d, %d] - [%d, %d]\n", l, r, bot, top); 490 491 if(l == 1 && r == 1 && top == n + 1) { 492 // puts("VIS"); 493 // 就是这玩意，起始位置，左上角 494 for(int i = bot ; i <= top ; ++ i) { 495 dp[l][i] = 1; 496 // printf("dp[%d][%d] = %d\n", l, i, dp[l][i]); 497 } 498 return ; 499 } 500 501 int width = r - l + 1, height = top - bot + 1; 502 503 vector<int> bottom(width), right(height); // 底下和右侧的答案 504 vector<int> a(width); // 上面的 505 vector<int> b(height); // 左面的 506 for(int i = 0 ; i < width ; ++ i) { 507 a[i] = dp[l + i][top + 1]; 508 } 509 for(int i = 0 ; i < height ; ++ i) { 510 b[i] = dp[l - 1][top - i]; 511 } 512 513 // printf("a: "); 514 // for(int x: a) printf("%d ", x); puts(""); 515 // printf("b: "); 516 // for(int x: b) printf("%d ", x); puts(""); 517 518 { 519 // left -> bottom 520 // for(int i = 0 ; i < width ; ++ i) { 521 // for(int j = 0 ; j < height ; ++ j) { 522 // upd(bottom[i], (ll) b[j] * C(i + height - 1 - j, i) % mod); 523 // } 524 // } 525 526 vector<int> x(height), y(width + height - 1); 527 for(int i = 0 ; i < height ; ++ i) { 528 x[i] = (ll) b[i] * invfac[height - 1 - i] % mod; 529 } 530 for(int i = 0 ; i < width + height - 1 ; ++ i) { 531 y[i] = fac[i]; 532 } 533 mul(x, y); 534 for(int i = 0 ; i < width ; ++ i) { 535 upd(bottom[i], (ll) x[height - 1 + i] * invfac[i] % mod); 536 } 537 } 538 539 { 540 // top -> bottom 541 // for(int i = 0 ; i < width ; ++ i) { 542 // for(int j = 0 ; j <= i ; ++ j) { 543 // upd(bottom[i], (ll) a[j] * C(i - j + height - 1, height - 1) % mod); 544 // } 545 // } 546 vector<int> x(width), y(width); 547 for(int i = 0 ; i < width ; ++ i) { 548 x[i] = a[i]; 549 y[i] = C(i + height - 1, height - 1); 550 } 551 mul(x, y); 552 for(int i = 0 ; i < width ; ++ i) { 553 upd(bottom[i], x[i]); 554 } 555 } 556 557 { 558 // top -> right 559 // for(int i = 0 ; i < height ; ++ i) { 560 // for(int j = 0 ; j < width ; ++ j) { 561 // upd(right[i], (ll) a[j] * C(i + width - 1 - j, i) % mod); 562 // } 563 // } 564 vector<int> x(width), y(width + height - 1); 565 for(int i = 0 ; i < width ; ++ i) { 566 x[i] = (ll) a[i] * invfac[width - 1 - i] % mod; 567 } 568 for(int i = 0 ; i < width + height - 1 ; ++ i) { 569 y[i] = fac[i]; 570 } 571 mul(x, y); 572 for(int i = 0 ; i < height ; ++ i) { 573 upd(right[i], (ll) x[width - 1 + i] * invfac[i] % mod); 574 } 575 } 576 577 { 578 // left -> right 579 // for(int i = 0 ; i < height ; ++ i) { 580 // for(int j = 0 ; j <= i ; ++ j) { 581 // upd(right[i], (ll) b[j] * C(i - j + width - 1, width - 1) % mod); 582 // } 583 // } 584 vector<int> x(height), y(height); 585 for(int i = 0 ; i < height ; ++ i) { 586 x[i] = b[i]; 587 y[i] = C(i + width - 1, width - 1); 588 } 589 mul(x, y); 590 for(int i = 0 ; i < height ; ++ i) { 591 upd(right[i], x[i]); 592 } 593 } 594 595 // printf("right: "); 596 // for(int x: right) printf("%d ", x); puts(""); 597 // printf("bottom: "); 598 // for(int x: bottom) printf("%d ", x); puts(""); 599 600 for(int i = l ; i <= r ; ++ i) { 601 upd(dp[i][bot], bottom[i - l]); 602 // printf("dp_1[%d][%d] = %d\n", i, bot, dp[i][bot]); 603 } 604 for(int i = top ; i >= bot ; -- i) { 605 606 // 由于右下角会被算两次，所以这里先鸽一次 607 if(i == bot) continue; 608 609 upd(dp[r][i], right[top - i]); 610 // printf("dp_2[%d][%d] = %d\n", r, i, dp[r][i]); 611 } 612 } 613 614 void sol(int l, int r, int bot) { 615 // 左右边界，底边的高度 616 if(l > r) { 617 return ; 618 } 619 int mid = (l + r) >> 1; 620 int x = mid, y = mid; 621 while(x - 1 >= l && a[x - 1] == a[mid]) -- x; 622 while(y + 1 <= r && a[y + 1] == a[mid]) ++ y; 623 sol(l, x - 1, a[mid] + 1); 624 625 // printf("mid = %d\n", mid); 626 // printf("x = %d, y = %d\n", x, y); 627 sol_rect(l, y, bot, a[mid]); 628 sol(y + 1, r, bot); 629 } 630 631 int main() { 632 NttInitAll(22); 633 init(N - 1); 634 int t; scanf("%d", &t); 635 while(t --) { 636 scanf("%d", &n); 637 for(int i = 1 ; i <= n ; ++ i) { 638 scanf("%d", &a[i]); 639 } 640 a[n + 1] = n + 1; 641 reverse(a + 1, a + 1 + n + 1); 642 for(int i = 0 ; i <= n + 1 ; ++ i) { 643 dp[i].clear(); 644 } 645 646 // for(int i = 1 ; i <= n + 1 ; ++ i) { 647 // printf("%d ", a[i]); 648 // } puts(""); 649 650 sol(1, n + 1, 1); 651 printf("%d\n", ((ll) dp[n + 1][1] % mod + mod) % mod); 652 } 653 }
