hdu6607 min25筛+杜教筛+伯努利数求k次方前缀和

推导过程类似https://www.cnblogs.com/acjiumeng/p/9742073.html
前面部分min25筛,后面部分杜教筛,预处理min25筛需要伯努利数

//#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<bits/stdc++.h>
//#include <bits/extc++.h>
#define fi first
#define se second
#define db double
#define mp make_pair
#define pb push_back
#define mt make_tuple
//#define pi acos(-1.0)
#define ll long long
#define vi vector<int>
#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 sqr(x) ((x)*(x))
#define pll pair<ll,ll>
#define pil pair<int,ll>
#define pli pair<ll,int>
#define pii pair<int,int>
#define ull unsigned long long
#define bpc __builtin_popcount
#define base 1000000000000000000ll
#define fin freopen("a.txt","r",stdin)
#define fout freopen("a.txt","w",stdout)
#define fio ios::sync_with_stdio(false);cin.tie(0)
#define mr mt19937 rng(chrono::steady_clock::now().time_since_epoch().count())
inline ll gcd(ll a,ll b){return b?gcd(b,a%b):a;}
inline void sub(ll &a,ll b){a-=b;if(a<0)a+=mod;}
inline void add(ll &a,ll b){a+=b;if(a>=mod)a-=mod;}
template<typename T>inline T const& MAX(T const &a,T const &b){return a>b?a:b;}
template<typename T>inline T const& MIN(T const &a,T const &b){return a<b?a:b;}
inline ll mul(ll a,ll b,ll c){return (a*b-(ll)((ld)a*b/c)*c+c)%c;}
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=mul(ans,a,c);a=mul(a,a,c),b>>=1;}return ans;}

using namespace std;
//using namespace __gnu_pbds;

const ld pi=acos(-1);
const ull ba=233;
const db eps=1e-5;
const ll INF=0x3f3f3f3f3f3f3f3f;
const int N=200000+10,maxn=10000000+10,inf=0x3f3f3f3f;

bool mark[maxn];
int prime[maxn],cnt,K,phi[maxn],cnt1,f[maxn];
ll g[N],id[2][N],val[N],sum[N],up,n;
ll inv2=qp(2,mod-2),inv6=qp(6,mod-2),inv[111],c[111][111],b[111];
map<ll,ll>phii;
ll getb(int n)
{
    if(b[n]!=-1)return b[n];
    if(n==0)return b[0]=1;
    ll ans=0;
    for(int i=0;i<n;i++)
        add(ans,c[n+1][i]*getb(i)%mod);
    ans=-ans*inv[n+1]%mod;
    ans=(ans%mod+mod)%mod;
    return b[n]=ans;
}
ll get(ll n,ll k)
{
    ll ans=0;
    for(int i=1;i<=k+1;i++)
        add(ans,c[k+1][i]*b[k+1-i]%mod*qp((n+1)%mod,i)%mod);
    return ans*inv[k+1]%mod;
}
void pre()
{
    phi[1]=1;
    for(int i=2;i<maxn;i++)
    {
        if(!mark[i])prime[++cnt1]=i,phi[i]=i-1;
        for(int j=1;j<=cnt1&&i*prime[j]<maxn;j++)
        {
            mark[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<maxn;i++)
    {
        f[i]=1ll*i*i%mod*phi[i]%mod;
        f[i]+=f[i-1];
        if(f[i]>=mod)f[i]-=mod;
    }
    inv[1]=1;
    for(ll i=2;i<111;i++)
        inv[i]=1ll*(mod-mod/i)*inv[mod%i]%mod;
    for(int i=0;i<111;i++)
    {
        c[i][0]=c[i][i]=1;
        for(int j=1;j<i;j++)
            c[i][j]=(c[i-1][j]+c[i-1][j-1])%mod;
    }
    memset(b,-1,sizeof b);
    getb(110);
}
void init()
{
    up=sqrt(n);
    for(int i=1;prime[i]<=up;i++)sum[i]=(sum[i-1]+qp(prime[i],K+1))%mod,cnt=i;
    int m=0;
    for(ll i=1,j;i<=n;i=j+1)
    {
        j=n/(n/i);
        val[++m]=n/i;
        g[m]=get(n/i,K+1);//insert val
        sub(g[m],1ll);
        if(n/i<=up)id[0][n/i]=m;
        else id[1][i]=m;
    }
    for(int j=1;j<=cnt;j++)for(int i=1;i<=m&&1ll*prime[j]*prime[j]<=val[i];i++)
    {
        ll te=val[i]/prime[j];
        int k=(te<=up)?id[0][te]:id[1][n/te];
        sub(g[i],1ll*(sum[j]-sum[j-1]+mod)*(g[k]-sum[j-1]+mod)%mod);
    }
}
ll getf(ll n)
{
    if(n<maxn)return f[n];
    if(phii.find(n)!=phii.end())return phii[n];
    ll ans=n%mod*((n+1)%mod)%mod*inv2%mod;ans=ans*ans%mod;
    for(ll i=2,j;i<=n;i=j+1)
    {
        j=n/(n/i);
        ll tj=j%mod,ti=i%mod;
        ll te=tj*(tj+1)%mod*(2ll*tj+1)%mod*inv6%mod-(ti-1)*ti%mod*(2ll*ti-1)%mod*inv6%mod;
        te=(te%mod+mod)%mod;
        sub(ans,te*getf(n/i)%mod);
    }
    return phii[n]=ans;
}
int main()
{
    pre();
    int T;scanf("%d",&T);
    while(T--)
    {
        scanf("%lld%d",&n,&K);
        init();
        ll ans=0;
        for(ll i=1,j;i<=n;i=j+1)
        {
            j=n/(n/i);
            int k1=(j<=up)?id[0][j]:id[1][n/j];
            int k2=(i-1<=up)?id[0][i-1]:id[1][n/(i-1)];
//            printf("%lld %lld\n",(g[k1]-g[k2]+mod)%mod,getf(n/i));
            add(ans,getf(n/i)*(g[k1]-g[k2]+mod)%mod);
        }
        printf("%lld\n",ans);
    }
    return 0;
}
/********************

********************/

原文地址:https://www.cnblogs.com/acjiumeng/p/11372892.html

时间: 2024-11-08 04:19:50

hdu6607 min25筛+杜教筛+伯努利数求k次方前缀和的相关文章

奇怪的数学题(51nod1847)——min_25筛+杜教筛+第二类斯特林数

题面 51nod1847 解析   定义$f(i)$为$i$的次大因数,特别地$f(1)=0$  那么我们就是去求$$\sum_{i=1}^{n}\sum_{j=1}^{n}f^{m}(gcd(i, j))$$ 这种东西的套路就是枚举$gcd$然后用欧拉函数处理, \begin{align*}\sum_{i=1}^{n}\sum_{j=1}^{n}f^{m}(gcd(i, j)) &= \sum_{i=1}^{n}\sum_{j=1}^{\left \lfloor \frac{n}{i}\rig

【数论】狄利克雷卷积及其快速计算方法及杜教筛

目录(假的 狄利克雷卷积基础知识 数论函数 狄利克雷卷积定义 狄利克雷卷积性质 常用卷积 卷积计算方法 最暴力的暴力 稍好的暴力 优美的暴力 莫比乌斯反演(待填坑) 杜教筛 经典杜教筛 第二种杜教筛 第三种杜教筛 背景 本人即将去CTS&APIO2019,由于一些特殊原因,发现自己数论突然变得很菜. 就决定在去的前一天,翻出来以前的数论学习资料看一看.翻到了czgj的校内狄利克雷卷积课件,发现其中提到了的任意数列\(f(n)\)和\(g(n)\)的狄利克雷卷积\((f*g)(n)\)(从1到n,

数论入门——莫比乌斯函数,欧拉函数,狄利克雷卷积,线性筛,莫比乌斯反演,杜教筛

一个菜鸡对数论的一点点理解... 莫比乌斯函数 定义函数\(\mu(n)\)为: 当n有平方因子时,\(\mu(n)=0\). 当n没有平方因子时,\(\mu(n)=(-1)^{\omega(n)}\),\(\omega(n)\)表示n不同质因子的个数. 性质1: \(\sum_{d|n}\mu(d)=[n=1]\) 证明:我们把n分解质因数,则原式\(=(-1+1)^{\omega(n)}=0\). 因为对于不同的质因子,只有选和不选两种方案,这是一个组合数相加的形式,偶数加奇数减,根据二项式

【bzoj 4176】 Lucas的数论 莫比乌斯反演(杜教筛)

Description 去年的Lucas非常喜欢数论题,但是一年以后的Lucas却不那么喜欢了. 在整理以前的试题时,发现了这样一道题目“求Sigma(f(i)),其中1<=i<=N”,其中 表示i的约数个数.他现在长大了,题目也变难了. 求如下表达式的值: 一行一个整数ans,表示答案模1000000007的值. Sample Input 2 Sample Output 8 HINT 对于100%的数据n <= 10^9. 题解: 解锁新技能:杜教筛. 再复习一下: 若$F(n)=\s

杜教筛 学习总结

看了看唐老师的blog,照猫画虎的做了几道题目,感觉对杜教筛有些感觉了 但是稍微有一点难度的题目还是做不出来,放假的时候争取都A掉(挖坑ing) 这篇文章以后等我A掉那些题目之后再UPD上去就好啦 由于懒得去写怎么用编辑器写公式,所以公式就准备直接copy唐老师的啦 首先积性函数和完全积性函数什么的就不再多说了 列举常见的积性函数: 1.约数个数函数和约数个数和函数 2.欧拉函数phi 3.莫比乌斯函数mu 4.元函数e 其中e(n)=[n==1] 5.恒等函数I 其中I(n)=1 6.单位函数

【bzoj4176】Lucas的数论 莫比乌斯反演+杜教筛

题目描述 去年的Lucas非常喜欢数论题,但是一年以后的Lucas却不那么喜欢了. 在整理以前的试题时,发现了这样一道题目“求Sigma(f(i)),其中1<=i<=N”,其中 表示i的约数个数.他现在长大了,题目也变难了. 求如下表达式的值: 其中f(ij)表示ij的约数个数. 他发现答案有点大,只需要输出模1000000007的值. 输入 第一行一个整数n. 输出 一行一个整数ans,表示答案模1000000007的值. 样例输入 2 样例输出 8 题解 莫比乌斯反演+杜教筛 首先有个神奇

杜教筛进阶+洲阁筛讲解+SPOJ divcnt3

Part 1:杜教筛进阶在了解了杜教筛基本应用,如$\sum_{i=1}^n\varphi(i)$的求法后,我们看一些杜教筛较难的应用.求$\sum_{i=1}^n\varphi(i)*i$考虑把它与$id$函数狄利克雷卷积后的前缀和.$$\sum_{i=1}^n\sum_{d|i}\varphi(d)*d*\frac id=\sum_{i=1}^ni^2$$枚举$T=\frac id$,原式化为$$\sum_{T=1}^nT\sum_{d=1}^{\lfloor\frac nT\rfloor}

【51Nod 1237】最大公约数之和 V3 莫比乌斯反演+杜教筛

题意 求$\sum_{i=1}^{n}\sum_{j=1}^{n}(i,j)$ 枚举约数 $$ \begin{align} ans &=\sum_{d=1}^{n}\sum_{i=1}^{n}\sum_{j=1}^{n}[(i,j)=d] \ &=\sum_{d=1}^{n}\sum_{i=1}^{\lfloor \frac{n}{d} \rfloor}\sum_{j=1}^{\lfloor \frac{n}{d} \rfloor}[(i,j)=1] \ \end{align} $$ 利用

【Luogu3768】简单的数学题(莫比乌斯反演,杜教筛)

[Luogu3768]简单的数学题(莫比乌斯反演,杜教筛) 题面 洛谷 \[求\sum_{i=1}^n\sum_{j=1}^nijgcd(i,j)\] $ n<=10^9$ 题解 很明显的把\(gcd\)提出来 \[\sum_{d=1}^nd\sum_{i=1}^n\sum_{j=1}^nij[gcd(i,j)==d]\] 习惯性的提出来 \[\sum_{d=1}^nd^3\sum_{i=1}^{n/d}\sum_{j=1}^{n/d}ij[gcd(i,j)==1]\] 后面这玩意很明显的来一发