HDU 6051 - If the starlight never fade | 2017 Multi-University Training Contest 2

/*
HDU 6051 - If the starlight never fade [ 原根,欧拉函数 ]  |  2017 Multi-University Training Contest 2
题意：给定 m,p, p 是素数设 f(i) 是 满足 (x+y)^i ≡ x^i mod p 的 (x,y) 对数 且 1 ≤ x ≤ p-1 , 1 ≤ y ≤ m 求 ∑[1≤i≤p-1] i*f(i)限制： m ≤ p-1, 2 ≤ p ≤ 1e9
分析：设 g 为 p 的原根，则x,y可表示为 x = g^a, y = g^b(x+y)^i ≡ x^i (mod p)(g^a + g^b)^i ≡ g^ai (mod p)(1 + g^(b-a))^i ≡ 1 (mod p)设 g^k = 1 + g^(b-a)，则 g^ki ≡ 1 (mod p)则 k 满足 ki % (p-1) == 0 ，即 k 是 (p-1)/gcd(i, p-1) 的倍数 由于 0 < k < p-1 , 则k的取值有  (p-1) / ((p-1)/gcd(i, p-1)) - 1 = gcd(i, p-1)-1 个回带 1 + y/x = g^k x = y * (g^k-1)^(-1)x = y * (g^k-1)^(φ(p)-1)则 当y固定时， x, k 一一对应，x的取值也有 gcd(i, p-1)-1 个ans = ∑[1≤i≤p-1] i*f(i)= ∑[1≤i≤p-1] i * m * (gcd(i, p-1)-1)= m * ( ∑[1≤i≤p-1] i * gcd(i, p-1) - p*(p-1)/2)求解 ∑[1≤i≤n] i * gcd(i, n)= ∑[1≤i≤n] i ∑[k|n] k * [gcd(i, n) == k]= ∑[k|n] ∑[1≤i≤n] i * k * [gcd(i, n) == k]= ∑[k|n] k^2 ∑[1≤i≤n/k] i * [gcd(i, n/k) == 1]求解 ∑[1≤i≤n] i * [gcd(i, n) == 1]= (∑[1≤i≤n] i * [gcd(i, n) == 1] + ∑[1≤i≤n] (n-i) * [gcd(i, n-i) == 1]) / 2= ∑[1≤i≤n] (i+n-i)/2 * [gcd(i, n) == 1]= (n * φ(n) + [n==1]) / 2
*/
#include <bits/stdc++.h>
using namespace std;
#define LL long long
const LL MOD = 1e9+7;
LL phi(LL n) {LL ans = n;for (LL i = 2; i * i <= n; i++) {if (n % i == 0) {ans -= ans / i;while (n % i == 0) n /= i;}}if (n > 1) ans -= ans/n;return ans;
}
LL Cal(LL x, LL n)
{LL res = 1;res *= ( (n/x)*phi(n/x) + bool(n/x == 1) ) / 2;res %= MOD;res *= x*x % MOD;res %= MOD;return res;
}
LL p, m;
int main()
{int t; scanf("%d", &t);for (int tt = 1; tt <= t; tt++){scanf("%lld%lld", &m, &p);LL ans = 0;for (LL i = 1; i*i <= p-1; i++){if (i*i == p-1){ans += Cal(i, p-1);ans %= MOD;}else if ((p-1)%i == 0){ans += Cal(i, p-1) + Cal((p-1)/i, p-1);ans %= MOD;}}ans += MOD -  p*(p-1)/2 % MOD;ans = ans % MOD * m % MOD;printf("Case #%d: %lld\n", tt, ans);}
}