[数论 反演] HDU6053. TrickGCD

可以先容斥一下，用所有方案减去gcd为1的方案
那么只要能求gcd为1的方案就可以了
推一下式子

$$\sum_{i_1=1}^{a_1}\sum_{i_2=1}^{a_2}\sum_{i_3=1}^{a_3}…\sum_{i_n=1}^{a_n}[ \gcd(i_1,i_2,i_3,..,i_n) =1]$$

$$=\sum_{i_1=1}^{a_1}\sum_{i_2=1}^{a_2}\sum_{i_3=1}^{a_3}…\sum_{i_n=1}^{a_n}\sum_{t|(i_1,i_2,...,i_n)}\mu(t)$$

$$=\sum_{t=1}^{a_{min}}\mu(t)\prod_{i=1}^{n}\lfloor{a_i\over t}\rfloor$$

但是这样还不够，继续推

$$=\sum_{t=1}^{a_{min}}\mu(t)\prod_{i=1}^{\lfloor{a_{max}\over t}\rfloor}i^{num(it,(i+1)t-1)}$$

$num(l,r)$表示有多上数大于等于l，小于等于r

这样就$O(n\log^2n)$了

#include <cstdio>
#include <iostream>
#include <algorithm>
#include <vector>
#include <cstring>
#include <string>
#include <map>
#include <set>
#include <cmath>
#define fi first
#define se second
#define pb push_back

using namespace std;

typedef long long ll;
typedef unsigned int U;
typedef unsigned long long Ull;
typedef pair<int,int> ii;

const int N=100010,P=1e9+7;

int mu[N],f[N],p[N];

inline char nc(){
  static char buf[100000],*p1=buf,*p2=buf;
  return p1==p2&&(p2=(p1=buf)+fread(buf,1,100000,stdin),p1==p2)?EOF:*p1++;
}

inline void rea(int &x){
  char c=nc(); x=0; int f=1;
  for(;c>'9'||c<'0';f=c=='-'?-f:f,c=nc());for(;c>='0'&&c<='9';x=x*10+f*(c-'0'),c=nc());
}

inline void add(int &x,int y){
  (x+=(P+y)%P)%=P;
}

int t,n,a[N],lg2[N],c[N];
int b[N];

inline void first(const int x){
  mu[1]=1;
  ll i,j;
  for(i=2;i<=x;i++){
    if(!p[i]) p[++p[0]]=i,mu[i]=-1;
    for(j=1;i*p[j]<=x;j++)
      if(p[p[j]*i]=1,i%p[j]) mu[i*p[j]]=-mu[i];
      else { mu[i*p[j]]=0; break;}
  }
  for(int i=1;i<=x;i++)
    lg2[i]=lg2[i-1]+((1<<lg2[i-1]+1)==i);
  for(int i=1;i<=x;i++)
    f[i]=(f[i-1]+P+mu[i])%P;
}

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;
}

int main(){
  freopen("1.in","r",stdin);
  freopen("1.out","w",stdout);
  rea(t); first(100000); int tt=0;
  while(t--){
    rea(n); int mn=1<<30,ans=1,mx=1;
    memset(b,0,sizeof(b));
    for(int i=1;i<=n;i++)
      rea(a[i]),mn=min(mn,a[i]),ans=1LL*ans*a[i]%P,mx=max(mx,a[i]),b[a[i]]++;
    for(int i=1;i<=mx;i++) b[i]+=b[i-1];
    int t=lg2[n];
    for(int i=1;i<=mn;i++){
      int cur=(P+mu[i])%P;
      for(int j=1;i*j<=mx;j++)
    cur=1LL*cur*Pow(j,b[min(i*(j+1)-1,mx)]-b[i*j-1])%P;
      add(ans,-cur);
    }
    printf("Case #%d: %d\n",++tt,ans);
  }
  return 0;
}