题意:求\(\sum_{i=1}^n\sum_{j=1}^mgcd(i,j)^k%1e9+7\)

题解:考虑枚举gcd,原式可化简为\(\sum_{d=1}^{n}d^k\sum_{i=1}^n\sum_{j=1}^m[gcd(i,j)==d]\)后面部分很明显是最基础的莫比乌斯反演,

那么有\(\sum_{d=1}^{n}d^k\sum_{x=1}^{\lfloor \frac{n}{d} \rfloor}\mu(x)*{\lfloor \frac{n}{d*x} \rfloor}*{\lfloor \frac{m}{d*x} \rfloor}\)

考虑枚举t=dx,(这里的套路很重要!!!*),那么有\(\sum_{t=1}^n{\lfloor \frac{n}{t} \rfloor}*{\lfloor \frac{m}{t} \rfloor}\sum_{d|t}d^k*\mu({\frac{t}{d}})\)

后面是一个积性函数可以O(n)预处理,前面可以分块,这里假设后面的积性函数是\(f(n)=\sum_{d|n}d^k*\mu({\frac{n}{d}})\)

\(f(n)=\prod_{i=1}^kf(p_i^{x_i})=\prod_{i=1}^k\mu(1)*p_i^{k*x_i}+\mu(p_i)*p_i^{k*(x_i-1)}=\prod_{i=1}^kp_i^{k*(x_i-1)}*(p_i^k-1)\)

这里第二步其他项没有的情况是\(\mu\)函数质因子两个以上就是0了,然后f(n)可以在线性筛的时候处理,可以先把素数的k次幂处理出来

/**************************************************************    Problem: 4407    User: walfy    Language: C++    Result: Accepted    Time:19252 ms    Memory:103832 kb****************************************************************/ //#pragma GCC optimize(2)//#pragma GCC optimize(3)//#pragma GCC optimize(4)//#pragma GCC optimize("unroll-loops")//#pragma comment(linker, "/stack:200000000")//#pragma GCC optimize("Ofast,no-stack-protector")//#pragma GCC target("sse,sse2,sse3,ssse3,sse4,popcnt,abm,mmx,avx,tune=native")#include
    
     #define fi first#define se second#define db double#define mp make_pair#define pb push_back#define pi acos(-1.0)#define ll long long#define vi vector
     
      #define mod 1000000007#define ld long double#define C 0.5772156649#define ls l,m,rt<<1#define rs m+1,r,rt<<1|1#define pll pair
      
       #define pil pair
       
        #define pli pair
        
         #define pii pair
         
          //#define cd complex
          
           #define ull unsigned long long#define base 1000000000000000000#define Max(a,b) ((a)>(b)?(a):(b))#define Min(a,b) ((a)<(b)?(a):(b))#define fin freopen("a.txt","r",stdin)#define fout freopen("a.txt","w",stdout)#define fio ios::sync_with_stdio(false);cin.tie(0)template
           
            inline T const& MAX(T const &a,T const &b){return a>b?a:b;}template
            
             inline T const& MIN(T const &a,T const &b){return a
             
              =mod)a-=mod;}inline void sub(ll &a,ll b){a-=b;if(a<0)a+=mod;}inline ll gcd(ll a,ll b){return b?gcd(b,a%b):a;}inline ll qp(ll a,ll b){ll ans=1;while(b){if(b&1)ans=ans*a%mod;a=a*a%mod,b>>=1;}return ans;}inline ll qp(ll a,ll b,ll c){ll ans=1;while(b){if(b&1)ans=ans*a%c;a=a*a%c,b>>=1;}return ans;} using namespace std; const double eps=1e-8;const ll INF=0x3f3f3f3f3f3f3f3f;const int N=5000000+10,maxn=400000+10,inf=0x3f3f3f3f; int prime[N],cnt,k;ll f[N],qk[N];bool mark[N];void init(){ f[1]=1; for(int i=2;i
              
               m)swap(n,m); ll ans=0; for(int i=1,j;i<=n;i=j+1) { j=min(n/(n/i),m/(m/i)); ll te=1ll*(f[j]-f[i-1])*(n/i)%mod*(m/i)%mod; te=(te+mod)%mod; add(ans,te); } printf("%lld\n",ans); } return 0;}/******************** ********************/