Power oj 2810 Grisaia 杜教筛

2810: Grisaia

Time Limit: 12000 MS Memory Limit: 1048576 KB
Total Submit: 83 Accepted: 15 Page View: 131
Submit Status Discuss

×

Submit Problem 2810:Grisaia

×Sorry, youd don't have permission to submit solution.

 

SubmitCancel

Description

Kazami Kazuki is a talented student.

One day, she met a challengeable problem: calculate the value of

ans=n∑i=1i∑j=1(n mod(i×j))ans=∑i=1n∑j=1i(n mod(i×j))

She worked it out easily. Is it easy for you too?

Input

The first line contains an integer TT representing the number of test cases. In each test case, there is an integer nn in one line. • 1≤T≤51≤T≤5 • 1≤n≤10111≤n≤1011 • It is guaranteed there is at most one test case satisfying that n>109n>109 .

Output

For each test case, output the answer in one line.

• Sample Input
• Raw

2 3 7

2 3 7

• Sample Output
• Raw

10 145

10 145

Source

The 2018 Sichuan Provincial Collegiate Programming Contest

杜教筛神仙题

先推一波公式

把n提出去

只看后面的式子

枚举  i*j=k

后半部分、

可以看成 k的小的因子

我们暂时把他当成

一会儿再稍作处理

式子就变成了

显然前面的可以数论分块处理

然而后面的要 好好处理 

设  

我们就努力的去寻找的另一半

使得

两函数的狄利克雷卷积容易计算

再根据杜教筛公式

有

问题又来了我们不会求

的前缀和

再次狄利克雷卷积杜教筛

设

有

然后我们就将求了出来

因为我们要求

对于一些完全平方数  要加上  k  （完全平方数有奇数个因子   ）

 

也就是加上  1到  刚好比n小的平方数  的平方和

然后再除以二

然后就没了

全程数论分块

部分值开__int128即可AC

10S飘过

#include<bits/stdc++.h>
const int maxn=4e7+10;
using namespace std;
typedef long long ll;
typedef __int128 lll;
bool vis[maxn];
int mu[maxn];
ll sum_muii[maxn];
ll d[maxn];
int a[maxn];
int cnt,prim[maxn];
unordered_map<ll,lll> w1;
unordered_map<ll,lll> w2;
inline void print(lll x)
{
    if(x<0)
    {
        putchar('-');
        x=-x;
    }
    if(x>9)
        print(x/10);
    putchar(x%10+'0');
}
void init()
{
    mu[1]=1;
    d[1]=1;
    for(ll i=2;i<maxn;i++)
    {
        if(!vis[i])
        {
            prim[++cnt]=i;
            d[i]=2*i;
            a[i]=1;
            mu[i]=-1;
        }
        for(int j=1;j<=cnt&&prim[j]*i<maxn;j++)
        {
            vis[i*prim[j]]=1;
            if(i%prim[j]==0)
            {
                d[i*prim[j]]=d[i]/(a[i]+1)*(a[i]+2)*prim[j];
                a[i*prim[j]]=a[i]+1;
                break;
            }
            else
            {
                d[i*prim[j]]=d[i]*d[prim[j]];
                a[i*prim[j]]=1;
                mu[i*prim[j]]=-mu[i];
            }
        }
    }
    for(ll i=1;i<maxn;i++)
    {
        sum_muii[i]=sum_muii[i-1]+mu[i]*i;
    }
    for(ll i=1;i<maxn;i++)
    {
        d[i]=d[i-1]+d[i];
    }
}
inline lll djsmuii(ll x)//mu[i]*i 筛
{
    if(x<maxn)
        return sum_muii[x];
    if(w1[x])
        return w1[x];
    lll ans=1;
    for(ll l=2,r;l<=x;l=r+1)
    {
        r=x/(x/l);
        ans-=(lll)(r+l)*(r-l+1)/2*djsmuii(x/l);
    }
    w1[x]=ans;
    return ans;
}
inline lll djsknn(ll x) //k*k的因子数目筛
{
    if(x<maxn)
        return d[x];
    if(w2[x])
        return w2[x];
    lll ans=(lll)x*(x+1)/2;
    for(ll l=2,r;l<=x;l=r+1)
    {
        r=x/(x/l);
        ans-=djsknn(x/l)*(djsmuii(r)-djsmuii(l-1));
    }
    w2[x]=ans;
    return ans;
}
inline lll ask(ll x)
{

    lll nth=sqrt(x+0.9);
    lll tp=nth*(nth+1)*(2*nth+1)/6;
    return (djsknn(x)+tp)/2;
}
lll solve(ll x)
{
    lll ans=(lll)(1+x)*x*x/2;
    for(ll l=1,r;l<=x;l=r+1)
    {
        r=x/(x/l);
        ans-=(ask(r)-ask(l-1))*(x/l);
    }
    return ans;
}
int main()
{
    init();
    int t;
    cin>>t;
    while(t--)
    {
        //        __int128 n;
        ll n;
        scanf("%lld",&n);
        //n=read();
        //cin>>n;
        //        printf("%lld\n",solve(n));
        print(solve(n));
        puts("");
    }
    return 0;
}