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);
    }
}