ACM-ICPC 2018 沈阳赛区网络预赛 G. Spare Tire

原题传送门

题意：
给定一个函数a(n) ，和两个整数，n, m;
求闭区间[1, n]内，所有的a(x)的和, x满足gcd（x, m) == 1

思路：
（下午打比赛的时候思路还是蛮接近的，但就是想不透那一点，赛后大佬提醒才想清楚）

首先，打个表就可以发现，a(n) = n * n + n，为了确保无误，代回去验算一下发现是满足的，然后是求[1, n]以内跟m互素的数字，由于求互素的答案比较难，那么我们可以选择计算总答案再减去不互素的答案，总答案就是a(n)的前缀和

接下来便是求[1, n]以内跟m互素的答案贡献，容斥一下，先对m进行素数分解（在这里可以预处理1e4的素数进行分解，也可以直接√n的做法分解），那么就可以得到小于等于8个素因子，根据莫比乌斯的定义，我们只需要枚举这些素因子的任意组合即可，（因为具有平方因子的数字的莫比乌斯值为0，贡献也是0，而这些素因子的任意组合是没有平方因子的），如果是奇数个素数，对答案贡献系数就是-1，偶数个素数则贡献系数为1。 最后是根据当前的素因子组合的乘积x，计算a(x) + a(2x) + a(3x) ……这个是可以直接O（1）计算的，就是式子推得有点烦

a(nx) = nx * nx + nx = x * x * n * n + n * x;
S(nx) = x * x * n * (n + 1) * (2n + 1) / 6 + x * n * (n + 1) / 2;

(一开始偷懒想用矩阵快速幂计算，就不用推式子了，结局很明显，T的明明白白）

AC代码：

#include <bits/stdc++.h>
#define pb push_back
using namespace std;
typedef long long ll;

const int N = 1e4 + 10;
const int mod = 1e9 + 7;
int cnt, prime[N];
bool vis[N];

void init()
{
    for(int i = 2; i < N ; i ++)
    {
        if(!vis[i])
            prime[cnt ++] = i;
        for(int j = 0, k; j < cnt && (k = i * prime[j]) < N; j ++)
        {
            vis[k] = 1;
            if(i % prime[j] == 0)   break;
        }
    }
}

ll kuai(ll x, ll n)
{
    ll a = 1;
    while(n)
    {
        if(n & 1)   a = a * x % mod;
        x = x * x % mod;
        n >>= 1;
    }
    return a;
}

int main()
{
    init();
    ll n, m;
    while(~scanf("%lld%lld", &n, &m))
    {
        ll ans = n * (n + 1) % mod * (n + 2) % mod * kuai(3, mod - 2) % mod;
        if(m != 1)
        {
            vector<int> v;
            for(int i = 0; i < cnt && m != 1; i ++)
            {
                if(m % prime[i] == 0)
                {
                    v.pb(prime[i]);
                    while(m % prime[i] == 0)    m /= prime[i];
                }
            }
            if(m - 1)   v.pb(m);

            for(int i = 1, k = v.size(); i < (1 << k); i ++)
            {
                ll a = 1, p = 1;
                for(int j = 0; j < k; j ++)
                    if(i & (1 << j))
                        a *= v[j], p *= -1;

                if(n < a)   continue ;

                ll x = n / a;
                ll S = a * a % mod * x  % mod * (x + 1) % mod * (2 * x + 1) % mod * kuai(6, mod - 2) % mod;
                S = (S + a * x % mod * (x + 1) % mod * kuai(2, mod - 2)) % mod;
                ans = (ans + S * p) % mod + mod;
                ans %= mod;
            }

        }
        printf("%lld\n", ans);
    }
    return 0;
}