Dirichlet 前缀和

发布于 # algorithm

和后缀和。

0 序

bk=ikaib_k = \sum_{i|k}a_i

bk=kdadb_k = \sum_{k|d}a_d

都一样,没啥区别。

1 内容

前缀和的写法:

for( auto &p: primes )
    for( int j = 1; 1LL * j * p <= N; j ++ )
        b[ j * p ] += a[j];

后缀和翻过来就行。

2 证明

一种比较野蛮的办法是参见 Sum over Subsets DP 的证明,这个东西实际上就是以质数为轴的高维前缀和。

另一种是循环不等式,没看懂。

这里考虑证明前缀和的形式,后缀和同理:

pip_i 表示从小到大第 ii 个质数,假设 Si,j={xx=kjk=i=1ppiai(aiN)}S_{i,j} = \{ x | x = k \cdot j \land k = \sum_{i=1}^p p_i^{a_i} ( a_i \in N )\}

我们令在第 kk 轮结束后,bi=jSi,jaib_i = \sum_{j \in S_{i,j}} a_i

k=0k = 0 显然成立。

k0k \neq 0 时,如果 pkip_k | i ,令 i=ipki' = \frac{i}{p_k},则有

bi=bi+bi=jSk1,iaj+jSk,ikajSk1,iSk,i=,Sk1,iSk,i=Sk,ibi=bi+bi=jSk,iaj\begin{aligned} b_i &= b_i + b_{i'} \\ &= \sum_{j \in S_{ k - 1, i }} a_j + \sum_{j \in S_{k,\frac{i}{k}}} a_j \\ &\because S_{ k - 1, i} \cap S_{ k, i' } = \emptyset, S_{ k - 1, i } \cup S_{k,i'}=S_{k,i} \\ &\therefore b_i = b_i + b_{i'} = \sum_{j \in S_{k,i}} a_j \\ \end{aligned}

则如果在 k1(k>0)k-1(k > 0) 时成立,则在 kk 时成立。

容易发现这么做复杂度和埃氏筛一样,O(nloglogn)O(n \log \log n)

3 例 Codeforces 1614 D Divan and Kostomuksha

题目大意

给定数组 aa,求任意重排后最大化 i=1nf(i)\sum_{i=1}^n f(i),其中 f(1)=a1,f(i)=gcd(f(i1),ai)f(1) = a_1, f(i) = \gcd( f( i - 1 ), a_i )

n1e5n \leq 1e5。两档,ai5×106,ai2×107a_i \leq 5 \times 10^6, a_i \leq 2 \times 10^7

做法

p_i = \sum_j \[ a_j \mod i = 0 \], gig_i 表示当前构造序列中所有数字的 gcd=ki(kN)\gcd = ki (k \in N) 时,能够得到的最大 if(i)\sum_i f(i)m=maxam =\max{a}

容易得到: gi=maxidfd+(pipd)×ig_i = \max_{i|d} f_d + ( p_i - p_d ) \times i

很容易得到 imi=O(nlog(n))\sum_i \frac{m}{i} = O(n\log(n)) 的做法,能过第一档,考虑优化。

瓶颈在两处,求 pip_igig_i

pip_i 直接使用 Dirichlet 后缀。

注意到可以 gig_i 递推式可以改善:

sis_i 表示为 ii 的质因子集合,则有 gi=maxjsifd+(pipd)×ig_i = \max_{j \in s_i} f_d + ( p_i - p_d ) \times i,其中 d=ijd = \frac{i}{j}

容易发现可以直接套用埃氏筛。

总复杂度 O(nloglogn)O(n \log \log n)

Code

/*
 * d2.cpp
 * Copyright (C) 2022 Woshiluo Luo <[email protected]>
 *
 * 「Two roads diverged in a wood,and I—
 * I took the one less traveled by,
 * And that has made all the difference.」
 *
 * Distributed under terms of the GNU AGPLv3+ license.
 */

#include <cstdio>
#include <cstdlib>
#include <cstring>

#include <vector>
#include <algorithm>

typedef const int cint;
typedef long long ll;
typedef unsigned long long ull;

inline bool isdigit( const char cur ) { return cur >= '0' && cur <= '9'; }/*{{{*/
template <class T>
T Max( T a, T b ) { return a > b? a: b; }
template <class T>
T Min( T a, T b ) { return a < b? a: b; }
template <class T>
void chk_Max( T &a, T b ) { if( b > a ) a = b; }
template <class T>
void chk_Min( T &a, T b ) { if( b < a ) a = b; }
template <typename T>
T read() {
    T sum = 0, fl = 1;
    char ch = getchar();
    for (; isdigit(ch) == 0; ch = getchar())
        if (ch == '-') fl = -1;
    for (; isdigit(ch); ch = getchar()) sum = sum * 10 + ch - '0';
    return sum * fl;
}
template <class T>
T pow( T a, int p ) {
    T res = 1;
    while( p ) {
        if( p & 1 )
            res = res * a;
        a = a * a;
        p >>= 1;
    }
    return res;
}/*}}}*/

const int N = 2e7 + 10;

int pool[N];
ll f[N];

std::vector<int> primes;
bool is_prime[N];

void pre() {
    memset( is_prime, true, sizeof(is_prime) );
    for( int i = 2; i < N; i ++ ) {
        if( is_prime[i] )
            primes.push_back(i);
        for( auto &x: primes ) {
            if( 1LL * x * i >= N )
                break;
            is_prime[ x * i ] = false;
            if( i % x == 0 )
                break;
        }
    }
}

int main() {
#ifdef woshiluo
    freopen( "d2.in", "r", stdin );
    freopen( "d2.out", "w", stdout );
#endif
    pre();

    cint n = read<int>();
    for( int i = 1; i <= n; i ++ )
        pool[ read<int>() ] ++;
    for( auto &p: primes ) {
        for( int j = ( N - 1 ) / p; j >= 1; j -- ) {
            pool[j] += pool[ j * p ];
        }
    }

    ll ans = 0;
    for( int i = N - 1; i >= 1; i -- ) {
        chk_Max( f[i], 1LL * pool[i] * i );
        for( auto &p: primes ) {
            if( 1LL * i * p >= N )
                break;
            cint nxt = i * p;
            chk_Max( f[i], f[nxt] + 1LL * ( pool[i] - pool[nxt] ) * i );
        }
        chk_Max( ans, f[i] );
    }

    printf( "%lld\n", ans );
}
Comment seems to stuck. Try to refresh?✨

Woshiluo's NoteBook

「Jump up HIGH!!」