[CQOI2015]选数

题目链接

问题分析

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

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

我们令

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

就有

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

所以

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

[Ans=f(k)=sum_{i=1}^{lfloorfrac{H}{k} floor}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(lfloorfrac{n}{i} floor)\ ]

令(g=id_0)，那么

[S(n)=1-sum_{i=2}^nS(lfloorfrac{n}{i} floor) ]

其中(*)代表狄利克雷卷积，(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
", 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;
}