[反演] Project Euler 608. Divisor Sums

$D(m,n)=\displaystyle\sum_{d|m}\sum_{k=1}^n\sigma_{\small 0}(kd)$

求 $D(200!,10^{12})$

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

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

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

那么就可以推式子

$$D(m,n)=\sum_{d|m}\sum_{k=1}^n\sum_{a|d}\sum_{b|k}[(a,b)=1]$$

$$=\sum_{a|m}\sum_{b=1}^n[(a,b)=1]\sigma_{\small 0}({m\over a})\lfloor{n\over b}\rfloor$$

$$=\sum_{t|m,t\le n}\mu(t)\sum_{a|{m\over t}}\sigma_{\small 0}({m\over at})\sum_{b=1}^{\lfloor{n\over t}\rfloor}\lfloor{n\over bt}\rfloor$$

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

$\sum_{b=1}^{n}\lfloor{n\over b}\rfloor$ 分块算，记忆化一下复杂度大概是 $O(n^{2\over 3})$

$\sum_{a|n}\sigma_{\small 0}({n\over a})$ 就是 $n$ 的所有约数的约数个数的和

稍微推一下就可以知道这个东西等价于 $\prod{(c_i+1)(c_i+2)\over 2}$，$n=\prod p_i^{c_i}$

这个在dfs过程中能 $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;
}