Dirichlet 前缀和

和后缀和。

0 序

$$
b_k = \sum_{i|k}a_i
$$

$$
b_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 的证明,这个东西实际上就是以质数为轴的高维前缀和。

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

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

令 $p_i$ 表示从小到大第 $i$ 个质数,假设 $S_{i,j} = \{ x | x = k \cdot j \land k = \sum_{i=1}^p p_i^{a_i} ( a_i \in N )\}$

我们令在第 $k$ 轮结束后,$b_i = \sum_{j \in S_{i,j}} a_i$。

在 $k = 0$ 显然成立。

在 $k \neq 0$ 时,如果 $p_k | i$ ,令 $i’ = \frac{i}{p_k}$,则有

$$
\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}
$$

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

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

3 例 Codeforces 1614 D Divan and Kostomuksha

题目大意

给定数组 $a$,求任意重排后最大化 $\sum_{i=1}^n f(i)$,其中 $f(1) = a_1, f(i) = \gcd( f( i – 1 ), a_i )$。

$n \leq 1e5$。两档,$a_i \leq 5 \times 10^6, a_i \leq 2 \times 10^7$。

做法

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

容易得到: $g_i = \max_{i|d} f_d + ( p_i – p_d ) \times i$

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

瓶颈在两处,求 $p_i$ 和 $g_i$。

$p_i$ 直接使用 Dirichlet 后缀。

注意到可以 $g_i$ 递推式可以改善:

令 $s_i$ 表示为 $i$ 的质因子集合,则有 $g_i = \max_{j \in s_i} f_d + ( p_i – p_d ) \times i$,其中 $d = \frac{i}{j}$。

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

总复杂度 $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 );
}

发表评论

您的电子邮箱地址不会被公开。

message
account_circle
Please input name.
email
Please input email address.
links

此站点使用Akismet来减少垃圾评论。了解我们如何处理您的评论数据