[CQOI2015]选数

题目链接

问题分析

题意还是很好理解的，就是求

\[Ans=f(k)=\sum_{a_1=L}^H\sum_{a_2=L}^H\cdots\sum_{a_N=L}^H[\gcd(a_1,a_2,\cdots,a_N)=K] \]

我们令

\[F(k)=\sum_{a_1=L}^H\sum_{a_2=L}^H\cdots\sum_{a_N=L}^H[K|\gcd(a_1,a_2,\cdots,a_N)]=Num^N(k)\\ (Num(k)=\lfloor\frac{H}{k}\rfloor-\max(0,\lfloor\frac{L-1}{k}\rfloor)) \]

就有

\[F(k)=\sum_{k|d}f(d) \]

所以

\[f(k)=\sum_{k|d}\mu(\frac{d}{k})F(d) \]

\[Ans=f(k)=\sum_{i=1}^{\lfloor\frac{H}{k}\rfloor}\mu(i)Num^N(ik) \]

\(Num\)可以整除分块，而\(\mu\)由于数据范围大需要使用杜教筛。

考虑

\[f(n)=\mu(n)\\ S(n)=\sum_{i=1}^nf(i)\\ g(1)S(n)=\sum_{i=1}^n(g*f)(i)-\sum_{i=2}^ng(i)S(\lfloor\frac{n}{i}\rfloor)\\ \]

令\(g=id_0\)，那么

\[S(n)=1-\sum_{i=2}^nS(\lfloor\frac{n}{i}\rfloor) \]

其中\(*\)代表狄利克雷卷积，\(id_0(k)=k^0=1\)。

这样就好了，时间复杂度\(O(n^\frac{2}{3})\)。

参考程序

这个写法不优良，时间复杂度是 \(O(n^{\frac{3}{4}})\) 的。还要带上一个 map 的大常数 log ?

#include <bits/stdc++.h>
#define LL long long
using namespace std;

const LL INF = 10000000000000000;
const LL Threshold = 2000000;
const LL Mod = 1000000007;
LL Mu[ Threshold + 10 ], Vis[ Threshold + 10 ], Prime[ Threshold + 10 ], Num, Sum[ Threshold + 10 ];
LL N, K, L, H, Ans;
map< LL, LL > Map;

void EulerSieve();
void Cal();

int main() {
    EulerSieve();
    scanf( "%lld%lld%lld%lld", &N, &K, &L, &H );
    Cal();
    printf( "%lld\n", Ans );
    return 0;
}

void EulerSieve() {
    Mu[ 1 ] = 1;
    for( LL i = 2; i <= Threshold; ++i ) {
        if( !Vis[ i ] ) {
            Mu[ i ] = -1;
            Prime[ ++Num ] = i;
        }
        for( LL j = 1; j <= Num && i * Prime[ j ] <= Threshold; ++j ) {
            Vis[ i * Prime[ j ] ] = 1;
            if( i % Prime[ j ] == 0 ) break;
            Mu[ i * Prime[ j ] ] = -Mu[ i ];
        }
    }
    for( LL i = 1; i <= Threshold; ++i ) Sum[ i ] = ( Sum[ i - 1 ] + Mu[ i ] ) % Mod;
    return;
}

LL Power( LL k, LL t ) {
    if( t == 0 ) return 1;
    LL T = Power( k, t >> 1 );
    T = T * T % Mod;
    if( t & 1 ) T = k % Mod * T % Mod;
    return T;
}

LL S( LL t ) {
    if( t <= Threshold ) return Sum[ t ];
    if( Map.find( t ) != Map.end() ) return Map[ t ];
    LL Ans = 1;
    for( LL i = 2, j; i <= t; i = j + 1 ) {
        j = t / ( t / i );
        Ans = ( ( Ans - ( j - i + 1 ) % Mod * S( t / i ) % Mod ) % Mod + Mod ) % Mod;
    }
    Map[ t ] = Ans;
    return Ans;
}

void Cal() {
    LL l = ( L - 1 ) / K, r = H / K; 
    for( LL i = 1, j; i <= r; i = j + 1 ) {
        j = min( ( l / i ) ? l / ( l / i ) : INF, r / ( r / i  ) );
        Ans = ( Ans + ( ( S( j ) - S( i - 1 ) ) % Mod + Mod ) % Mod * Power( r / i - l / i, N ) % Mod ) % Mod;
    }
    return;
}