JZOJ4720. 地下的太阳

题目大意

有$n*m$个火炉，排成$m$层，每层恰好$n$个火炉。
初始时第一层第$i$个炉的火力为$i$，然后火力从上而下传递。设$energy[i][j]$为第$i$层，第$j$个火炉的能量，那么满足

$$energy[i][j] = \sum_{k=1}^{n}(energy[i-1][k]+j^i)$$
然后从当前最后一层的有火力的火炉中，把编号最前的 $s$个火炉和最后的$s$个火炉单独取出（取过的火炉没有火力，即不会再被取），设火炉火力值集合分别为 $A$和$B$，然后把 $A$和$B$集合的火炉分别按原来编号升序排序，得到的能量为:
$$\sum_{i=1}^{s} \sum_{j=1}^{s}(~A_i * B_j * (|i-j|+1)^2~)$$
一共进行 $T$操作，每次操作输出最后得到的火力，答案对$1004535809$取模。

Data Constraint
$T\leq10 , n\leq10^{18} , m\leq5000 , s\leq30000$

题解

先化简一下式子

$$energy[i][j] = \sum_{k=1}^{n}(energy[i-1][k]+j^i)=n*j^i+\sum_{k=1}^{n}energy[i-1][k]$$
令 $$Sum_i=\sum_{j=1}^{n}energy[i][j]$$
则 $$energy[i][j] = n*j^i+Sum_{i-1}$$
又 $$Sum_i=\sum_{j=1}^{n}energy[i][j]=\sum_{j=1}^{n}(n*j^i+Sum_{i-1})=n*(Sum_{i-1}+\sum_{j=1}^{n}j^i)$$
然后我们可以 $O(m^2)$预处理 $S_i=\sum_{j=1}^{n}j^i$，计算出 $Sum_m$，再 $O(slogn)$取出 $A,B$的元素。

---

自然数幂和

计算$S_i=\sum_{j=1}^{n}j^i$就是求自然数幂和，这里简单讲一种$O(n^2)$的方法。
已知

$$(a+b)^{k+1}=\sum_{i=0}^{k+1}C_{k+1}^i*a^ib^{k+1-i}$$
$$\therefore (n+1)^{k+1}-n^{k+1}=C_{k+1}^1n^k+C_{k+1}^2n^{k-1}+...+C_{k+1}^kn+1$$
分别令 $n=1,2,3,...,n$代入后将 $n$个式子相加，得到
$$(n+1)^{k+1}-1=C_{k+1}^1\sum_{i=1}^n i^k+C_{k+1}^2\sum_{i=1}^n i^{k-1}+...+C_{k+1}^k\sum_{i=1}^n i+n$$
移项得 $$\sum_{i=1}^n i^k=\frac{1}{k+1}[(n+1)^{k+1}-(C_{k+1}^2\sum_{i=1}^n i^{k-1}+...+C_{k+1}^k\sum_{i=1}^n i+n+1)]$$
记 $S(n,k)=\sum_{i=1}^n i^k$，原式即为
$$S(n,k)=\frac{1}{k+1}[(n+1)^{k+1}-(C_{k+1}^2S(n,k-1)+...+C_{k+1}^kS(n,1)+n+1)]$$
然后就能递归求解了。

---

现在取出了$A$和$B$，剩下的问题就是如何快速统计答案。观察式子：

$$\sum_{i=1}^{s} \sum_{j=1}^{s}(~A_i * B_j * (|i-j|+1)^2~)$$
$$(|i-j|+1)^2=(i-j)^2+1+2|i-j|=i^2+j^2-2ij+1+2|i-j|$$
前四项的计算都比较简单，以 $i^2$为例：
将式子变形：
$$\sum_{i=1}^{s}\sum_{j=1}^{s}A_i*B_j*i^2=\sum_{i=1}^{s}A_i*i^2\sum_{j=1}^{s}B_j$$
这样一来就比较容易计算了，维护对应次数的和即可，其他三项同理。
剩下的绝对值，先讨论 $i\geq j$的情况。此时
$$\sum_{i=1}^{s}\sum_{j=1}^{i}A_i*B_j*2(i-j)$$
拆开括号有两项： $2i-2j$，以 $2i$为例
$$\sum_{i=1}^{s}A_i*i^2\sum_{j=1}^{i}B_j$$
维护 $B$的前缀和就很好计算了。
$j\geq i$的情况也同理。

SRC

#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<iostream>
#include<algorithm>
using namespace std ;

#define N 5000 + 10
#define M 30000 + 10
typedef long long ll ;
const int MO = 1004535809 ;

ll fac[N] , _fac[N] ;
ll S[N] , Sum[N] ;
ll A[M] , B[M] ;
ll sA[M] , sa1[M] , sB[M] , sb1[M] ;
int Case , m , len ;
ll n , ans ;

ll Power( ll x , int k ) {
    ll s = 1 ;
    x %= MO ;
    while ( k ) {
        if ( k & 1 ) s = s * x % MO ;
        x = x * x % MO ;
        k /= 2 ;
    }
    return s ;
}

void Pre() {
    fac[0] = _fac[0] = 1 ;
    for (int i = 1 ; i <= 5000 ; i ++ ) {
        fac[i] = fac[i-1] * i % MO ;
        _fac[i] = Power( fac[i] , MO - 2 ) ;
    }
}

ll Calc( int n , int m ) { return fac[n] * _fac[m] % MO * _fac[n-m] % MO ; }

void GetS( int k ) {
    if ( k == 1 ) {
        S[k] = n % MO * ((n + 1) % MO) % MO * _fac[2] % MO ;
        return ;
    }
    GetS( k - 1 ) ;
    S[k] = Power( (n + 1) % MO , k + 1 ) ;
    ll sum = (n + 1) % MO ;
    for (int i = 2 ; i <= k ; i ++ )
        sum = (sum + Calc(k+1,i) * S[k+1-i] % MO) % MO ;
    S[k] = (S[k] - sum + MO) % MO ;
    S[k] = S[k] * Power( k + 1 , MO - 2 ) % MO ;
}

int main() {
    scanf( "%lld%d" , &n , &m ) ;
    scanf( "%d" , &Case ) ;
    Pre() ;
    GetS(m) ;
    Sum[1] = S[1] ;
    for (int i = 2 ; i <= m ; i ++ )
        Sum[i] = (n % MO) * ((Sum[i-1] + S[i]) % MO) % MO ;
    int last = 0 ;
    while ( Case -- ) {
        A[0] = B[0] = ans = 0 ;
        scanf( "%d" , &len ) ;
        for (ll i = last + 1 ; i <= last + len ; i ++ ) A[++A[0]] = (Sum[m-1] + n % MO * Power(i,m) % MO) % MO ;
        for (ll i = n - last - len + 1 ; i <= n - last ; i ++ )
            B[++B[0]] = (Sum[m-1] + n % MO * Power(i,m) % MO) % MO ;
        ll suma = 0 , suma1 = 0 , suma2 = 0 ;
        ll sumb = 0 , sumb1 = 0 , sumb2 = 0 ;
        for (int i = 1 ; i <= len ; i ++ ) {
            suma = (suma + A[i]) % MO ;
            suma1 = (suma1 + A[i] * i % MO) % MO ;
            suma2 = (suma2 + A[i] * i % MO * i % MO) % MO ;
            sumb = (sumb + B[i]) % MO ;
            sumb1 = (sumb1 + B[i] * i % MO) % MO ;
            sumb2 = (sumb2 + B[i] * i % MO * i % MO) % MO ;
            sA[i] = (sA[i-1] + A[i]) % MO ;
            sa1[i] = (sa1[i-1] + 2ll * A[i] % MO * i % MO) % MO ;
            sB[i] = (sB[i-1] + B[i]) % MO ;
            sb1[i] = (sb1[i-1] + 2ll * B[i] % MO * i % MO) % MO ;
        }
        ans = (ans + suma * sumb % MO) % MO ;
        ans = (ans + suma * sumb2 % MO) % MO ;
        ans = (ans + suma2 * sumb % MO) % MO ;
        ans = (ans - 2ll * suma1 % MO * sumb1 % MO + MO) % MO ;
        for (int i = 1 ; i <= len ; i ++ ) {
            ans = (ans + 2ll * i * A[i] % MO * sB[i] % MO) % MO ;
            ans = (ans - A[i] * sb1[i] % MO + MO) % MO ;
            ll tp = (sB[len] - sB[i] + MO) % MO ;
            ans = (ans - 2ll * i * A[i] % MO * tp % MO + MO) % MO ;
            tp = (sb1[len] - sb1[i] + MO) % MO ;
            ans = (ans + A[i] * tp % MO) % MO ;
        }
        printf( "%lld\n" , ans ) ;
        last += len ;
    }
    return 0 ;
}

以上.