「LOJ 2476」蒜头的奖杯

题目

点这里看题目，简要题意见下：

给定 6 个长度均为 \(n\) 的数列 \(A,B,C,D,E,F,G\)，求：

\[\sum_{i=1}^n\sum_{j=1}^n\sum_{k=1}^nA_iB_jC_kD_{\gcd(i,j)}E_{\gcd(i,k)}F_{\gcd(j,k)} \]

对于 \(100\%\) 的数据，满足 \(1\le n\le 10^5\)，且任意数列中的任意一个数都是 \([0,10^{18}]\) 范围内的整数。

分析

变量很多，想法之一就是设置枚举变量的“优先次序”：

\[\begin{aligned} &\sum_{i=1}^n\sum_{j=1}^n\sum_{k=1}^nA_iB_jC_kD_{\gcd(i,j)}E_{\gcd(i,k)}F_{\gcd(j,k)}\\ =&\sum_{i=1}^nA_i\sum_{j=1}^nB_jD_{\gcd(i,j)}\sum_{k=1}^nC_kE_{\gcd(i,k)}F_{\gcd(j,k)}\\ \end{aligned} \]

其中的 \(\gcd\) 仍然比较难以处理，我们先将一部分用 Mobius 反演展开。设 \(E'=E*\mu,F'=F*\mu\)：

\[\newcommand{\lcm}{\operatorname{lcm}} \newcommand{\div}[2]{\lfloor\frac{#1}{#2}\rfloor} \begin{aligned} &\sum_{i=1}^nA_i\sum_{j=1}^nB_jD_{\gcd(i,j)}\sum_{k=1}^nC_kE_{\gcd(i,k)}F_{\gcd(j,k)}\\ =&\sum_{i=1}^nA_i\sum_{j=1}^nB_jD_{\gcd(i,j)}\sum_{k=1}^nC_k\sum_{p|i,p|k}E'_p\sum_{q|j,q|k}F'_q\\ =&\sum_{p=1}^n\sum_{q=1}^nE'_pF'_q\sum_{p|k,q|k,k\le n}C_k\sum_{p|i,i\le n}\sum_{q|j,j\le n}A_iB_jD_{\gcd(i,j)}\\ =&\sum_{\lcm(p,q)\le n}E'_pF'_q\sum_{\lcm(p,q)|k,k\le n}C_k\sum_{i=1}^{\div{n}{p}}\sum_{j=1}^{\div{n}{q}}A_{ip}B_{jq}D_{\gcd(ip,jq)} \end{aligned} \]

注意到满足 \(\lcm(p,q)\le n\) 的 \((p,q)\) 对是可以枚举的（？），我们修改式子的形式：

\[\begin{aligned} &\sum_{\lcm(p,q)\le n}E'_pF'_q\sum_{\lcm(p,q)|k,k\le n}C_k\sum_{i=1}^{\div{n}{p}}\sum_{j=1}^{\div{n}{q}}A_{ip}B_{jq}D_{\gcd(ip,jq)}\\ =&\sum_{d=1}^n\sum_{p=1}^{\div n d}\sum_{q=1}^{\div n {pd}}[\gcd(p,q)=1]E'_{pd}F'_{qd}\sum_{dpq|k,k\le n}C_k\sum_{i=1}^{\div{n}{pd}}\sum_{j=1}^{\div{n}{qd}}A_{ipd}B_{jqd}D_{d\gcd(ip,jq)} \end{aligned} \]

这样变换过后，所需的下标相应地 scale 了，因此我们调整一下。

设 \(m=\div n d\)，且对于 \(X=A,B,C,D,E,F\)，令 \(x_k=X_{kd}\)：

\[\begin{aligned} &\sum_{d=1}^n\sum_{p=1}^{\div n d}\sum_{q=1}^{\div n {pd}}[\gcd(p,q)=1]E'_{pd}F'_{qd}\sum_{dpq|k,k\le n}C_k\sum_{i=1}^{\div{n}{pd}}\sum_{j=1}^{\div{n}{qd}}A_{ipd}B_{jqd}D_{d\gcd(ip,jq)}\\ =&\sum_{d=1}^n\sum_{p=1}^{m}\sum_{q=1}^{\div{m}{p}}[\gcd(p,q)=1]e'_{p}f'_{q}\sum_{pq|k,k\le m}c_k \sum_{i=1}^{\div{m}{p}}\sum_{j=1}^{\div{m}{q}}a_{ip}b_{jq}d_{\gcd(ip,jq)}\end{aligned} \]

考虑复杂性的问题——理论上 \(d,p,q\) 都是可以直接枚举的。枚举完 \(d\) 之后，我们有充裕的时间去处理出 \(a,b,c,d,e',f'\)，相应地也就可以计算出 \(c^*_r=\sum_{r|k,k\le m}c_k\)。因此，主要的问题在于如何快速地处理出关于 \(i,j\) 的和式。

我们先设计比较明确的 \(p,q\) 部分：我们可以按 \(\min\{p,q\}\) 的顺序枚举 \((p,q)\) 对，这样就实际上只需要枚举 \(O(\sqrt m)\) 个 \(\min\{p,q\}\)。如果我们枚举好了 \(d\) 和 \(\min\{p,q\}\)，则这部分的时间仅有 \(O(\sum_k \sqrt{nk^{-1}})<O(n)\)。

我们尝试在知道了 \(d\) 和 \(\min\{p,q\}\) 的情况下尽快计算出后半部分。不妨将它看作是关于 \(p,q\) 的二元函数 \(g(p,q)=\sum_{i=1}^{\div m p}\sum_{j=1}^{\div m q}a_{ip}b_{jq}d_{\gcd(ip,jq)}\)，那么我们相当于知道了其中一个变量的取值，所以只需要处理好 \(g(\min\{p,q\},x)\) 和 \(g(x,\min\{p,q\})\)。

不妨来处理 \(g(\min\{p,q\},x)\) 的情况。仍然是先将 \(d\) 用 Mobius 反演改造一下，设 \(d'=d*\mu\)，那么我们将可以得到：

\[\begin{aligned} &\sum_{i=1}^{\div{m}{p}}\sum_{j=1}^{\div{m}{q}}a_{ip}b_{jq}d_{\gcd(ip,jq)}\\ =&\sum_{i=1}^{\div m p}\sum_{j=1}^{\div m q}a_{ip}b_{jq}\sum_{r|ip,r|jq}d'_{r}\\ =&\sum_{p|i,i\le m}\sum_{q|j,j\le m}\sum_{r|i,r|j}a_ib_jd'_r\\ =&\sum_{q|j}b_j\sum_{r|j}d'_r\sum_{r|i}a_i[p|i] \end{aligned} \]

唯一的变量 \(q\) 被转移到了最外层和式的枚举条件里，这样我们就可以从内到外，用 Dirichlet 卷积计算出每个 \(q\) 对应的值。对于 \(g(x,\min\{p,q\})\) 的处理同理。

再研究一下现在的复杂度，为 \(O(\sum_{d}((\frac{n}{d})^{\frac 3 2}\log\log n+\sum_p\sum_q[pq\le \frac n d]))=O(n^{\frac 3 2}\log\log n)\)，而且实现起来貌似常数还比较小。

小结：

1. 学习一下如何入手一个比较复杂的式子，比如定主元，枚举之类的；
2. 在化式子的同时，关注计算复杂性。这一部分在某些情况下可以指导化式子的方向。
3. 关注变量与已知量，调整结构将变量放在靠外的位置。

代码

#include <cstdio>

#define rep( i, a, b ) for( int i = (a) ; i <= (b) ; i ++ )
#define per( i, a, b ) for( int i = (a) ; i >= (b) ; i -- )

typedef unsigned long long ull;

const int MAXN = 1e5 + 5, MAXS = 330;

template<typename _T>
void read( _T &x ) {
    x = 0; char s = getchar(); bool f = false;
    while( s < '0' || '9' < s ) { f = s == '-', s = getchar(); }
    while( '0' <= s && s <= '9' ) { x = ( x << 3 ) + ( x << 1 ) + ( s - '0' ), s = getchar(); }
    if( f ) x = -x;
}

template<typename _T>
void write( _T x ) {
    if( x < 0 ) putchar( '-' ), x = -x;
    if( 9 < x ) write( x / 10 );
    putchar( x % 10 + '0' );
}

bool cop[MAXN][MAXS], vis[MAXN][MAXS];

ull mu[MAXN];
int prime[MAXN], pn;
bool isPrime[MAXN];

ull A[MAXN], B[MAXN], C[MAXN], D[MAXN], E[MAXN], F[MAXN];
ull A_[MAXN], B_[MAXN], C_[MAXN], D_[MAXN], E_[MAXN], F_[MAXN], tp[MAXN], tq[MAXN];

int N;

bool Chk( const int x, const int y ) {
    if( x < y ) return Chk( y, x );
    if( vis[x][y] ) return cop[x][y];
    vis[x][y] = true;
    if( ! y ) return cop[x][y] = x == 1;
    return cop[x][y] = Chk( y, x % y );
}

void EulerSieve( const int n ) {
    mu[1] = 1;
    for( int i = 2 ; i <= n ; i ++ ) {
        if( ! isPrime[i] ) 
            mu[prime[++ pn] = i] = -1;
        for( int j = 1 ; j <= pn && 1ll * i * prime[j] <= n ; j ++ ) {
            isPrime[i * prime[j]] = true;
            if( ! ( i % prime[j] ) ) break;
            mu[i * prime[j]] = - mu[i];
        }
    }
}

void Transform( ull *f, const int n ) {
    static ull tmp[MAXN];
    rep( i, 1, n ) tmp[i] = 0;
    for( int p = 1 ; p <= n ; p ++ )
        for( int q = p ; q <= n ; q += p )
            tmp[q] += f[p] * mu[q / p];
    rep( i, 1, n ) f[i] = tmp[i];
}

int main() {
    read( N ), EulerSieve( N );
    rep( i, 1, N ) read( A[i] );
    rep( i, 1, N ) read( B[i] );
    rep( i, 1, N ) read( C[i] );
    rep( i, 1, N ) read( D[i] );
    rep( i, 1, N ) read( E[i] );
    rep( i, 1, N ) read( F[i] );
    Transform( E, N );
    Transform( F, N );
    for( int i = 1 ; i <= pn ; i ++ )
        for( int j = N / prime[i] ; j ; j -- )
            C[j] += C[j * prime[i]];
    ull ans = 0;
    rep( d, 1, N ) {
        const int M = N / d;
        rep( i, 1, M )
            A_[i] = A[i * d], D_[i] = D[i * d], 
            B_[i] = B[i * d], E_[i] = E[i * d],
            C_[i] = C[i * d], F_[i] = F[i * d];
        Transform( D_, M );
        for( int p = 1 ; 1ll * p * p <= M ; p ++ ) {
            rep( i, 1, M ) 
                tq[i] = ( i % p == 0 ) * A_[i],
                tp[i] = ( i % p == 0 ) * B_[i];
            for( int i = 1 ; i <= pn && prime[i] <= M ; i ++ )
                for( int j = M / prime[i] ; j ; j -- )
                    tq[j] += tq[j * prime[i]],
                    tp[j] += tp[j * prime[i]];
            rep( r, 1, M ) tq[r] *= D_[r], tp[r] *= D_[r];
            for( int i = 1 ; i <= pn && prime[i] <= M ; i ++ )
                for( int j = 1 ; j * prime[i] <= M ; j ++ )
                    tq[j * prime[i]] += tq[j],
                    tp[j * prime[i]] += tp[j];
            rep( j, 1, M ) tq[j] *= B_[j], tp[j] *= A_[j];
            for( int i = 1 ; i <= pn && prime[i] <= M ; i ++ )
                for( int j = M / prime[i] ; j ; j -- )
                    tq[j] += tq[j * prime[i]],
                    tp[j] += tp[j * prime[i]];
            for( int q = p ; 1ll * q * p <= M ; q ++ )
                if( Chk( p, q ) ) {
                    ans += E_[p] * F_[q] * C_[p * q] * tq[q];
                    if( p ^ q ) ans += E_[q] * F_[p] * C_[p * q] * tp[q];
                }
        }
    }
    write( ans ), putchar( '\n' );
    return 0;
}