[反演] Project Euler 608. Divisor Sums

D(m,n)=d|mk=1nσ0(kd) D ( m , n ) = ∑ d | m ∑ k = 1 n σ 0 ( k d )

D(200!,1012) D ( 200 ! , 10 12 )

做的第一道PE题,我终于会反演啦

kd k d 的每个约数可以表示成 k k 的约数乘上 d 的约数,但是这样一个约数会被表示很多次

假设 a|k a | k b|d b | d ,那么如果 gcd(a,b)=1 gcd ( a , b ) = 1 就唯一表示了 a×b a × b

那么就可以推式子

D(m,n)=d|mk=1na|db|k[(a,b)=1] D ( m , n ) = ∑ d | m ∑ k = 1 n ∑ a | d ∑ b | k [ ( a , b ) = 1 ]

=a|mb=1n[(a,b)=1]σ0(ma)nb = ∑ a | m ∑ b = 1 n [ ( a , b ) = 1 ] σ 0 ( m a ) ⌊ n b ⌋

=t|m,tnμ(t)a|mtσ0(mat)b=1ntnbt = ∑ t | m , t ≤ n μ ( t ) ∑ a | m t σ 0 ( m a t ) ∑ b = 1 ⌊ n t ⌋ ⌊ n b t ⌋

因为 μ(t)0 μ ( t ) ≠ 0 对答案才有贡献,所以只要枚举 m m 的素数的子集,使得这些素数乘积小于等于 n 就行了,200以内的素数只有46个,而第二个限制就使得有效的情况数不会很多

nb=1nb ∑ b = 1 n ⌊ n b ⌋ 分块算,记忆化一下复杂度大概是 O(n23) O ( n 2 3 )

a|nσ0(na) ∑ a | n σ 0 ( n a ) 就是 n n 的所有约数的约数个数的和

稍微推一下就可以知道这个东西等价于 (ci+1)(ci+2)2 n=pcii n = ∏ p i c i

这个在dfs过程中能 O(1) O ( 1 ) 处理

跑个20s就跑出来啦

#include <cstdio>
#include <iostream>
#include <algorithm>
#include <map>
#include <ctime>

using namespace std;

typedef long long ll;

const int P=1e9+7,inv2=P+1>>1;

int p[210],mu[210],inv[310];

inline int Pow(int x,int y){
  int ret=1;
  for(;y;y>>=1,x=1LL*x*x%P) if(y&1) ret=1LL*ret*x%P;
  return ret;
}

inline void Pre(const int &n){
  mu[1]=1;
  for(int i=2;i<=n;i++){
    if(!p[i]) p[++*p]=i,mu[i]=-1;
    for(int j=1;j<=*p && p[j]*i<=n;j++){
      p[p[j]*i]=1;
      if(i%p[j]) mu[i*p[j]]=-mu[i];
      else break;
    }
  }
  for(int i=1;i<=300;i++) inv[i]=Pow(i,P-2);
}

int c[210],prod=1;

ll n;
int m,ans;

inline int f(int x){
  return 1LL*(x+1)*(x+2)/2%P;
}

inline int invf(int x){
  return 2LL*inv[x+1]*inv[x+2]%P;
}

map<ll,int> M;

inline int g(ll n){
  if(M.count(n)) return M[n];
  int ret=0;
  for(ll i=1,j;i<=n;i=j+1){
    j=n/(n/i);
    int a=(j-i+1)%P,b=(n/i)%P;
    ret=(ret+1LL*a*b)%P;
  }
  return M[n]=ret;
}

ll tot;

void dfs(int x,int num,ll cc){
  if(cc>n) return ;
  if(x>*p){
    ans=(ans+1LL*num*prod*g(n/cc))%P; return ;
  }
  dfs(x+1,num,cc);
  prod=1LL*prod*invf(c[x])%P*f(c[x]-1)%P; c[x]--;
  dfs(x+1,-num,cc*p[x]);
  prod=1LL*prod*invf(c[x])%P*f(c[x]+1)%P; c[x]++;
}

int main(){
  cin>>m>>n; Pre(m);
  for(int i=1;i<=*p;i++){
    int cur=m;
    while(cur/p[i]) c[i]+=(cur/=p[i]);
    prod=1LL*prod*f(c[i])%P;
  }
  dfs(1,1,1);
  printf("%d\n",(ans+P)%P);
  printf("%d\n",clock());
  return 0;
}
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值