「SDOI2014」数表

「SDOI2014」数表

设$\sigma(x)$为$x$的约束和，求

$$\sum_{i=1}^n \sum_{j=1}^m \sum_{k=1}^a k[\sigma(gcd(i, j)) = k]$$

数据组数$<= 20000$
$n, m <= 100000, a<=10^9$

---

很明显$\sigma(gcd(i, j))=k$并不好求，所以我们将他拆开

$$ans = \sum_{k=1}^{lim} \sigma(k) \bigg( \sum_{i=1}^n \sum_{j=1}^m [gcd(i, j) = k] \bigg)$$

很明显后面括号内的东西我们可以反演
设

$$f(k) = \sum_{i=1}^n \sum_{j=1}^m [gcd(i, j) = k]$$
$$g(k) = \sum_{i=1}^n \sum_{j=1}^m [gcd(i, j) \%k = 0] = \lfloor\dfrac{n}{k}\rfloor \lfloor \dfrac{m}{k} \rfloor$$
$$\therefore g(x) = \sum_{x|d} f(d)$$
$$\Rightarrow f(x) = \sum_{x|d} \mu\big(\dfrac{d}{x}\big) g(d)$$
$$f(x)= \sum_{x|d} \mu\big(\dfrac{d}{x}\big) \lfloor\dfrac{n}{d}\rfloor \lfloor \dfrac{m}{d} \rfloor$$

带入原式中

$$ans = \sum_{k=1}^{lim} \sigma(k) \bigg( \sum_{k|d}\mu\big(\dfrac{d}{k}\big)\lfloor\dfrac{n}{d}\rfloor \lfloor \dfrac{m}{d} \rfloor \bigg)$$

交换一下

$$ans = \sum_{k=1}^{lim} \sum_{k|d} \sigma(k) \mu\big(\dfrac{d}{k}\big)\lfloor\dfrac{n}{d}\rfloor \lfloor \dfrac{m}{d} \rfloor$$

改为枚举$d$与$d$的约数$k$

$$ans = \sum_{d=1}^{lim} \sum_{k|d} \sigma(k) \mu\big(\dfrac{d}{k}\big)\lfloor\dfrac{n}{d}\rfloor \lfloor \dfrac{m}{d} \rfloor$$

交换一下

$$ans = \sum_{d=1}^{lim} \lfloor\dfrac{n}{d}\rfloor \lfloor \dfrac{m}{d} \rfloor \bigg( \sum_{k|d} \sigma(k) \mu\big(\dfrac{d}{k}\big) \bigg)$$

我们可以线性筛求出$\sigma(x), \mu(x)$
我们每次要求的是$\sigma(x) <= a, x <= max\{gcd(n, m)\}$的$ans$
将询问按$a$从小到大排序，$\sigma(x)$按值从小到大排序，每次询问完处理$a_{last} < \sigma(x) <= a_{this}$ 的$\sigma(x)$。

可以遇见的是，我们最多要处理$n$个$\sigma(x)$，每次查询的复杂度是$O(\sqrt{n} log n)$。

所以总复杂度为$O(nlog^2n + q\sqrt{n}logn)$
(并不知道为什么处理$\sigma(x)$的复杂度是$O(nlog^2n)$，我觉得应该是$O(n\sqrt{n}logn)$呀。。求大神解答)

code：

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N = 200000;
const int M = 120000;

bool is_prime[N];
int prime[N], tot = 0, mobius[N], max_n = 0;
ll num[N], qsum[N], dsum[N];

inline ll ksm( ll a, int b )
{
    ll ret = 1; 
    while( b )
    {
        if( b & 1 ) ret *= a;
        b >>= 1; a *= a;
    }
    return ret;
}

void init( int size )
{
    memset( is_prime, true, sizeof( is_prime ) );
    is_prime[1] = false; mobius[1] = 1; 
    dsum[1] = num[1] = qsum[1] = 1;

    for( int i = 2; i <= size; i ++ )
    {
        if( is_prime[i] )
        {
            prime[++tot] = i;
            num[i] = 1; mobius[i] = -1;    
            dsum[i] = qsum[i] = i+1;
        }
        for( int j = 1; j <= tot; j ++ )
        {
            int k = i * prime[j];
            if( k > size ) break;
            is_prime[k] = false;

            if( i % prime[j] == 0 )
            {
                mobius[k] = 0;
                num[k] = num[i] + 1;
                qsum[k] = qsum[i] + ksm( prime[j], num[i] + 1 );
                dsum[k] = dsum[i] / qsum[i] * qsum[k];
                break;
            }

            mobius[k] = - mobius[i];
            num[k] = 1; qsum[k] = prime[j] + 1;
            dsum[k] = dsum[i] * dsum[prime[j]];
        }
    }
}

struct tQuery {
    int n, m, a, q_id; ll ans;

//    bool operator < ( const tQuery &rhs ) const { return a < rhs.a; }
} q[N];

struct tDsum {
    int num; ll sum;
    inline void init( int _num, ll _sum ) { num = _num; sum = _sum; }

    bool operator < ( const tDsum &rhs ) const
    {
//      if( num > max_n ) return false;
//      if( rhs.num > max_n ) return true;
        return sum < rhs.sum;
    }
} d[N];

bool cmp_a( tQuery x, tQuery y ) { return x.a < y.a; }
bool cmp_id( tQuery a, tQuery b ) { return a.q_id < b.q_id; }  

int now_p = 1;

int t[N];
void add(int x,int val)
{
    for(int i=x;i<=max_n;i+=i&-i)t[i]+=val;
}
int query(int x)
{
    int tmp=0;
    for(int i=x;i;i-=i&-i)tmp+=t[i];
    return tmp;
}

inline void build_seg( int ed )
{   
    while( d[now_p].sum <= ed && now_p < M )
    { 
        if( d[now_p].num > max_n ) { now_p ++; continue; }

//      printf( "now_p.sum:%d\n", d[now_p].sum );
        for( int i = d[now_p].num; i <= max_n; i += d[now_p].num )
        {
            int tmp = mobius[i/d[now_p].num];
            add( i, d[now_p].sum * tmp );  
//            seg.modify( 1, 1, M, i, d[now_p].sum * tmp );
//          printf( "[ADD] j:%d, val:%d\n", i, d[now_p].sum * tmp );
        }
        now_p ++;
    }  
}

int main()
{
//  freopen( "g8.in", "r", stdin ); /
//  freopen( "g8.out", "w", stdout );

    int ttt; scanf( "%d", &ttt );
    for( int i = 1; i <= ttt; i ++ )
    {
        scanf( "%d%d%d", &q[i].n, &q[i].m, &q[i].a );
        q[i].q_id = i; if( q[i].n > q[i].m ) swap( q[i].n, q[i].m );
        max_n = max( max_n, q[i].m );
    }

    init( M ); 

//    for( int i = 1;i <= 20; i ++ )
//      printf( "%lld ", dsum[i] ); printf( "dsum\n" );

    for( int i = 1; i <= M; i ++ )
        d[i].init( i, dsum[i] );

    sort( d+1, d+M+1 );
//    printf( "ok" );
    sort( q+1, q+ttt+1, cmp_a );
//    printf( "ok2" );

    for( int i = 1; i <= ttt; i ++ )
    {
//      printf( "i:%d\n", i );
        build_seg( q[i].a ); 

        int lim = min( q[i].n, q[i].m ); 
        int n = q[i].n, m = q[i].m, next = 0;

        for( int j = 1; j <= lim; j = next + 1 )
        {
//          printf( "j:%d\n", j ); 
            next = min( n/(n/j), m/(m/j) ); if( next > lim ) next = lim;
            q[i].ans += 1ll * (n/j) * (m/j) * ( query( next ) - query( j-1 ) );
//            printf( "n/j:%d, m/j:%d, seg:%d, ans:%lld\n", n/j, m/j, seg.query( 1, 1, M, j, next ), q[i].ans );
        }
    }

    sort( q+1, q+ttt+1, cmp_id );
    for( int i = 1; i <= ttt; i ++ )
        printf( "%lld\n", q[i].ans % (1ll<<31) );

    return 0;
}