该题难在内循环是j，然而j又是分母，所以没办法快速当做整体求和
所以把分母j换做外循环


把每个j的里面分成n/j份（每份值同，（i/j）^j），然后累加前面系数i
我们把（i/j）^j记录到p[k]，k=i/j=1,2,3,4,5......
每次由于j变，()^j 相对上次加1就是*k，且每次p[]项个数增加（前系数求和增大）

由于最后一份里面p[]的个数与前几份不一样所以单独处理 
for i：1~n
    q[]=1
for j:1~n
    每份宽度limit=n/j,系数左=j，系数右=2j-1
    for i：1~limit循环每一份，最后一份单独处理
        p[i]=前次*i
        ans+=p[i]*前系数累加
        系数左边界+=j,系数右边界+=j
    //处理最后一份：
    p[limit]=p[limit]*i
    ans+=p[limit]*前系数累加
代码： 
#include <bits/stdc++.h>
using namespace std;
 
const long long N=3e6+16,mod=1e9+7;
long long  n,l,r,ans;
long long p[N];
 
int main(int argc, char** argv) {
    cin>>n;
    for(int i=1;i<=N;i++) p[i]=1;
    for(int j=1;j<=n;j++){
        long long limit=n/j;
        long long  l=j,r=j+j-1;
        for(int i=1;i<limit;i++){
            p[i]=p[i]*i%mod;
            ans=ans+(l+r)*j/2%mod*p[i]%mod;
            ans%=mod;
            l+=j,r+=j;
        }
        p[limit]=p[limit]*limit%mod;
        ans=ans+(l+n)*(n-l+1)/2%mod*p[limit]%mod;
        ans%=mod;
    }
    ans%=mod;
    cout<<ans<<endl;
    return 0;
}