2017 Multi-University Training Contest - 第二场 07 If the starlight never fade （数论）

题目链接：
HDU 6051

题意：
$f(i)$ 表示满足$(x+y)^{i}≡x^{i}\% p$ 中 $(x,y)$ 的对数，其中，$1≤x≤p−1,1≤y≤m.$
给你一个非负数 $m$ 和素数 $p$，让你求$\sum_{i=1}^{p-1}if\left(i\right)$，答案$\mod1e9+7$。

题解：
数学题。

令 $g$ 为 $p$ 的原根，$x=g^{a},y=g^{b}$。
那么有$\left(g^{a} + g^{b} \right)^{i} = g^{ai} \% p$。即$\left(1 + g^{b-a} \right)^{i} = 1 \% p$。
令$1 + g^{b-a} = g^{k}$，则有$g^{ki}=1 \% p$，即$ki \% \left(p - 1 \right) = 0$。

因为 $0 < k < p - 1$，所以 $k$ 的取值有$gcd\left(i,p - 1 \right) - 1$个，递推回去，可知$x = y * inv \left(g^{k} - 1 \right)$,故$x$的取值也是$gcd\left(i,p - 1 \right) - 1$个。

所以我们最后只需要求$m \sum_{i=1}^{p-1} \left\{ i*gcd(i,p-1)-i \right\}$。
其中$\sum_{i=1}^{p-1}i$，我们可以直接求出。

而另一部分$\sum_{i=1}^{p-1}i*gcd(i,p-1)$，我们可以通过枚举$gcd$的值，化为$\sum_{d|(p-1)}d^{2}\sum_{k=1}^{\frac{p-1}{d}}k*\left [ gcd\left ( k,\frac{p-1}{d} \right )=1 \right ]$。
然后根据公式$\sum_{i=1}^{n} i*\left [ gcd\left ( i,n \right )=1\right ]= \left ( n\varphi \left ( n \right )+\left [ n=1 \right ]\right )/2$求解。
其中，$\varphi \left ( n \right )$为欧拉函数。

AC代码：

#include<bits/stdc++.h>
const int mod = 1e9+7;
using namespace std;
typedef long long ll;

int mod_mul(int a,int b)
{
    ll c = 1LL * a * b;
    return c - c / mod * mod;
}

int euler(int n) 
{
    int ret = n;
    int i;
    for(i = 2 ; i * i <= n ; i ++)
    {
        if(n % i == 0)
        {
            ret = ret - ret / i;
            while(n % i == 0)
            {
                n /= i;
            }
        }
    }
    if(n > 1)
    {
        ret = ret - ret / n;
    }
    return ret;
}

int solve(int n)
{
    int ans = 0;
    int tmp = 0;
    for(int i = 1 , j ; i * i <= n; i++)
    {
        j = n / i;
        if(i * j == n)
        {
            tmp = mod_mul( (1LL * j * euler(j) + ( j == 1 ) ) / 2 % mod, mod_mul( i , i ) );
            if((ans += tmp) >= mod) ans -= mod;


            if(i != j)
            {
                tmp = mod_mul( (1LL * i * euler(i) + ( i == 1 ) ) / 2 % mod, mod_mul( j , j ) );
                if((ans += tmp) >= mod) ans -= mod;
            }

        }
    }
    ans -= mod_mul( mod_mul( n, n + 1), 500000004);
    if(ans < 0) ans += mod;
    return ans;
}
int main()
{
    int t;
    int m,p,cas=0;
    scanf("%d",&t);
    while(t--)
    {
        cas++;
        scanf("%d%d",&m,&p);
        int ans = mod_mul( solve(p - 1), m );
        printf("Case #%d: ",cas);
        printf("%d\n",ans);
    }
    return 0;
}