简单的数学题 - 数论

题目大意：给定正整数nn，求有多少个整数四元组$a,b,c,d\in [0,n-1]$满足$ab=cd(\mathrm{m}\mathrm{o}\mathrm{d}n)$。由于$n$非常大，将以质因数分解的形式给出$n={\prod }_{i=1}^m{p}_i^{c_i}$。$m\le 5\times 10^5,p,c\le 10^9$。
题解：由中国剩余定理的结论知我们只需要求出$n=p^c$的答案然后乘起来即可。
考虑这个怎么做，不妨设$cn{t}_i={\sum }_{(a,b)}[ab\mathrm{m}\mathrm{o}\mathrm{d}n=i]$,那么答案就是${\sum }_{i=0}^{n-1}cn{t}_i^2$。显然$cn{t}_0=n^2-{\sum }_{i=1}^{n-1}cn{t}_i$，考虑某个$cn{t}_i(i>0)$怎么算。
不妨设$i=p^{k}q$，满足$\gcd (p,q)=1$。显然$0\le k<c,q>0$。
那么$ab=i$，等价于$a=p^{k^{\prime }}{a}^{\prime },b=p^{k-k^{\prime }}{b}^{\prime },\gcd (a^{\prime },p)=\gcd (b^{\prime },p)=1,a^{\prime }{b}^{\prime }=q(\mathrm{m}\mathrm{o}\mathrm{d}{p}^{c-k})$，并且每求出这样的一组$(a^{\prime },b^{\prime },k^{\prime })$，都会有$p^{k^{\prime }}\times p^{k-k^{\prime }}=p^k$组$(a,b,k)$与之对应，而$(a^{\prime },b^{\prime },k^{\prime })$的组数显然就是$\varphi (p^{c-k})$，与$k^{\prime }$无关，因此对于每个$k^{\prime }$，答案就是$p^k\varphi (p^{c-k})$，因此$cn{t}_i=(k+1)p^k\varphi (p^{c-k})=(k+1)\varphi (p^c)$。
然后考虑对于每个$k$有多少个$i$，显然就是$\varphi (p^{c-k})$
首先考虑${\sum }_{i=1}^{n-1}cn{t}_i^2={\sum }_{k=0}^{c-1}(k+1{)}^2{\varphi }^2(p^c)\varphi (p^{c-k})={\varphi }^2(p^c)(p-1)p^c{\sum }_{k=1}^c{k}^2{\left(\frac{1}{p}\right)}^k$
后面那个怎么求：
$$\begin{gathered} S_2=\sum _{k=1}^n{k}^2{q}^k \\ {S}_2-n^2{q}^n=\sum _{k=2}^n(k-1{)}^2{q}^{k-1} \\ q{S}_2-n^2{q}^{n+1}=\sum _{k=2}^n(k^2-2k+1)q^k \\ q{S}_2-n^2{q}^{n+1}=\sum _{k=2}^n{k}^2{q}^k-2\sum _{k=2}^{n}k{q}^k+\sum _{k=2}^n{q}^k \\ q{S}_2-n^2{q}^{n+1}=S_2-q+2(S_1-q)+S_0-q \\ {S}_2=\frac{n^2{q}^{n+1}-q+2(S_1-q)+S_0-q}{q-1} \end{gathered}$$

其中$S_1={\sum }_{k=1}^{n}k{q}^k,S_0={\sum }_{k=1}^n{q}^k$，求法同理。
剩余还有一个$cn{t}_0^2$，过程类似，略。

// luogu-judger-enable-o2
#include<bits/stdc++.h>
#define rep(i,a,b) for(int i=a;i<=b;i++)
#define Rep(i,v) rep(i,0,(int)v.size()-1)
#define lint long long
#define mod 1000000007
#define ull unsigned lint
#define db long double
#define pb push_back
#define mp make_pair
#define fir first
#define sec second
#define gc getchar()
#define debug(x) cerr<<#x<<"="<<x
#define sp <<" "
#define ln <<endl
using namespace std;
typedef pair<int,int> pii;
typedef set<int>::iterator sit;
inline int inn()
{
    int x,ch;while((ch=gc)<'0'||ch>'9');
    x=ch^'0';while((ch=gc)>='0'&&ch<='9')
        x=(x<<1)+(x<<3)+(ch^'0');return x;
}
inline int fast_pow(int x,int k,int ans=1) { for(;k;k>>=1,x=(lint)x*x%mod) (k&1)?ans=(lint)ans*x%mod:0;return ans; }
inline int inv(int x) { return fast_pow(x,mod-2); }
inline lint squ(int x) { return (lint)x*x; }
inline int solve(int p,int c)
{
    int n=fast_pow(p,c),t=n-fast_pow(p,c-1);if(t<0) t+=mod;
    int q=inv(p),v=fast_pow(q,c),z=inv(q-1),c0=q*(v-1ll)%mod*z%mod;
    int c1=((lint)c*v%mod*q-c0)%mod*z%mod;if(c1<0) c1+=mod;
    int c2=((lint)c*c%mod*v%mod*q%mod-q-2ll*(c1-q)%mod+(lint)q*(c0-v)%mod)*z%mod;if(c2<0) c2+=mod;
    return ((lint)t*t%mod*(p-1)%mod*n%mod*c2+squ((lint)n*n%mod-t*(p-1ll)%mod*n%mod*c1%mod))%mod;
}
int main()
{
    int ans=1;
    for(int T=inn(),p,c;T;T--) p=inn(),c=inn(),ans=(lint)ans*solve(p,c)%mod;
    return !printf("%d\n",ans);
}