洛谷3678：简单的数学题（画柿子+杜教筛）

最新推荐文章于 2022-11-11 12:42:19 发布

Fat_tu

最新推荐文章于 2022-11-11 12:42:19 发布

阅读量452

点赞数

分类专栏：数论杜教筛文章标签：数学

本文链接：https://blog.csdn.net/q582116859/article/details/78867907

版权

数论同时被 2 个专栏收录

19 篇文章 0 订阅

订阅专栏

杜教筛

3 篇文章 0 订阅

订阅专栏

最近看了《魔法少女伊莉雅~雪下的誓言》，
这部番除了百合与萌以外还有其他东西，
就和当初看《fate》的我一样。

题面
题意：求

\sum i = 1 n \sum j = 1 n i j (i, j)

$\sum_{i=1}^n\sum_{j=1}^nij(i,j)$

$n≤10^{10}$

这就是数论题，自己数学课上半个小时推出来，看Kscla的博客一个下午才懂。

根据我两页草稿纸的结果，原式为

\sum i = 1 n i 3 \sum j = 1 ⌊ n i ⌋ j 2 φ (j)

$\sum_{i=1}^ni^3\sum_{j=1}^{\lfloor\frac{n}{i}\rfloor}j^2\varphi(j)$
对右边进行杜教筛就可以了。
（我今天本来想学洲阁筛的，结果沉迷到这题了）

#include <iostream>
#include <fstream>
#include <algorithm>
#include <cmath>
#include <ctime>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <map>

using namespace std;
#define mmst(a, b) memset(a, b, sizeof(a))
#define mmcp(a, b) memcpy(a, b, sizeof(b))

typedef long long LL;

const int N=10100000,nn=10000000;

LL n,p,ans;
int phi[N],prime[N],num;
bool b[N];

map<LL,LL> mp;

LL get1(LL x)
{
    LL y=x+1;
    if(!(y&1))
    y=y/2;
    else
    x=x/2;
    x%=p;
    y%=p;
    return x*x%p*y%p*y%p;
}

LL get2(LL x)
{
    LL y=x+1,z=2*x+1;
    bool f2=false,f3=false;
    if ( !f2 && !(x%2) ) x/=2,f2=true;
    if ( !f3 && !(x%3) ) x/=3,f3=true;
    if ( !f2 && !(y%2) ) y/=2,f2=true;
    if ( !f3 && !(y%3) ) y/=3,f3=true;
    if ( !f2 && !(z%2) ) z/=2,f2=true;
    if ( !f3 && !(z%3) ) z/=3,f3=true;
    x%=p;
    y%=p;
    z%=p;
    return x*y%p*z%p;
}

LL work(LL x)
{
    if(x<=nn)
    return phi[x];
    if(mp.count(x))
    return mp[x];
    LL res=get1(x),last;
    for(LL i=2;i<=x;i=last+1)
    {
        last=x/(x/i);
        res=(res-(get2(last)-get2(i-1)+p)%p*work(x/i)%p+p)%p;
    }
    mp[x]=res;
    return res;
}

int main()
{
    cin>>p>>n;

    phi[1]=1;
    for(int i=2;i<=nn;i++)
    {
        if(!b[i])
        {
            phi[i]=i-1;
            prime[++num]=i;
        }
        for(int j=1;j<=num&&prime[j]*i<=nn;j++)
        {
            b[i*prime[j]]=1;
            if(i%prime[j]==0)
            {
                phi[i*prime[j]]=phi[i]*prime[j];
                break;
            }
            phi[i*prime[j]]=phi[i]*(prime[j]-1);
        }
    }

    for(int i=1;i<=nn;i++)
    phi[i]=((LL)phi[i-1]+(LL)phi[i]*i%p*i%p)%p;

    LL last;
    for(LL i=1;i<=n;i=last+1)
    {
        last=n/(n/i);
        ans=(ans+(get1(last)-get1(i-1)+p)%p*work(n/i)%p)%p;
    }

    cout<<ans<<endl;

    return 0;
}