【ZOJ】3881 From the ABC conjecture【暴力容斥】

传送门：【ZOJ】3881 From the ABC conjecture

复杂度大概O(N0.67)O(N^{0.67})，我也不会算www，首先转换一下（我们是根据积性函数打表找规律得到的，也可以推出来）使得：

g(N)=∏pi ϵ N(pia+1)g(N)=\prod_{pi~\epsilon ~N} (pi^a+1)

暴力展开发现贡献为：

h(N)=∑x=1N∑y=1N(x⋅y≤N and gcd(x,y)=1)⋅xh(N)=\sum_{x=1}^{N}\sum_{y=1}^{N} ( x \cdot y \leq N ~and~ gcd(x,y)=1) \cdot x

然后是一个暴力dfs容斥的过程：

ans=∑x=1N∑y=1N(x⋅y≤N)⋅x−∑k=1N√∑y=1N(x⋅y≤N and gcd(x,y)=k)⋅xans=\sum_{x=1}^{N}\sum_{y=1}^{N} ( x \cdot y \leq N ) \cdot x - \sum_{k=1}^{\sqrt N}\sum_{y=1}^{N}(x \cdot y \leq N ~and~gcd(x,y)=k)\cdot x

=∑x=1N∑y=1N(x⋅y≤N)⋅x−∑k=1N√k∑y=1N/k2(x/k⋅y/k≤N/k2 and gcd(x/k,y/k)=1)⋅x/k=\sum_{x=1}^{N}\sum_{y=1}^{N} ( x \cdot y \leq N ) \cdot x - \sum_{k=1}^{\sqrt N}k\sum_{y=1}^{N/k^2}(x/k \cdot y/k \leq N/k^2 ~and~gcd(x/k,y/k)=1)\cdot x/k

可以发现右边是一个原式的递归式，于是我们可以暴力枚举N−−√\sqrt N项来暴力，当然计算的值需要记忆化一下。

my code:my~~code:

[code]#include <stdio.h>
#include <string.h>
#include <map>
#include <algorithm>
using namespace std ;

typedef long long LL ;

#define clr( a , x ) memset ( a , x , sizeof a )

const int MAXN = 20000005 ;
const int mod = 1e9 + 9 ;
const int v2 = 5e8 + 5 ;

map < LL , int > mp ;
map < LL , int > :: iterator it ;
bool vis[MAXN] ;
int dp[MAXN] ;
LL n ;

int dfs ( LL n , LL ans = 0 ) {
    for ( LL i = 1 , j ; i <= n ; i = j + 1 ) {
        j = min ( n , n / ( n / i ) ) ;
        LL ni = n / i , x = i + j , y = j - i + 1 ;
        if ( ni >= mod ) ni %= mod ;
        if ( x >= mod ) x %= mod ;
        if ( y >= mod ) y %= mod ;
        ans += ni * x % mod * y % mod * v2 % mod ;
    }
    for ( LL i = 2 ; i * i <= n ; ++ i ) {
        LL x = n / ( i * i ) ;
        if ( x < MAXN ) {
            if ( !vis[x] ) {
                vis[x] = true ;
                dp[x] = dfs ( x ) ;
            }
            ans -= dp[x] * i % mod ;
        } else {
            it = mp.find ( x ) ;
            if ( it != mp.end () ) ans -= it->second * i % mod ;
            else {
                int tmp = dfs ( x ) ;
                mp.insert ( map < LL , int > :: value_type ( x , tmp ) ) ;
                ans -= tmp * i % mod ;
            }
        }
    }
    return ( ans % mod + mod ) % mod ;
}

int main () {
    while ( ~scanf ( "%lld" , &n ) ) printf ( "%d\n" , dfs ( n ) ) ;
    return 0 ;
}