Problem

51nod

Solution

每次遇到这种神仙题，总有人告诉我这是出了很多遍套路题……
Orz 知乎大佬：张一钊的回答

通过 min-max 容斥我们可以知道：

$$lcm(f_S)=\prod_{T\subseteq S}\gcd(f_T)^{(-1)^{|T|-1}}=\prod_{T\subseteq S}f_{\gcd(f_T)}^{(-1)^{|T|-1}}$$

构造数列 g 使得 $f_n=\prod_{d|n} g_d$，则：

$$lcm(f_S)=\prod_{T\subseteq S}\biggl(\prod_{d|\gcd(f_T)}g_d \biggr)^{(-1)^{|T|-1}}$$

$$=\prod_{T\subseteq S} \prod_{d}g_d^{\sum_{T\subseteq S,d|\gcd(f_T)}{(-1)^{|T|-1}}}$$

我们单独考虑一下那个很长的指数，我们用中括号表示一个布尔表达式

$$\sum_{T\subseteq S,d|\gcd(f_T)}{(-1)^{|T|-1}}=[\exists_{x\in S} d|x]$$

为什么是这样的？首先如果 S 集合中存在一个数是 d 的倍数，那么我们必定能从 S 集合中找到一个极大子集合使得它的 gcd 是 d 的倍数，显然这个极大子集合是唯一的，考虑反证法，如果存在另一个等大的子集合，那么我们完全可以把这两个子集合合并起来，得到一个更大的子集合，它的 gcd 依然是 d 的倍数，与前面的假设相矛盾。由于我们删去一个数，只会使得集合的 gcd 乘上一个倍数，所以这个极大子集合的所有非空子集均合法。

假设极大子集合的大小为 l，那么就相当于对组合数进行了奇偶性分组，由于 i 从 1 开始所以：

$$\sum_{i=1}^l (-1)^{i-1}\binom l i=1$$

再看原来的式子，我们就可以得到

$$=\prod_{\exists d|x \in S}g(d)$$

怎么构造 g 呢，把之前我们定义 g 的式子改一下即有

$$g(n)=f(n)*\biggl(\prod_{d|n,d\not=n} g(d) \biggr)^{-1}$$

喵妙啊！

然后我们套路地改分解质因数为枚举倍数就可以做到 $O(n\ln n)$啦！

Code

#include <cstdio>
#define rg register
using namespace std;
typedef long long ll;
const int maxn=50010,maxm=1000010,mod=1e9+7;
template <typename Tp> inline int getmin(Tp &x,Tp y){return y<x?x=y,1:0;}
template <typename Tp> inline int getmax(Tp &x,Tp y){return y>x?x=y,1:0;}
template <typename Tp> inline void read(Tp &x)
{
    x=0;int f=0;char ch=getchar();
    while(ch!='-'&&(ch<'0'||ch>'9')) ch=getchar();
    if(ch=='-') f=1,ch=getchar();
    while(ch>='0'&&ch<='9') x=x*10+ch-'0',ch=getchar();
    if(f) x=-x;
}
int n,m,tmp,ans=1,f[maxm],g[maxm],vis[maxm];
inline int pls(int x,int y){return x+y>=mod?x+y-mod:x+y;}
int power(int x,int y)
{
    int res=1;
    for(;y;y>>=1,x=(ll)x*x%mod)
      if(y&1)
        res=(ll)res*x%mod;
    return res;
}
int main()
{
    #ifndef ONLINE_JUDGE
    freopen("in.txt","r",stdin);
    #endif
    read(n);
    for(rg int i=1;i<=n;i++){read(tmp);getmax(m,tmp);vis[tmp]=1;}
    f[1]=g[1]=1;
    for(rg int i=2;i<=m;i++) f[i]=g[i]=pls(f[i-1],f[i-2]);
    for(rg int i=1;i<=m;i++)
    {
        tmp=power(g[i],mod-2);
        for(rg int j=i<<1;j<=m;j+=i)
          g[j]=(ll)g[j]*tmp%mod;
    }
    for(rg int i=1;i<=m;i++)
      for(rg int j=i;j<=m;j+=i)
        if(vis[j])
        {ans=(ll)ans*g[i]%mod;break;}
    printf("%d\n",ans);
    return 0;
}