zoukankan      html  css  js  c++  java
  • 「笔记」莫比乌斯反演

    Updated on 2020.8.6
    巨幅更新,对积性函数和狄利克雷卷积部分进行重构。
    新增对一类特殊求和式的讲解。

    Updated on 2020.9.9
    添加了几个例题。


    前置知识

    小碎骨

    艾佛森括号 ([P] = egin{cases} 1 & ext{If P is true}\ 0 & ext{Otherwise} end{cases})
    此处 (P) 是一个可真可假的命题。

    引理1

    [forall a,b,cin mathbb{Z},leftlfloordfrac{a}{bc} ight floor = leftlfloor{dfrac{leftlfloordfrac{a}{b} ight floor}{c}} ight floor ]

    证明

    [dfrac{a}{b} = leftlfloor{dfrac{a}{b}} ight floor + r(0le r < 1) ]

    [leftlfloordfrac{a}{bc} ight floor = leftlfloordfrac{a}{b} imesdfrac{1}{c} ight floor = leftlfloor{dfrac{1}{c} imesleft({leftlfloor{dfrac{a}{b} + r} ight floor} ight)} ight floor = leftlfloor{dfrac{leftlfloordfrac{a}{b} ight floor}{c} + dfrac{r}{c}} ight floor = leftlfloor{dfrac{leftlfloor{dfrac{a}{b}} ight floor}{c}} ight floor ]

    数论分块

    内容独立了出来,详细内容见 数论分块 - Luckyblock

    对于一类含有(leftlfloorfrac{n}{i} ight floor)的求和式 ((n) 为常数),由于(leftlfloorfrac{n}{i} ight floor)单调不增,故存在多个区间([l,r]), 使得(leftlfloorfrac{n}{i} ight floor = leftlfloorfrac{n}{j} ight floor(i,jin [l,r])) 成立。

    对于任意一个(i),最大的满足上式的 (j=leftlfloor{dfrac{n}{leftlfloor{dfrac{n}{i}} ight floor}} ight floor)


    积性函数

    定义

    (gcd(x,y) = 1)(f(xy)=f(x)f(y)),则 (f(n)) 为积性函数。

    性质

    (f(x))(g(x))均为积性函数,则以下函数也为积性函数:

    [egin{aligned} & h(x) = f(x^p)\ & h(x) = f^p(x)\ & h(x) = f(x)g(x)\ & h(x) = sum_{dmid x} f(d)gleft(dfrac{x}{d} ight) end{aligned}]

    常见积性函数

    • 单位函数 (e(n) = [n = 1])
    • 幂函数 (operatorname{Id}_{k}(n) = n^k)(operatorname{id}_1(n)) 通常简记为(operatorname{id}(n))
    • 常数函数 (1(n) = 1)
    • 因数个数 (operatorname{d}(n) = sumlimits_{dmid n} 1)
    • 除数函数 (sigma_{k}(n) = sumlimits_{dmid n} d^k)
      (k=1) 时为因数和函数,通常简记为 (sigma(n))
      (k=0) 时为因数个数函数 (sigma_{0}(n))
    • 欧拉函数 (varphi(n) = sumlimits_{i=1}^{n} [gcd(i,n) = 1])
    • 莫比乌斯函数 (mu(n) = egin{cases}1 &n=1\0 &n ext{含有平方因子}\(-1)^k &k ext{为} n ext{的本质不同质因子个数} end{cases})

    不是很懂上面写的什么玩意?
    不用深究,有个印象继续往下看就好。


    莫比乌斯函数

    定义

    (mu) 为莫比乌斯函数,定义为

    [mu(n) = egin{cases} 1 &n=1\0 &n ext{含有平方因子}\(-1)^k &k ext{为} n ext{的本质不同质因子个数} end{cases}]

    解释

    (n = prodlimits_{i=1}^{k} p_{i}^{c_i}),其中(p_i)为质因子,(c_ige 1)

    1. (n=1)时,(mu (n) = 1)
    2. (n ot ={1})时 ,
      • (exist iin [1,k], c_i > 1) 时,(mu (n) = 0)
        当某质因子出现次数大于(1)时,(mu (n) = 0)
      • (forall iin [1,k], c_i = 1) 时,(mu (n) = (-1)^k)
        当每个质因子只出现一次时,即(n = prodlimits_{i=1}^{k}p_i)({p_i})中元素唯一。
        (mu (n) = (-1)^k),此处(k)为质因子的种类数。

    性质

    莫比乌斯函数是积性函数,且具有以下性质

    [large sum_{dmid n} mu (d) = [n=1] ]

    证明,设 (n = prodlimits_{i=1}^{k}{p_i^{c_i}}, n' = prodlimits_{i=1}^{k}{p_i})

    • 根据莫比乌斯函数定义,则有:(sumlimits_{dmid n}{mu(d)} = sumlimits_{dmid n'}{mu(d)})
    • (n') 的某因子 (d),有 (mu (d) = (-1)^i),则它由 (i) 个 本质不同的质因子组成。
      由于质因子总数为 (k),满足上式的因子数为 (C_{k}^{i})
    • 对于原求和式,转为枚举 (mu(d)) 的值。
      (sumlimits_{dmid n'}{mu(d)} = sumlimits_{i=0}^{k}{C_{k}^{i} imes (-1)^i} = sumlimits_{i=0}^{k}{C_{k}^{i} imes (-1)^i imes 1^{k-i}})
      根据二项式定理,上式 (= (1+(-1))^k)
      易知该式在 (k=0),即 (n=0) 时为 (1),否则为 (0)

    反演常用结论

    一个反演常用结论:

    [[gcd(i,j) = 1] = sumlimits_{dmid gcd(i,j)} {mu (d)} ]

    证明 1:
    (n = gcd(i,j)),则右(= sumlimits_{dmid n} {mu (d)} = [n = 1] = [gcd(i,j) = 1]=) 左。

    证明 2:
    暴力展开:([gcd(i,j) = 1] = e(gcd(i,j)) = sumlimits_{dmid gcd(i,j)}mu(d))

    线性筛求莫比乌斯函数

    (mu) 为积性函数,因此可以线性筛莫比乌斯函数。

    int pnum, mu[kMaxn], p[kMaxn];
    bool vis[kMaxn];
    
    void Euler(int lim_) {
      vis[1] = true, mu[1] = 1ll;
      for (int i = 2; i <= lim_; ++ i) {
        if (! vis[i]) {
          p[++ pnum] = i;
          mu[i] = - 1;
        }
        for (int j = 1; j <= pnum && i * p[j] <= lim_; ++ j) {
          vis[i * p[j]] = true;
          if (i % p[j] == 0) { //勿忘平方因子
            mu[i * p[j]] = 0;
            break;
          }
          mu[i * p[j]] = - mu[i];
        }
      }
    }
    

    狄利克雷(Dirichlet)卷积

    建议阅读 算法学习笔记(35): 狄利克雷卷积 By: Pecco

    定义两个数论函数 (f,g) 的狄利克雷卷积为

    [large(fast g) (n) = sum_{dmid n} f(d)gleft(dfrac{n}{d} ight) ]

    性质

    1. 显然满足 交换律,结合律,分配律。
      • (f ast g = g ast f)
      • ((f ast g) ast h = fast (gast h))
      • (fast (g+h) = fast g + fast h)
    2. (e) 为狄利克雷卷积的单位元,有((fast e)(n) = f(n))
    3. (f,g) 为积性函数,则 (fast g) 为积性函数。

    关于单位元 (e)

    有:

    [e = mu ast 1=sumlimits_{dmid n} mu (d) ]

    证明

    [egin{aligned} (fast e)(n) = sum_{dmid n} f(d)e(dfrac{n}{d}) = sum_{dmid n} f(d)left[dfrac{n}{d} = 1 ight] end{aligned}]

    • 对于([dfrac{n}{d} = 1]),当且仅当 (dfrac{n}{d}=1),即 (d=n) 时为 (1),否则为(0)
    • 则当 (d=n) 时,(f(d)left[dfrac{n}{d} = 1 ight] = f(n))
      (d ot ={n}) 时,(f(d)left[dfrac{n}{d} = 1 ight] = 0)

    综上,((fast e)(n) = sumlimits_{dmid n} f(d)left[dfrac{n}{d} = 1 ight] = f(n)),满足单位元的性质。
    (e = mu ast 1) 成立。

    除数函数与幂函数

    幂函数 (operatorname{Id}_{k}(n) = n^k)
    除数函数 (sigma_{k}(n) = sumlimits_{dmid n} d^k)

    显然有:

    [(operatorname{Id}_kast 1)(n) = sum_{dmid n} operatorname{Id_k}(d) = sum_{dmid n} d^k = sigma_k(n) ]

    (k=0) 时,(operatorname{Id_0}=1)(sigma_0) 为因数个数函数,有:

    [(1ast 1)(n) = sum_{dmid n}1 = sigma_0(n) ]

    欧拉函数与恒等函数

    [egin{aligned} varphi ast 1 =& operatorname{Id}\ varphi =& operatorname{Id}ast mu end{aligned}]

    对于一式,当 (n=p^m) 时((p) 为质数),有:

    [(varphi ast 1)(p^m) = sum_{dmid n}varphi(d) = varphi(1) +sum_{i=1}^{m}varphi(p^i) = 1 +sum_{i=1}^{m}(p^i-p^{i-1}) = p^m ]

    (p^i) 的因子有 (p^{i-1}) 个,为 (1sim p^{i-1}),故 (varphi(p^i) = p^i-p^{i-1})

    ((varphi ast 1)(n)) 为积性函数,则对于任意正整数 (n),有:

    [(varphi ast 1)(n) = (varphi ast 1)left(prod p^m ight) = prodleft(varphi ast 1 ight)(p^m) = prod p^m = n ]

    得证。

    对于 2 式,在 1 式基础上两侧同时 (ast mu) 即得。
    左侧变为 (varphi ast 1ast mu = varphi ast e = varphi)

    计算

    [(fast g) (n) = sum_{dmid n} f(d)gleft(dfrac{n}{d} ight) ]

    考虑枚举 (n) 的因子,将贡献累加即可。
    显然可以使用埃氏筛筛出所有前缀狄利克雷卷积,复杂度 (O(nklog n)),其中 (k) 是计算函数值的复杂度。


    莫比乌斯反演

    反演是个啥?反演

    公式

    (f(n),g(n))为两个数论函数。
    如果有

    [large f(n) = (gast 1)(n) = sumlimits_{dmid n}{g(d)} ]

    那么有

    [large g(n) = (mu ast f)(n)=sumlimits_{dmid n} {mu(d)fleft(dfrac{n}{d} ight)} ]

    证明

    方法一:运用卷积。

    原问题为:已知 (f = gast 1),证明 (g = fast mu)
    易知如下转化:(fast mu = gast 1 ast mu Longrightarrow fast mu = gast e = g)

    方法二:对原式进行数论变换。

    1. (sumlimits_{dmid n}g(d)) 替换(fleft(dfrac{n}{d} ight))

      [sum_{dmid n}{mu(d)sum_{kmid frac{n}{d}}g(k)} ]

    2. 变换求和顺序。

      [sum_{kmid n}g(k)sum_{dmid frac{n}{k}}{mu(d)} ]

    3. 依据 (sumlimits_{dmid n}{mu(d)} = [n=1]),仅当 (dfrac{n}{k} = 1) 时,(sumlimits_{dmid frac{n}{k}}{mu(d)} = 1),否则为 (0)
      此时(k=n),故上式等价于 (sumlimits_{kmid n} {[n=k]cdot g(k)} = g(n))

    举例

    狄利克雷(Dirichlet)卷积 部分可以知道一些积性函数的关系。
    尝试对它们进行反演:

    [e = mu ast 1iffmu = muast e = mu ]

    [sigma_k = operatorname{Id}_kast 1iff operatorname{Id}_k = mu ast sigma_k ]

    [operatorname{Id}=varphi ast 1iff varphi = operatorname{Id}ast mu ]


    关于一类求和式

    [sum_{i=1}^{n}sum_{j=1}^{m}f(gcd(i,j)) ]

    一般套路

    考虑构造出函数 (g),满足下列形式:

    [f(n) = gast 1 = sum_{dmid n}g(d) ]

    则求和式变为:

    [sum_{i=1}^{n}sum_{j=1}^{m}sum_{dmid gcd(i,j)}g(d) ]

    考虑算术基本定理,发现若满足 (dmid gcd (i,j)),则 (dmid i)(dmid j) 成立。
    考虑 (g(d)) 在何时做出贡献,调整枚举顺序:

    [sum_{d=1}g(d)sum_{i=1}^n[dmid i]sum_{j=1}^m[dmid j] ]

    (sumlimits_{i=1}^{n}[dmid i]) 等价于 (1sim n)(d) 的倍数的个数,则上式等价于:

    [sum_{d=1}g(d)leftlfloordfrac{n}{d} ight floorleftlfloordfrac{m}{d} ight floor ]

    数论分块求解即可。

    例 1

    [sum_{i=1}^{n}sum_{j=1}^{m}gcd(i,j) ]

    发现此题的 (f) 等价于 (operatorname{Id}),则上式等价于:

    [egin{aligned} &sum_{i=1}^{n}sum_{j=1}^{m}operatorname{Id}(gcd(i,j))\ =& sum_{i=1}^{n}sum_{j=1}^{m}sum_{dmid gcd(i,j)}varphi(d)\ =& sum_{d=1}varphi(d)leftlfloordfrac{n}{d} ight floorleftlfloordfrac{m}{d} ight floor end{aligned}]

    例 2

    [egin{aligned} &sumlimits_{i=1}^{n}sumlimits_{j=1}^{m}[gcd(i,j) = 1]\ =& sumlimits_{i=1}^{n}sumlimits_{j=1}^{m}e(gcd(i,j))\ =& sumlimits_{i=1}^{n}sumlimits_{j=1}^{m}sum_{dmid gcd(i,j)}mu (d)\ =& sum_{d=1}mu(d)leftlfloordfrac{n}{d} ight floorleftlfloordfrac{m}{d} ight floor end{aligned}]

    感悟

    卷点什么东西,把 (g) 卷出来。
    (g) 不一定是特殊意义的函数。


    例题

    [HAOI2011] Problem b

    (n) 组询问,每次给定参数 (a,b,c,d,k),求:

    [sumlimits_{x=a}^{b}sumlimits_{y=c}^{d}[gcd(x,y) = k] ]

    (1le n,k,a,b,c,dle 5 imes 10^4)(ale b,cle d)

    (f(n,m) = sumlimits_{i=1}^{n}sumlimits_{j=1}^{m}[gcd(i,j) = k])
    根据容斥原理,则原式等价于 (f(b,d) - f(a-1,d) - f(b,d-1) + f(a-1,d-1))
    (f) 变成了上述一类求和式的形式,考虑化简 (f)

    易知原式等价于

    [sumlimits_{i=1}^{leftlfloorfrac{n}{k} ight floor}sumlimits_{j=1}^{leftlfloorfrac{m}{k} ight floor}[gcd(i,j) = 1] ]

    代入反演常用结论 ([gcd(i,j) = 1] = sumlimits_{dmid gcd(i,j)} {mu (d)}),上式化为

    [sumlimits_{i=1}^{leftlfloorfrac{n}{k} ight floor}sumlimits_{j=1}^{leftlfloorfrac{m}{k} ight floor}sum_{dmid gcd(i,j)}{mu(d)} ]

    变换求和顺序,先枚举(dmid gcd(i,j)),可得

    [sum_{d=1}mu(d)sum_{i=1}^{leftlfloorfrac{n}{k} ight floor}[dmid i]sum_{j=1}^{leftlfloorfrac{m}{k} ight floor}[dmid j] ]

    对于上式后半的解释:当(dmid i)(dmid j) 时,(dmid gcd(i,j))

    易知(1sim leftlfloordfrac{n}{k} ight floor)(d) 的倍数有 (leftlfloordfrac{leftlfloordfrac{n}{k} ight floor}{d} ight floor = leftlfloordfrac{n}{kd} ight floor) 个(由引理 1 可知),原式变为

    [sum_{d=1}mu(d)leftlfloordfrac{n}{kd} ight floorleftlfloordfrac{m}{kd} ight floor ]

    预处理 (mu) 后,显然可以数论分块求解,复杂度(O(n + Tsqrt{n}))


    代码

    //知识点:莫比乌斯反演
    /*
    //By:Luckyblock
    */
    #include <cstdio>
    #include <cctype>
    #include <algorithm>
    #define ll long long
    const int MARX = 6e4 + 10;
    //=============================================================
    int N, a, b, c, d, k;
    int cnt, Mobius[MARX], Prime[MARX], Sum[MARX];
    bool vis[MARX];
    //=============================================================
    inline int read()
    {
        int f = 1, w = 0; char ch = getchar();
        for(; !isdigit(ch); ch = getchar()) if(ch == '-') f = -1;
        for(; isdigit(ch); ch = getchar()) w = (w << 3) + (w << 1) + (ch ^ '0');
        return f * w;
    }
    void GetMobius()
    {
        Mobius[1] = 1;
        int MAX = MARX - 10;
        for(int i = 2; i <= MAX; i ++)
        {
          if(! vis[i]) Mobius[i] = - 1, Prime[++ cnt] = i;
          for(int j = 1; j <= cnt && i * Prime[j] <= MAX; j ++)
          {
            vis[i * Prime[j]] = true;
            if(i % Prime[j] == 0) break;
            Mobius[i * Prime[j]] = - Mobius[i];
          }
        }
        for(int i = 1; i <= MAX; i ++) Sum[i] = Sum[i - 1] + Mobius[i];
    }
    ll Calc(int x, int y)
    {
        ll ans = 0ll; int max = std :: min(x, y);
        for(int l = 1, r; l <= max; l = r + 1)
          r = std :: min(x / (x / l), y / (y / l)),
          ans += (1ll * x / (1ll * l * k)) * (1ll * y / (1ll * l * k)) * (Sum[r] - Sum[l - 1]);
        return ans;
    }
    ll Solve()
    {
        a = read(), b = read(), c = read(), d = read(), k = read();
        return Calc(b, d) - Calc(b, c - 1) - Calc(a - 1, d) + Calc(a - 1, c - 1);
    }
    //=============================================================
    int main()
    { 
        N = read(); GetMobius();
        while(N --) printf("%lld
    ", Solve());
        return 0;
    }
    

    [国家集训队]Crash的数字表格

    给定 (n,m) , 求:

    [sum_{i=1}^nsum_{j=1}^{m} operatorname{lcm}(i,j)mod 20101009 ]

    (1le n,mle 10^7)

    易知原式等价于:

    [sum_{i=1}^{n}sum_{j=1}^{m}dfrac{ij}{gcd (i,j)} ]

    考虑枚举 (gcd(i,j)),设枚举量为 (d)
    (d=gcd(i,j)) 的充要条件是满足 (d|i, d|j)(gcd(dfrac{i}{d},dfrac{j}{d}) = 1),则原式等价于:

    [sum_{i=1}^nsum_{j=1}^m sum_{d=1} dfrac{ij}{d}[d|i][d|j][gcd(frac{i}{d}, frac{j}{d})=1] ]

    先枚举 (d),则原式等价于:

    [sum_{d=1}sum_{i=1}^{n}[dmid i]sum_{j=1}^m [dmid j][gcd(dfrac{i}{d}, dfrac{j}{d}=1)] dfrac{ij}{d} ]

    这个 (d) 很烦人,把 (i,j) 中的 (d) 提出来,变为枚举 (frac{i}{d})(frac{j}{d})
    消去 (dmid i)(dmid j) 的限定条件,则原式等价于:

    [egin{aligned} &sum_{d=1}sum_{i=1}^{leftlfloorfrac{n}{d} ight floor}sum_{j=1}^{leftlfloorfrac{m}{d} ight floor}[gcd(i,j)=1]dfrac{id imes jd}{d}\ =& sum_{d=1}sum_{i=1}^{leftlfloorfrac{n}{d} ight floor}sum_{j=1}^{leftlfloorfrac{m}{d} ight floor}[gcd(i,j)=1]ijd\ =& sum_{d=1}d sum_{i=1}^{leftlfloorfrac{n}{d} ight floor}sum_{j=1}^{leftlfloorfrac{m}{d} ight floor}[gcd(i,j)=1]ij end{aligned}]

    单独考虑后半部分,设 (f(x,y) = sumlimits_{i=1}^{x} sumlimits_{j=1}^{y}[gcd(i,j)=1]ij)
    发现 (f(x,y)) 的形式与 [HAOI2011] Problem b 中的式子类似,代入 ([gcd(i,j) = 1] = sumlimits_{dmid gcd(i,j)} {mu (d)}) 套路一波:

    [egin{aligned} f(x,y) =& sumlimits_{i=1}^{x} sumlimits_{j=1}^{y}[gcd(i,j)=1]ij\ =& sum_{i=1}^{x} sum_{j=1}^{y}sum_{dmid gcd(i,j)}mu(d) ij\ =& sum_{d=1}mu(d)sum_{i=1}^{x}[dmid i] sum_{j=1}^{y}[dmid j] ij\ =& sum_{d=1}mu(d) d^2sum_{i=1}^{leftlfloor frac{x}{d} ight floor}sum_{j=1}^{leftlfloorfrac{y}{d} ight floor}ij end{aligned}]

    前半部分 (sumlimits_{d=1}mu(d) d^2),可以考虑筛出 (mu(d)) 后求前缀和。
    后半部分是等差数列乘等差数列的形式,设 (g(p,q) = sumlimits_{i=1}^{p} sumlimits_{j=1}^{q}ij)(g_{p,q}) 可以通过下式 (O(1)) 计算:

    [g(p,q) = sum_{i=1}^{p} i sum_{j=1}^{q}j= dfrac{(1 + p) imes p}{2} imes dfrac{(1+q) imes q}{2} ]

    则对于 (f(x,y)),有:

    [f(x,y) = sum_{d=1}mu(d) d^2cdot g(leftlfloor frac{x}{d} ight floor, leftlfloorfrac{y}{d} ight floor) ]

    数论分块求解即可。

    再看回原式,原式等价于:

    [sum_{d=1}dcdot f(leftlfloorfrac{n}{d} ight floor, leftlfloorfrac{m}{d} ight floor) ]

    又是一个可以数论分块求解的形式。
    线性筛预处理后 数论分块套数论分块,复杂度 (O(n + m)),瓶颈是线性筛。


    一些注意的点

    处理时会出现求平方的运算,需要特别注意取模问题,ll 都会爆,被坑惨了。

    在预处理前缀和的这个地方:

    sum[i] = (sum[i - 1] + 1ll * i * i % kMod * (mu[i] + kMod)) % kMod; //仅令 mu + kMod
    

    注意先对平方取模,在乘上 (mu),否则就会爆掉。
    以及可以仅令 (mu + mod)

    以及这个地方:

    int g(int n_, int m_) {
      return (1ll * n_ * (n_ + 1ll) / 2ll % kMod) * (1ll * m_ * (m_ + 1ll) / 2ll % kMod) % kMod;
    }
    

    平方计算,注意随时取模。

    代码

    //知识点:莫比乌斯反演
    /*
    By:Luckyblock
    */
    #include <algorithm>
    #include <cctype>
    #include <cmath>
    #include <cstdio>
    #include <cstring>
    #define ll long long
    const ll kMod = 20101009;
    const int kMaxn = 1e7 + 10;
    //=============================================================
    int pnum, p[kMaxn];
    ll mu[kMaxn], sum[kMaxn];
    bool vis[kMaxn];
    //=============================================================
    inline int read() {
      int f = 1, w = 0;
      char ch = getchar();
      for (; !isdigit(ch); ch = getchar())
        if (ch == '-') f = -1;
      for (; isdigit(ch); ch = getchar()) w = (w << 3) + (w << 1) + (ch ^ '0');
      return f * w;
    }
    void Getmax(int &fir_, int sec_) {
      if (sec_ > fir_) fir_ = sec_;
    }
    void Getmin(int &fir_, int sec_) {
      if (sec_ < fir_) fir_ = sec_;
    }
    void Euler(int lim_) {
      vis[1] = true, mu[1] = 1ll;
      for (int i = 2; i <= lim_; ++ i) {
        if (! vis[i]) {
          p[++ pnum] = i;
          mu[i] = - 1;
        }
        for (int j = 1; j <= pnum && i * p[j] <= lim_; ++ j) {
          vis[i * p[j]] = true;
          if (i % p[j] == 0) { //勿忘平方因子
            mu[i * p[j]] = 0;
            break;
          }
          mu[i * p[j]] = - mu[i];
        }
      }
      sum[1] = 1ll;
      for (int i = 1; i <= lim_; ++ i) {
        sum[i] = (sum[i - 1] + 1ll * i * i % kMod * (mu[i] + kMod)) % kMod; //仅令 mu + kMod
      }
    }
    int g(int n_, int m_) {
      return (1ll * n_ * (n_ + 1ll) / 2ll % kMod) * (1ll * m_ * (m_ + 1ll) / 2ll % kMod) % kMod;
    }
    int f(int n_, int m_) {
      int lim = std :: min(n_, m_), ret = 0;
      for (int l = 1, r; l <= lim; l = r + 1) {
        r = std :: min(n_ / (n_ / l), m_ / (m_ / l));
        ret = (ret + 1ll * (sum[r] - sum[l - 1] + kMod) * g(n_ / l, m_ / l) % kMod) % kMod;
      }
      return ret;
    }
    int Sum(ll l, ll r) {
      return (1ll * (r - l + 1ll) * (l + r) / 2ll) % kMod;
    }
    //=============================================================
    int main() { 
      int n = read(), m = read();
      int lim = std :: min(n, m), ans = 0;
      Euler(lim);
      for (int l = 1, r; l <= lim; l = r + 1) {
        r = std :: min(n / (n / l), m / (m / l));
        ans = (ans + 1ll * Sum(l, r) * f(n / l, m / l) % kMod) % kMod;
      }
      printf("%d", ans);
      return 0;
    }
    /*
    7718820 8445880
    */
    

    SP5971 LCMSUM - LCM Sum

    (T) 次询问,每次询问给定 (n),求:

    [sum_{i=1}^{n}operatorname{lcm}(i,n) ]

    (1<Tle 3 imes 10^5)(1le nle 10^6)

    法一:无脑暴力

    先拆 (operatorname{lcm}),原式等价于:

    [nsum_{i=1}^{n}dfrac{i}{gcd(i,n)} ]

    套路的枚举 (gcd(i,n)),调换枚举顺序,原式等价于:

    [egin{aligned} &nsum_{i=1}sum_{d|i} [gcd(i,n) = d]dfrac{i}{d}\ =& nsum_{i=1}sum_{d|i} [gcd(dfrac{i}{d},dfrac{n}{d}) = 1]dfrac{i}{d}\ =& nsum_{d|n}sum_{i=1}^{n}[d|i][gcd(dfrac{i}{d},dfrac{n}{d}) = 1]dfrac{i}{d} end{aligned}]

    (i,n) 中的 (d) 提出来,变为枚举 (frac{i}{d}),消去整除的条件,原式等价于:

    [nsum_{d|n}sum_{i=1}^{leftlfloorfrac{n}{d} ight floor}[gcd(i,dfrac{n}{d}) = 1]i ]

    代入 ([gcd(i,j) = 1] = sumlimits_{dmid gcd(i,j)} {mu (d)}),原式等价于:

    [nsum_{d|n}sum_{i=1}^{leftlfloorfrac{n}{d} ight floor} i sum_{t|gcd(i,frac{n}{d})}mu (t) ]

    值得注意的是 (t) 的上界为 (frac{n}{d})(dtle n)
    调换枚举顺序,先枚举 (t),原式等价于:

    [nsum_{d|n}sum_{t|frac{n}{d}} mu(t) sum_{i=1}^{leftlfloorfrac{n}{d} ight floor} [t|i] i ]

    套路地消去整除的条件,把 (i) 中的 (t) 提出来,原式等价于:

    [nsum_{d|n}sum_{t|frac{n}{d}} mu(t)t sum_{i=1}^{leftlfloorfrac{n}{dt} ight floor} i ]

    对于最后的一个求和项,设 (g(x) = sumlimits_{i=1}^{x}i = frac{x(x+1)}{2}),显然可以 (O(1)) 求解,原式等价于:

    [nsum_{d|n}sum_{t|frac{n}{d}} mu(t)tcdot g(leftlfloordfrac{n}{dt} ight floor) ]

    考虑枚举 (T = dt),显然 (Tle n)
    (mu(t)t)(d) 无关,可以直接考虑枚举 (t|T),原式等价于:

    [nsum_{T=1}^{n}g(leftlfloordfrac{n}{T} ight floor)sum_{t|T}mu(t)t ]

    前半块是一个数论分块的形式,可以 (O(sqrt{n})) 求解。
    考虑后半块,设 (f(n)=sumlimits_{d|n}mu(d)d),发现它是一个积性函数,可以线性筛筛出,具体地:

    [f(n)= egin{cases} 1-n &nin mathrm{primes} \ f(frac{x}{p}) &p^2mid n\ f(frac{x}{p})f(p) &p^2 mid n end{cases} ]

    其中 (p)(n) 的最小质因子。

    此时已经可以线性筛 + 数论分块求解,复杂度 (O(n+Tsqrt{n})),比较菜鸡,时限 500ms 过不了。
    考虑筛出 (f) 后再用埃氏筛预处理 (sumlimits_{T=1}^{n}g(leftlfloordfrac{n}{T} ight floor)f(T)),输出时乘上 (n),复杂度变为 (O(nlog^2 n + n))


    法二:

    同样先拆 (operatorname{lcm}),枚举 (gcd(i,n)),调换枚举顺序,原式等价于:

    [egin{aligned} &nsum_{i=1}^{n}dfrac{i}{gcd(i,n)}\ =& nsum_{i=1}sum_{d|i} [gcd(i,n) = d]dfrac{i}{d}\ =& nsum_{d|n}sum_{i=1}^{n}[d|i][gcd({i},{n}) = d]dfrac{i}{d} end{aligned}]

    (i,n) 中的 (d) 提出来,变为枚举 (frac{i}{d}),消去整除的条件,原式等价于:

    [egin{aligned} &nsum_{d|n}sum_{i=1}^{leftlfloorfrac{n}{d} ight floor}[gcd(id,n) = d]i\ =& nsum_{d|n}sum_{i=1}^{leftlfloorfrac{n}{d} ight floor}[gcd(i,dfrac{n}{d}) = 1]i end{aligned}]

    调整枚举对象,上式等价于:

    [nsum_{d|n}sum_{i=1}^{d}[gcd(i,d) = 1]i ]

    考虑 (sumlimits_{i=1}^{d}[gcd(i,d) = 1]i) 的实际意义,表示 ([1,d]) 中与 (d) 互质的数的和。
    (d>1) 时,与 (d) 互质的数总是成对存在,即若 (gcd(i,d)=1) 成立,则 (gcd(d-i,d)=1) 成立。
    每对这样的数的和为 (d),共有 (frac{varphi(d)}{2}) 对这样的数。
    则原式等价于:

    [nsum_{d|n}dfrac{varphi(d)d}{2} ]

    可以直接预处理答案。
    预处理时先线性筛出 (varphi),再埃氏筛枚举 (i) 的倍数,令它们的答案加上 (frac{varphi(i)i}{2}),最后输出时乘上 (n)
    复杂度 (O(nlog^2 n + T))


    法二代码

    //知识点:莫比乌斯反演
    /*
    By:Luckyblock
    */
    #include <algorithm>
    #include <cctype>
    #include <cstdio>
    #include <cstring>
    #define ll long long
    const int kMaxn = 1e6;
    //=============================================================
    ll phi[kMaxn + 10], ans[kMaxn + 10];
    int pnum, p[kMaxn + 10];
    bool flag[kMaxn + 10];
    //=============================================================
    inline int read() {
      int f = 1, w = 0;
      char ch = getchar();
      for (; !isdigit(ch); ch = getchar())
        if (ch == '-') f = -1;
      for (; isdigit(ch); ch = getchar()) w = (w << 3) + (w << 1) + (ch ^ '0');
      return f * w;
    }
    void GetMax(int &fir, int sec) {
      if (sec > fir) fir = sec;
    }
    void GetMin(int &fir, int sec) {
      if (sec < fir) fir = sec;
    }
    void GetPrime() {
      phi[1] = 1, flag[1] = true; //注意初始化
      for (int i = 2; i <= kMaxn; ++ i) {
        if (! flag[i]) {
          p[++ pnum] = i;
          phi[i] = i - 1ll;
        }
        for (int j = 1; j <= pnum && i * p[j] <= kMaxn; ++ j) {
          flag[i * p[j]] = true;
          if (i % p[j]) {
            phi[i * p[j]] = phi[i] * phi[p[j]];
          } else {
            phi[i * p[j]] = phi[i] * p[j];
            break;
          }
        }
      }
      for (int i = 1; i <= kMaxn; ++ i) {
        for (int j = 1; i * j <= kMaxn; ++ j) {
          ans[i * j] += (i == 1 ? 1 : 1ll * phi[i] * i / 2);
        }
      }
    }
    //=============================================================
    int main() { 
      GetPrime();
      int T = read();
      while (T --) {
        int n = read();
        printf("%lld
    ", 1ll * ans[n] * n);
      }
      return 0; 
    }
    

    [SDOI2015]约数个数和

    (T) 次询问,每次询问给定 (n,m)
    定义 >(operatorname{d}(i))(i) 的约数个数,求:

    [sum_{i=1}^{n}sum_{j=1}^moperatorname{d}(ij) ]

    (1<T,nle 5 imes 10^4)

    一个结论:

    [operatorname{d}(ij) = sum_{x|i}sum_{y|j}[gcd(x,y)=1] ]

    证明

    先考虑 (i = p^a)(j=p^b(pin mathrm{primes})) 的情况,有:

    [operatorname{d}(p^{a+b})=sum_{x|p^a}sum_{y|p^b}[gcd(x,y)=1] ]

    对于等式左侧,(p^{a+b}) 的约数个数为 (a+b+1)
    对于等式右侧,在保证 (gcd(x,y)=1) 成立的情况下,有贡献的数对 ((x,y)) 只能是下列三种形式:

    • (x>0,y-0)(x)(a) 种取值方法。
    • (x=0,y>0)(y)(b) 种取值方法。
    • (x=0,y=0)

    则等式右侧贡献次数为 (a+b+1) 次,等于 (p^{a+b}) 的约数个数。
    则当 (i = p^a)(j=p^b(pin mathrm{primes})) 时等式成立。

    又不同质因子间相互独立,上述情况可拓展到一般的情况。


    (operatorname{d}(i,j)) 进行化简,代入 ([gcd(i,j) = 1] = sumlimits_{dmid gcd(i,j)} {mu (d)}),原式等价于:

    [egin{aligned} operatorname{d}(ij) =& sum_{x|i}sum_{y|j}[gcd(x,y)=1]\ =& sum_{x|i}sum_{y|j}sumlimits_{dmid gcd(x,y)} {mu (d)} end{aligned}]

    调换枚举顺序,先枚举 (d),原式等价于:

    [sum_{d=1}[d|i][d|j]{mu (d)}sum_{x|i}[d|x]sum_{y|j}[d|y] ]

    把各项中的 (d) 提出来,消去整除的条件,原式等价于:

    [egin{aligned} &sum_{d=1}[d|i][d|j]{mu (d)}sum_{x|frac{i}{d}}sum_{y|frac{j}{d}}1\ =& sum_{d=1}[d|i][d|j]{mu (d)}cdot operatorname{d}(dfrac{i}{d})operatorname{d}(dfrac{j}{d}) end{aligned}]


    (operatorname{d}(ij)) 代回原式,原式等价于:

    [sum_{i=1}^{n}sum_{j=1}^m sum_{d=1}[d|i][d|j]{mu (d)}cdot operatorname{d}(dfrac{i}{d})operatorname{d}(dfrac{j}{d}) ]

    调换枚举顺序,先枚举 (d),原式等价于:

    [sum_{d=1}{mu (d)}sum_{i=1}^{n}[d|i]sum_{j=1}^m [d|j]cdot operatorname{d}(dfrac{i}{d})operatorname{d}(dfrac{j}{d}) ]

    (i,j) 中的 (d) 提出来,变为枚举 (frac{i}{d}, frac{j}{d}),消去整除的条件,原式等价于:

    [sum_{d=1}{mu (d)}sum_{i=1}^{leftlfloorfrac{n}{d} ight floor}operatorname{d}(i)sum_{j=1}^{leftlfloorfrac{m}{d} ight floor}operatorname{d}(j) ]

    考虑预处理 (S(x) = sumlimits_{i=1}^{x}operatorname{d}(i)),则原式等价于:

    [sum_{d=1}{mu (d)}Sleft({leftlfloorfrac{n}{d} ight floor} ight)Sleft({leftlfloorfrac{m}{d} ight floor} ight) ]

    线性筛预处理 (mu,operatorname{d}),数论分块求解即可,复杂度 (O(n+Tsqrt{n}))


    代码

    //知识点:莫比乌斯反演
    /*
    By:Luckyblock
    */
    #include <algorithm>
    #include <cctype>
    #include <cmath>
    #include <cstdio>
    #include <cstring>
    #define ll long long
    const int kMaxn = 5e4 + 10;
    //=============================================================
    int pnum, p[kMaxn];
    ll mu[kMaxn], num[kMaxn], d[kMaxn]; //num 为最小质因子的次数
    ll summu[kMaxn], sumd[kMaxn];
    bool vis[kMaxn];
    //=============================================================
    inline int read() {
      int f = 1, w = 0;
      char ch = getchar();
      for (; !isdigit(ch); ch = getchar())
        if (ch == '-') f = -1;
      for (; isdigit(ch); ch = getchar()) w = (w << 3) + (w << 1) + (ch ^ '0');
      return f * w;
    }
    void Getmax(int &fir_, int sec_) {
      if (sec_ > fir_) fir_ = sec_;
    }
    void Getmin(int &fir_, int sec_) {
      if (sec_ < fir_) fir_ = sec_;
    }
    void Euler(int lim_) {
      vis[1] = true;
      mu[1] = d[1] = 1ll;
      for (int i = 2; i <= lim_; ++ i) {
        if (! vis[i]) {
          p[++ pnum] = i;
          mu[i] = - 1;
          num[i] = 1;
          d[i] = 2;
        }
        for (int j = 1; j <= pnum && i * p[j] <= lim_; ++ j) {
          vis[i * p[j]] = true;
          if (i % p[j] == 0) { //勿忘平方因子
            mu[i * p[j]] = 0;
            num[i * p[j]] = num[i] + 1;
            d[i * p[j]] = 1ll * d[i] / num[i * p[j]] * (num[i * p[j]] + 1ll);
            break;
          }
          mu[i * p[j]] = - mu[i];
          num[i * p[j]] = 1;
          d[i * p[j]] = 2ll * d[i]; //
    
        }
      }
      for (int i = 1; i <= lim_; ++ i) {
        summu[i] = summu[i - 1] + mu[i];
        sumd[i] = sumd[i - 1] + d[i];
      }
    }
    
    //=============================================================
    int main() { 
      Euler(kMaxn - 10);
      int T = read();
      while (T --) {
        int n = read(), m = read(), lim = std :: min(n, m);
        ll ans = 0ll;
        for (int l = 1, r; l <= lim; l = r + 1) {
          r = std :: min(n / (n / l), m / (m / l));
          ans += 1ll * (summu[r] - summu[l - 1]) * sumd[n / l] * sumd[m / l]; //
        }
        printf("%lld
    ", ans);
      }
      return 0;
    }
    /*
    in
    1
    32 43
    */
    /*
    out
    15420
    */
    

    P3768 简单的数学题

    给定参数 (n)(p),求:

    [left(sum_{i=1}^nsum_{j=1}^nicdot jcdot gcd(i,j) ight)mod p ]

    (nleq10^{10})(5 imes10^8leq pleq1.1 imes10^9)(pin mathrm{primes})
    时限 4s。

    无脑套路暴力。

    考虑先枚举 (gcd(i,j)),原式等价于:

    [egin{aligned} &sum_{d=1}dsum_{i=1}^{n}sum_{j=1}^{n}[gcd(i,j)=d]ij\ =& sum_{d=1}dsum_{i=1}^{n}[d|i]sum_{j=1}^{n}[d|j][gcd(dfrac{i}{d},dfrac{j}{d})=1]ij end{aligned}]

    提出 (d),变为枚举 (frac{i}{d})(frac{j}{d}),消去整除的条件,原式等价于:

    [egin{aligned} &sum_{d=1}dsum_{i=1}^{leftlfloorfrac{n}{d} ight floor}sum_{j=1}^{leftlfloorfrac{n}{d} ight floor}[gcd(i,j)=1]idcdot jd\ =& sum_{d=1} d^3sum_{i=1}^{leftlfloorfrac{n}{d} ight floor}sum_{j=1}^{leftlfloorfrac{n}{d} ight floor}[gcd(i,j)=1]ij end{aligned}]

    代入 ([gcd(i,j) = 1] = sumlimits_{dmid gcd(i,j)} {mu (d)}),原式等价于:

    [sum_{d=1} d^3sum_{i=1}^{leftlfloorfrac{n}{d} ight floor}sum_{j=1}^{leftlfloorfrac{n}{d} ight floor}ijsum_{t|gcd(i,j)}mu (t) ]

    值得注意的是 (t) 的上界为 (frac{n}{d})(dtle n)
    调换枚举顺序,先枚举 (t),原式等价于:

    [sum_{d=1} d^3sum_{t=1}mu(t)sum_{i=1}^{leftlfloorfrac{n}{d} ight floor}[t|i] sum_{j=1}^{leftlfloorfrac{n}{d} ight floor}[t|j]ij ]

    和上面一样,提出 (t),套路地消去整除的条件,原式等价于:

    [sum_{d=1} d^3sum_{t=1}mu(t)t^2sum_{i=1}^{leftlfloorfrac{n}{dt} ight floor}sum_{j=1}^{leftlfloorfrac{n}{dt} ight floor}ij ]

    发现后面两个求和是等差数列乘等差数列的形式。
    (g(p,q) = sumlimits_{i=1}^{p} sumlimits_{j=1}^{q}ij)(g_{p,q}) 可以通过下式 (O(1)) 计算:

    [g(p,q) = sum_{i=1}^{p} i sum_{j=1}^{q}j= dfrac{(1 + p) imes p}{2} imes dfrac{(1+q) imes q}{2} ]

    代入原式,原式等价于:

    [sum_{d=1} d^3sum_{t=1}mu(t)t^2cdot gleft({leftlfloorfrac{n}{dt} ight floor},{leftlfloorfrac{n}{dt} ight floor} ight) ]

    考虑枚举 (T = dt),显然 (Tle n)
    再考虑枚举 (d|T),即可得到 (t = frac{T}{d}),原式等价于:

    [egin{aligned} &sum_{T=1}^{n}gleft({leftlfloorfrac{n}{T} ight floor},{leftlfloorfrac{n}{T} ight floor} ight)sum_{d|T}d^3mu{left(dfrac{T}{d} ight)}left(dfrac{T}{d} ight)^2\ =& sum_{T=1}^{n}T^2 gleft({leftlfloorfrac{n}{T} ight floor},{leftlfloorfrac{n}{T} ight floor} ight)sum_{d|T}dcdot mu{left(dfrac{T}{d} ight)} end{aligned}]

    对于后面这一坨,用 (sumlimits_{d|T}dcdot mu{left(frac{T}{d} ight)} = operatorname{Id} ast mu(T)= varphi(T)) 反演,则原式等价于:

    [sum_{T=1}^{n}T^2 varphi(T) cdot gleft({leftlfloorfrac{n}{T} ight floor},{leftlfloorfrac{n}{T} ight floor} ight) ]

    后半块可以数论分块,考虑前半块。
    发现前半段即为 (operatorname{Id}^2(T)varphi(T)),又是前缀和形式,考虑杜教筛。

    有:

    [f(n) = operatorname{Id}^2varphi(n) ]

    考虑找到一个函数 (g),构造函数 (h = fast g) 使其便于求值,有:

    [h(n) = sum_{d|n} d^2varphi(d)cdot gleft(dfrac{n}{d} ight) ]

    看到同时存在 (d)(frac{n}{d}),考虑把 (d^2) 消去。
    (g = operatorname{Id}^2),有:

    [egin{aligned} h(n) =& sum_{d|n} d^2varphi(d)cdot left(dfrac{n}{d} ight)^2\ =& n^2sum_{d|n} varphi(d)\ =& n^2 cdot varphi ast 1(n)\ end{aligned}]

    (varphi ast 1 = operatorname{Id}),则有:

    [h(n) = n^3 ]

    找到了合适的 (g),套杜教筛的公式。

    [egin{aligned} g(1)S(n) &= sum_{i=1}^{n}h(i) - sum_{i=2}^{n}g(i)Sleft(leftlfloordfrac{n}{i} ight floor ight)\ S(n) &= sum_{i=1}^{n} i^3 - sum_{i=2}^{n} i^2cdot Sleft(leftlfloordfrac{n}{i} ight floor ight) end{aligned}]

    前一项是自然数的立方和,有 (sumlimits_{i=1}^{n} i^3 = (frac{n(n+1)}{2})^2)。证明详见:自然数前n项平方和、立方和公式及证明 - 百度文库
    后一项直接等差数列求和 + 数论分块求解即可。


    「SDOI2017」数字表格

    (f_{i}) 表示斐波那契数列的第 (i) 项。
    (T) 组数据,每次给定参数 (n,m),求:

    [prod_{i=1}^{n}prod_{j=1}^{m}f_{gcd(i,j)} pmod {10^9 + 7} ]

    (1le Tle 10^3)(1le n,mle 10^6)
    5S,256MB。

    以下钦定 (nge m)
    大力化式子,先套路地枚举 (gcd(i,j)),用初中知识把两个 (prod) 化到指数位置,原式等于:

    [largeegin{aligned} &prod_{d = 1}^{n}prod_{i=1}^{n}prod_{j=1}^{m}f_{d}[gcd(i,j)=d]\ =&prod_{d = 1}^{n}f_{d}^{left(sumlimits_{i=1}^{n}sumlimits_{j=1}^{m}[gcd(i,j)=d] ight)} end{aligned}]

    分母套路一波,有:

    [egin{aligned} &sumlimits_{i=1}^{n}sumlimits_{j=1}^{m}[gcd(i,j) = d]\ =& sumlimits_{i=1}^{leftlfloorfrac{n}{d} ight floor}sumlimits_{j=1}^{leftlfloorfrac{m}{d} ight floor}[gcd(i,j) = 1]\ =& sumlimits_{i=1}^{leftlfloorfrac{n}{d} ight floor}sumlimits_{j=1}^{leftlfloorfrac{m}{d} ight floor}sum_{kmid gcd(i,j)}mu (k)\ =& sum_{k=1}mu(k)leftlfloordfrac{n}{kd} ight floorleftlfloordfrac{m}{kd} ight floor end{aligned}]

    代回原式,原式等于:

    [large egin{aligned} &prod_{d = 1}^{n}f_{d}^{left(sumlimits_{k=1}mu(k)leftlfloorfrac{n}{kd} ight floorleftlfloorfrac{m}{kd} ight floor ight)}\ =& prod_{d = 1}^{n}left(f_{d}^{{sumlimits_{k=1}mu(k)}} ight)^{leftlfloorfrac{n}{kd} ight floorleftlfloorfrac{m}{kd} ight floor} end{aligned}]

    考虑再暴力拆一波,原式等于:

    [large egin{aligned} &prod_{d = 1}^{n}left(f_{d}^{{sumlimits_{k=1}mu(k)}} ight)^{leftlfloorfrac{n}{kd} ight floorleftlfloorfrac{m}{kd} ight floor}\ =& prod_{d = 1}^{n}left(prod_{k=1}^{leftlfloorfrac{n}{d} ight floor}f_{d}^{mu(k)leftlfloorfrac{n}{kd} ight floorleftlfloorfrac{m}{kd} ight floor} ight) end{aligned}]

    做不动了,但发现变量仅有 (k,d,kd),考虑更换枚举对象改为枚举 (t = kd)(d),则原式等于:

    [largeprod_{t=1}^{n}left(prod_{d | t}^{n}f_{d}^{{mu(frac{t}{d})}} ight)^{leftlfloorfrac{n}{t} ight floorleftlfloorfrac{m}{t} ight floor} ]

    枚举对象变成了约数形式。从后面的式子推前面的式子是比较显然的,可以认为这种枚举 (t=kd) 的形式是一种逆向操作。

    设:

    [large g(t)=prod_{d | t}^{n}f_{d}^{{mu(frac{t}{d})}} ]

    (g(t)) 可以用类似埃氏筛的方法 (O(nlog ^2 n)) 地预处理出来。再把 (g) 代回原式,原式等于:

    [largeprod_{t=1}^{n}g(t)^{leftlfloorfrac{n}{t} ight floorleftlfloorfrac{m}{t} ight floor} ]

    可以考虑预处理 (g(t)) 的前缀积,数论分块枚举指数求解即可。

    总时间复杂度 (O(nlog ^2 n + Tsqrt n)),轻微卡常可以过。

    //知识点:莫比乌斯反演 
    /*
    By:Luckyblock
    */
    #include <algorithm>
    #include <cctype>
    #include <cstdio>
    #include <cstring>
    #define LL long long
    const LL mod = 1e9 + 7;
    const int kN = 1e6;
    //=============================================================
    LL n, m, ans;
    int p_num, p[kN + 10];
    bool vis[kN + 10];
    LL mu[kN + 10], f[kN + 10], g[kN + 10], prod[kN + 10];
    LL invf[kN + 10], invp[kN];
    //=============================================================
    inline int read() {
      int f = 1, w = 0;
      char ch = getchar();
      for (; !isdigit(ch); ch = getchar())
        if (ch == '-') f = -1;
      for (; isdigit(ch); ch = getchar()) {
        w = (w << 3) + (w << 1) + (ch ^ '0');
      }
      return f * w;
    }
    void Chkmax(int &fir_, int sec_) {
      if (sec_ > fir_) fir_ = sec_;
    }
    void Chkmin(int &fir_, int sec_) {
      if (sec_ < fir_) fir_ = sec_;
    }
    LL QPow(LL x_, LL y_) {
      x_ %= mod;
      LL ret = 1;
      for (; y_; y_ >>= 1ll) {
        if (y_ & 1) ret = ret * x_ % mod;
        x_ = x_ * x_ % mod;
      }
      return ret;
    }
    void Euler() {
      vis[1] = true, mu[1] = 1;
      for (int i = 2; i <= kN; ++ i) {
        if (! vis[i]) {
          p[++ p_num] = i;
          mu[i] = -1;
        }
        for (int j = 1; j <= p_num && i * p[j] <= kN; ++ j) {
          vis[i * p[j]] = true;
          if (i % p[j] == 0) {
            mu[i * p[j]] = 0;
            break;
          }
          mu[i * p[j]] = -mu[i];
        }
      }
    }
    void Prepare() {
      g[1] = g[2] = 1;
      f[1] = f[2] = 1;
      invf[1] = invf[2] = 1;
      for (int i = 3; i <= kN; ++ i) {
        g[i] = 1;
        f[i] = (f[i - 1] + f[i - 2]) % mod;
        invf[i] = QPow(f[i], mod - 2);
      }
      
      Euler();
      for (int d = 1; d <= kN; ++ d) {
        for (int j = 1; d * j <= kN; ++ j) {
          if (mu[j] == 1) {
            g[d * j] = g[d * j] * f[d] % mod;
          } else if (mu[j] == -1) {
            g[d * j] = g[d * j] * invf[d] % mod;
          }
        }
      }
      invp[0] = prod[0] = 1;
      for (int i = 1; i <= kN; ++ i) {
        prod[i] = prod[i - 1] * g[i] % mod;
        invp[i] = QPow(prod[i], mod - 2);
      }
    }
    //=============================================================
    int main() {
      Prepare();
      int T = read();
      while (T -- ){
        n = read(), m = read(), ans = 1;
        if (n < m) std::swap(n, m);
        for (LL l = 1, r = 1; l <= m; l = r + 1) {
          r = std::min(n / (n / l), m / (m / l));
          ans = (ans * QPow(prod[r] * invp[l - 1] % mod, 1ll * (n / l) * (m / l))) % mod;
        }
        printf("%lld
    ", ans);
      }
      return 0;
    }
    

    P5518 [MtOI2019]幽灵乐团 / 莫比乌斯反演基础练习题

    给定参数 (p),有 (T) 组数据,每次给定参数 (A,B,C),求:

    [prod_{i=1}^{A}prod_{j=1}^{B}prod_{k=1}^{C}left(dfrac{operatorname{lcm}(i,j)}{gcd(i,k)} ight)^{f(type)} ]

    其中 (f(type)) 的取值如下:

    [f(type) = egin{cases} 1 &type = 0\ i imes j imes k &type = 1\ gcd(i,j,k) &type = 2 end{cases}]

    (1le A,B,Cle 10^5)(10^7le ple 1.05 imes 10^9)(pin mathbb{P})(T=70)
    2.5S,128MB。

    先化下原式,原式等于:

    [prod_{i=1}^{A}prod_{j=1}^{B}prod_{k=1}^{C}left(dfrac{i imes j }{gcd(i,j) imes gcd(i,k)} ight)^{f(type)} ]

    发现每一项仅与两个变量有关,设:

    [egin{aligned} f_1(a,b,c) &= prod_{i=1}^{a}prod_{j=1}^{b}prod_{k=1}^{c} i^{f(type)}\ f_2(a,b,c) &= prod_{i=1}^{a}prod_{j=1}^{b}prod_{k=1}^{c} gcd(i,j)^{f(type)} end{aligned}]

    发现 (prod) 可以随意交换,则原式等价于:

    [dfrac{f_1(a,b,c) imes f_1(b,a,c)}{f_2(a,b,c) imes f_2(a,c,b)} ]

    考虑在 (type) 取值不同时,如何快速求得 (f_1)(f_2)
    一共有 (6) 个需要推导的式子,不妨就叫它们 (1sim 6) 面了(


    type = 0

    [egin{aligned} f_1(a,b,c) &= prod_{i=1}^{a}prod_{j=1}^{b}prod_{k=1}^{c} i\ f_2(a,b,c) &= prod_{i=1}^{a}prod_{j=1}^{b}prod_{k=1}^{c} gcd(i,j) end{aligned}]

    对于 1 面,显然有:

    [prod_{i=1}^{a}prod_{j=1}^{b}prod_{k=1}^{c} i = left(prod_{i=1}^{a}i ight)^{b imes c} ]

    预处理阶乘 + 快速幂即可,单次计算时间复杂度 (O(log n))


    再考虑 2 面,套路地枚举 (gcd),显然有:

    [large egin{aligned} &prod_{i=1}^{a}prod_{j=1}^{b}prod_{k=1}^{c} gcd(i,j)\ =&left(prod_{i=1}^{a}prod_{j=1}^{b}gcd(i,j) ight)^c\ =& left(prod_{d=1} d^{left(sumlimits_{i=1}^{a}sumlimits_{j=1}^{b}[gcd(i,j) = d] ight)} ight)^c end{aligned}]

    指数是个套路,可以看这里:P3455 [POI2007]ZAP-Queries。于是有:

    [egin{aligned} &sumlimits_{i=1}^{a}sumlimits_{j=1}^{b}[gcd(i,j) = d]\ =& sumlimits_{i=1}^{leftlfloorfrac{a}{d} ight floor}sumlimits_{j=1}^{leftlfloorfrac{b}{d} ight floor}[gcd(i,j) = 1]\ =& sumlimits_{i=1}^{leftlfloorfrac{a}{d} ight floor}sumlimits_{j=1}^{leftlfloorfrac{b}{d} ight floor}sum_{kmid gcd(i,j)}mu (k)\ =& sum_{k=1}mu(k)leftlfloordfrac{a}{kd} ight floorleftlfloordfrac{b}{kd} ight floor end{aligned}]

    代回原式,略做处理,则原式等于:

    [large egin{aligned} &left(prod_{d=1} d^{left(sumlimits_{k=1}mu(k)leftlfloorfrac{a}{kd} ight floorleftlfloorfrac{b}{kd} ight floor ight)} ight)^c\ =& left(prod_{d=1} left(d^{sumlimits_{k=1}mu(k)} ight)^{leftlfloorfrac{a}{kd} ight floorleftlfloorfrac{b}{kd} ight floor} ight)^c\ =& prod_{d=1} left(prod_{k=1}^{leftlfloorfrac{n}{d} ight floor}d^{left(mu(k)leftlfloorfrac{a}{kd} ight floorleftlfloorfrac{b}{kd} ight floor ight)} ight)^c end{aligned}]

    「SDOI2017」数字表格 一样,考虑枚举 (t=kd)(d),则原式等于:

    [large prod_{t=1}^{n}left(left(prod_{d|t} d^{mu{left(frac{t}{d} ight)}} ight)^{leftlfloorfrac{a}{t} ight floorleftlfloorfrac{b}{t} ight floor} ight)^c ]

    设:

    [large g_0(t) = prod_{d|t}d^{muleft(frac{t}{d} ight)} ]

    线性筛预处理 (mu) 后,(g_0(t)) 可以用埃氏筛预处理,时间复杂度 (O(nlog n))。再代回原式,原式等于:

    [large prod_{t=1}^{a}left(g_0(t)^{leftlfloorfrac{a}{t} ight floorleftlfloorfrac{b}{t} ight floor} ight)^c ]

    预处理 (g_0(t)) 的前缀积和前缀积的逆元,复杂度 (O(nlog n))
    数论分块 + 快速幂计算即可,单次时间复杂度 (O(sqrt nlog n))


    type = 1

    [egin{aligned} f_1(a,b,c) &= prod_{i=1}^{a}prod_{j=1}^{b}prod_{k=1}^{c} i^{i imes j imes k}\ f_2(a,b,c) &= prod_{i=1}^{a}prod_{j=1}^{b}prod_{k=1}^{c} gcd(i,j)^{i imes j imes k} end{aligned}]

    考虑 3 面,把 (prod k) 扔到指数位置,有:

    [large prod_{i=1}^{a}prod_{j=1}^{b}prod_{k=1}^{c} i^{i imes j imes k} = prod_{i=1}^{a}prod_{j=1}^{b}i^{left(i imes j imes sumlimits_{k = 1}^{c} k ight)} ]

    再把 (prod j) 也扔到指数位置,引入 (operatorname{sum}(n) = sum_{i=1}^{n} i = frac{n(n+1)}{2}),原式等于:

    [left(prod_{i=1}^{a}i^i ight)^{operatorname{sum}(b) imes operatorname{sum}(c)} ]

    预处理 (i^i) 的前缀积,复杂度 (O(nlog n))
    指数可以 (O(1)) 算出,然后快速幂,单次时间复杂度 (O(log n))

    根据费马小定理,指数需要对 (p - 1) 取模。注意 (p-1) 不是质数,计算 (operatorname{sum}) 时不能用逆元,但乘不爆 LL,直接算就行。


    再考虑 4 面,发现 (k)(gcd) 无关,则同样把 (prod k) 扔到指数位置,则有:

    [large egin{aligned} &prod_{i=1}^{a}prod_{j=1}^{b}prod_{k=1}^{c} gcd(i,j)^{i imes j imes k}\ =& left(prod_{i=1}^aprod_{j=1}^bgcd(i,j)^{i imes j} ight)^{operatorname{sum}(c)} end{aligned}]

    套路地枚举 (gcd),原式等于:

    [large left(prod_{d=1}d^{left(sumlimits_{i=1}^a sumlimits_{j=1}^b i imes j[gcd(i,j)=d] ight)} ight)^{operatorname{sum}(c)} ]

    大力化简指数,有:

    [large egin{aligned} &sumlimits_{i=1}^a sumlimits_{j=1}^b i imes j[gcd(i,j)=d]\ =& d^2 sumlimits_{i=1}^{leftlfloorfrac{a}{d} ight floor} sumlimits_{j=1}^{leftlfloorfrac{b}{d} ight floor} i imes j[gcd(i,j)=1\ =& d^2 sumlimits_{i=1}^{leftlfloorfrac{a}{d} ight floor} i sumlimits_{j=1}^{leftlfloorfrac{b}{d} ight floor} jsumlimits_{t|gcd(i,j)}mu(t)\ =& d^2 sumlimits_{i=1}^{leftlfloorfrac{a}{d} ight floor} i sumlimits_{j=1}^{leftlfloorfrac{b}{d} ight floor} jsumlimits_{k|gcd(i,j)}mu(k)\ =& d^2 sumlimits_{k=1}mu(k)sumlimits_{i=1}^{leftlfloorfrac{a}{d} ight floor} i[k|i] sumlimits_{j=1}^{leftlfloorfrac{b}{d} ight floor} j[k|j]\ =& d^2 sumlimits_{k=1}k^2mu(k)sumlimits_{i=1}^{leftlfloorfrac{a}{kd} ight floor} isumlimits_{j=1}^{leftlfloorfrac{b}{kd} ight floor} j\ =& d^2sumlimits_{k=1}k^2mu(k)operatorname{sum}{left(leftlfloorfrac{a}{kd} ight floor ight)} operatorname{sum}{left(leftlfloorfrac{b}{kd} ight floor ight)}\ end{aligned}]

    指数化不动了,代回原式,原式等于:

    [large left(prod_{d=1}d^{left(d^2sumlimits_{k=1}k^2mu(k)operatorname{sum}{left(leftlfloorfrac{a}{kd} ight floor ight)} operatorname{sum}{left(leftlfloorfrac{b}{kd} ight floor ight)} ight)} ight)^{operatorname{sum}(c)} ]

    同 2 面的情况,先展开一下,再枚举 (t=kd)(d),原式等于:

    [large egin{aligned} &left(prod_{d=1}left(prod_{k=1}^{leftlfloorfrac{n}{d} ight floor}d^{left(d^2 k^2mu(k) ight)} ight)^{left(operatorname{sum}{left(leftlfloorfrac{a}{kd} ight floor ight)} operatorname{sum}{left(leftlfloorfrac{b}{kd} ight floor ight)} ight)} ight)^{operatorname{sum}(c)}\ =& prod_{t=1}left(left(prod_{d|t}d^{left(d^2left(frac{t}{d} ight)^2muleft(frac{t}{d} ight) ight)} ight)^{operatorname{sum}{left(leftlfloorfrac{a}{t} ight floor ight)} operatorname{sum}{left(leftlfloorfrac{b}{t} ight floor ight)}} ight)^{operatorname{sum}(c)}\ =& prod_{t=1}left(left(prod_{d|t}d^{left(t^2muleft(frac{t}{d} ight) ight)} ight)^{operatorname{sum}{left(leftlfloorfrac{a}{t} ight floor ight)} operatorname{sum}{left(leftlfloorfrac{b}{t} ight floor ight)}} ight)^{operatorname{sum}(c)} end{aligned}]

    与二面相同,设:

    [large g_1(t) = prod_{d|t}d^{left(t^2muleft(frac{t}{d} ight) ight)} ]

    (g_1(t)) 可以用埃氏筛套快速幂筛出,时间复杂度 (O(nlog^2 n))。再代回原式,原式等于:

    [prod_{t=1}left(g_1(t)^{operatorname{sum}{left(leftlfloorfrac{a}{t} ight floor ight)} operatorname{sum}{left(leftlfloorfrac{b}{t} ight floor ight)}} ight)^{operatorname{sum}(c)} ]

    同样预处理 (g_1(t)) 的前缀积及其逆元,时间复杂度 (O(nlog n))
    整除分块 + 快速幂即可,单次时间复杂度 (O(sqrt nlog n))

    注意指数的取模。


    type = 2

    [egin{aligned} f_1(a,b,c) &= prod_{i=1}^{a}prod_{j=1}^{b}prod_{k=1}^{c} i^{gcd(i,j,k)}\ f_2(a,b,c) &= prod_{i=1}^{a}prod_{j=1}^{b}prod_{k=1}^{c} gcd(i,j)^{gcd(i,j,k)} end{aligned}]

    考虑 5 面,手段同上,大力反演化简一波,再调换枚举对象,则有:

    [large egin{aligned} &prod_{i=1}^{a}prod_{j=1}^{b}prod_{k=1}^{c} i^{gcd(i,j,k)}\ =&prod_{d=1}prodlimits_{i=1}^{a}i^{left(sumlimits_{j=1}^{b}sumlimits_{k=1}^{c}[gcd(i,j,k)=d] ight)}\ =& prod_{d=1}prodlimits_{i=1}^{a}i^{left(sumlimits_{j=1}^{leftlfloorfrac{b}{d} ight floor}sumlimits_{k=1}^{leftlfloorfrac{c}{d} ight floor}[gcd(frac{i}{d},j,k)=1] ight)}\ =& prod_{d=1}prodlimits_{i=1}^{leftlfloorfrac{a}{d} ight floor}(id)^{left(dsumlimits_{j=1}^{leftlfloorfrac{b}{d} ight floor}sumlimits_{k=1}^{leftlfloorfrac{c}{d} ight floor}[gcd(i,j,k)=1] ight)}\ =& prod_{d=1}prodlimits_{i=1}^{leftlfloorfrac{a}{d} ight floor}(id)^{left(dsumlimits_{j=1}^{leftlfloorfrac{b}{d} ight floor}sumlimits_{k=1}^{leftlfloorfrac{c}{d} ight floor}sumlimits_{x|gcd(i,j,k)}{mu(x)} ight)}\ =& prod_{d=1}prodlimits_{i=1}^{leftlfloorfrac{a}{d} ight floor}(id)^{left(dsumlimits_{x=1}mu(x)[x|i]sumlimits_{j=1}^{leftlfloorfrac{b}{d} ight floor}[x|j]sumlimits_{k=1}^{leftlfloorfrac{c}{d} ight floor}[x|k] ight)}\ =& prod_{d=1}prod_{x=1}prodlimits_{i=1}^{leftlfloorfrac{a}{d} ight floor}(id)^{left(d imes mu(x)[x|i]sumlimits_{j=1}^{leftlfloorfrac{b}{d} ight floor}[x|j]sumlimits_{k=1}^{leftlfloorfrac{c}{d} ight floor}[x|k] ight)}\ =& prod_{d=1}prod_{x=1}prodlimits_{i=1}^{leftlfloorfrac{a}{xd} ight floor}(ixd)^{left(d imes mu(x){leftlfloorfrac{b}{xd} ight floor}{leftlfloorfrac{c}{xd} ight floor} ight)}\ =& prod_{t = 1}prod_{d|T}prod_{i=1}^{leftlfloorfrac{a}{t} ight floor}(it)^{left(d imes muleft(frac{t}{d} ight){leftlfloorfrac{b}{t} ight floor}{leftlfloorfrac{c}{t} ight floor} ight)}\ =& prod_{t = 1}prod_{d|T}left(t^{leftlfloorfrac{a}{t} ight floor}prod_{i=1}^{leftlfloorfrac{a}{t} ight floor}i ight)^{d imes muleft(frac{t}{d} ight){leftlfloorfrac{b}{t} ight floor}{leftlfloorfrac{c}{t} ight floor}}\ end{aligned}]

    引入 (operatorname{fac}(n) = prod_{i=1}^{n} i),再根据枚举对象调整一下指数,原式等于:

    [large egin{aligned} &prod_{t = 1}prod_{d|t}left(t^{leftlfloorfrac{a}{t} ight floor} imes operatorname{fac}left(leftlfloorfrac{a}{t} ight floor ight) ight)^{left(d imes muleft(frac{t}{d} ight){leftlfloorfrac{b}{t} ight floor}{leftlfloorfrac{c}{t} ight floor} ight)}\ =& prod_{t = 1}left(prod_{d|t}left(t^{leftlfloorfrac{a}{t} ight floor} imes operatorname{fac}left(leftlfloorfrac{a}{t} ight floor ight) ight)^{d imes muleft(frac{t}{d} ight)} ight)^{{leftlfloorfrac{b}{t} ight floor}{leftlfloorfrac{c}{t} ight floor}}\ =& prod_{t = 1}left(left(t^{leftlfloorfrac{a}{t} ight floor} imes operatorname{fac}left(leftlfloorfrac{a}{t} ight floor ight) ight)^{sumlimits_{d|t}d imes muleft(frac{t}{d} ight)} ight)^{{leftlfloorfrac{b}{t} ight floor}{leftlfloorfrac{c}{t} ight floor}} end{aligned}]

    指数中出现了一个经典的狄利克雷卷积的形式,对其进行反演。
    ((operatorname{Id}ast mu) (n)= varphi (n)) 代入原式,原式等于:

    [large egin{aligned} &prod_{t = 1}left(t^{leftlfloorfrac{a}{t} ight floor} imes operatorname{fac}left(leftlfloorfrac{a}{t} ight floor ight) ight)^{varphi(t){leftlfloorfrac{b}{t} ight floor}{leftlfloorfrac{c}{t} ight floor}}\ =& prod_{t = 1}left(t^{varphi(t)leftlfloorfrac{a}{t} ight floor{leftlfloorfrac{b}{t} ight floor}{leftlfloorfrac{c}{t} ight floor}} imes operatorname{fac}left(leftlfloorfrac{a}{t} ight floor ight)^{varphi(t){leftlfloorfrac{b}{t} ight floor}{leftlfloorfrac{c}{t} ight floor}} ight) end{aligned}]

    预处理 (t^{varphi(t)}) 的前缀积及逆元,阶乘的前缀积及阶乘逆元,(pmod {p-1}) 下的 (varphi(t)) 的前缀和(指数
    ),时间复杂度 (O(nlog n))
    同样整除分块 + 快速幂即可,单次时间复杂度 (O(sqrt nlog n))


    然后是最掉 sans 的 6 面。有 (gcd(i,j,k) = gcd(gcd(i,j), k)),考虑先枚举 (gcd(i,j)),然后套路化式子,则有:

    [large egin{aligned} &prod_{i=1}^{a}prod_{j=1}^{b}prod_{k=1}^{c} gcd(i,j)^{gcd(i,j,k)}\ =& prod_{d=1}prod_{i=1}^{a}prod_{j=1}^{b}prod_{k=1}^{c} [gcd(i,j)=d] d^{gcd(d,k)}\ =& prod_{d=1} left(d^{left(sumlimits_{i=1}^{a}sumlimits_{j=1}^{b} [gcd(i,j)=d] ight)} ight)^{sumlimits_{k=1}^{c}gcd(d,k)} end{aligned}]

    先考虑最外面的指数,这也是个套路,可以参考 一个例子。用 (operatorname{Id} = varphi ast 1) 反演,显然有:

    [large egin{aligned} &sumlimits_{k=1}^{c}gcd(d,k)\ =& sumlimits_{k=1}^{c}sum_{x|gcd(d,k)}varphi(x)\ =& sum_{x=1}varphi(x)[x|d]sum_{k=1}^{c}[x|k]\ =& sum_{x|d}varphi(x)leftlfloorfrac{c}{x} ight floor end{aligned}]

    再考虑里面的指数,发现这式子在 2 面已经推了一遍了,于是直接拿过来用,有:

    [sumlimits_{i=1}^{a}sumlimits_{j=1}^{b}[gcd(i,j) = d]=sum_{y=1}mu(y)leftlfloordfrac{a}{yd} ight floorleftlfloordfrac{b}{yd} ight floor ]

    将化简后的两个指数代入原式,原式等于:

    [large egin{aligned} &prod_{d=1} left(d^{left(sumlimits_{y=1}mu(y)leftlfloorfrac{a}{yd} ight floorleftlfloorfrac{b}{yd} ight floor ight)} ight)^{sumlimits_{x|d}varphi(x)leftlfloorfrac{c}{x} ight floor}\ =& prod_{d=1} left(prodlimits_{y=1}d^{left(mu(y)leftlfloorfrac{a}{yd} ight floorleftlfloorfrac{b}{yd} ight floor ight)} ight)^{sumlimits_{x|d}varphi(x)leftlfloorfrac{c}{x} ight floor} end{aligned}]

    与 2、4 面同样套路地,考虑枚举 (t=yd)(d),再略作调整,原式等于:

    [large egin{aligned} &prod_{d=1} left(prodlimits_{y=1}d^{left(mu(y)leftlfloorfrac{a}{yd} ight floorleftlfloorfrac{b}{yd} ight floor ight)} ight)^{sumlimits_{x|d}varphi(x)leftlfloorfrac{c}{x} ight floor}\ =& prod_{t=1}prod_{d|t} d^{left(muleft(frac{t}{d} ight)leftlfloorfrac{a}{t} ight floorleftlfloorfrac{b}{t} ight floorsumlimits_{x|d}varphi(x)leftlfloorfrac{c}{x} ight floor ight)}\ =& prod_{t=1}left(prod_{d|t} d^{left(muleft(frac{t}{d} ight)sumlimits_{x|d}varphi(x)leftlfloorfrac{c}{x} ight floor ight)} ight)^{leftlfloorfrac{a}{t} ight floorleftlfloorfrac{b}{t} ight floor}\ =& prod_{t=1}left(prod_{d|t} prod_{x|d}d^{left(muleft(frac{t}{d} ight)varphi(x)leftlfloorfrac{c}{x} ight floor ight)} ight)^{leftlfloorfrac{a}{t} ight floorleftlfloorfrac{b}{t} ight floor} end{aligned}]

    发现要同时枚举 (d)(x),化不动了。
    从题解里学到一个比较神的技巧,考虑把 (d) 拆成 (x)(frac{d}{x}) 分别计算贡献再相乘,即分别计算下两式:

    [large egin{aligned} &prod_{t=1}left(prod_{d|t} prod_{x|d}x^{left(muleft(frac{t}{d} ight)varphi(x)leftlfloorfrac{c}{x} ight floor ight)} ight)^{leftlfloorfrac{a}{t} ight floorleftlfloorfrac{b}{t} ight floor}\ &prod_{t=1}left(prod_{d|t} prod_{x|d}{left(frac{d}{x} ight)}^{left(muleft(frac{t}{d} ight)varphi(x)leftlfloorfrac{c}{x} ight floor ight)} ight)^{leftlfloorfrac{a}{t} ight floorleftlfloorfrac{b}{t} ight floor} end{aligned}]


    先考虑 (x) 的情况,首先把枚举 (x) 调整到最外层。设 (operatorname{lim}=max(a,b,c)),则原式等于:

    [large egin{aligned} &prod_{x=1} prod_{t=1}^{operatorname{lim}}[x|t]left(prod_{d|t} [x|d]{x}^{left(muleft(frac{t}{d} ight)varphi(x)leftlfloorfrac{c}{x} ight floor ight)} ight)^{leftlfloorfrac{a}{t} ight floorleftlfloorfrac{b}{t} ight floor}\ =& prod_{x=1} prod_{t=1}^{leftlfloorfrac{operatorname{lim}}{x} ight floor}left(prod_{d|t} {x}^{left(muleft(frac{tx}{dx} ight)varphi(x)leftlfloorfrac{c}{x} ight floor ight)} ight)^{leftlfloorfrac{a}{tx} ight floorleftlfloorfrac{b}{tx} ight floor}\ =& prod_{x=1} prod_{t=1}^{leftlfloorfrac{operatorname{lim}}{x} ight floor}prod_{d|t} {x}^{left(varphi(x)leftlfloorfrac{c}{x} ight floorleftlfloorfrac{a}{tx} ight floorleftlfloorfrac{b}{tx} ight floormuleft(frac{t}{d} ight) ight)} end{aligned}]

    (prod {t}) 挪到指数位置,原式等于:

    [large egin{aligned} &prod_{x=1} {x}^{left(sumlimits_{t=1}^{leftlfloorfrac{operatorname{lim}}{x} ight floor}sumlimits_{d|t}varphi(x)leftlfloorfrac{c}{x} ight floorleftlfloorfrac{a}{tx} ight floorleftlfloorfrac{b}{tx} ight floormuleft(frac{t}{d} ight) ight)}\ =& prod_{x=1} {x}^{left(varphi(x)leftlfloorfrac{c}{x} ight floorsumlimits_{t=1}^{leftlfloorfrac{operatorname{lim}}{x} ight floor}leftlfloorfrac{a}{tx} ight floorleftlfloorfrac{b}{tx} ight floorsumlimits_{d|t}muleft(frac{t}{d} ight) ight)} end{aligned}]

    指数中又出现了一个经典的狄利克雷卷积的形式,对其进行反演。
    ((mu ast 1) (n)= epsilon (n)=[n=1]) 代入原式,原式等于:

    [large egin{aligned} &prod_{x=1} {x}^{left(varphi(x)leftlfloorfrac{c}{x} ight floorsumlimits_{t=1}^{leftlfloorfrac{operatorname{lim}}{x} ight floor}leftlfloorfrac{a}{tx} ight floorleftlfloorfrac{b}{tx} ight floor[t=1] ight)}\ =& prod_{x=1} {x}^{left(varphi(x)leftlfloorfrac{a}{x} ight floorleftlfloorfrac{b}{x} ight floorleftlfloorfrac{c}{x} ight floor ight)} end{aligned}]

    得到了一个非常优美的式子,而且发现这个式子是 5 面最终答案的一部分。同 5 面的做法,直接整除分块即可。


    再考虑 (frac{d}{x}) 的情况,同上先把枚举 (x) 放到最外层,并调整一下指数,则原式等于:

    [large egin{aligned} &prod_{x=1} prod_{t=1}^{operatorname{lim}}[x|t]left(prod_{d|t} [x|d]{left(frac{d}{x} ight)}^{left(muleft(frac{t}{d} ight)varphi(x)leftlfloorfrac{c}{x} ight floor ight)} ight)^{leftlfloorfrac{a}{t} ight floorleftlfloorfrac{b}{t} ight floor}\ =& prod_{x=1} prod_{t=1}^{leftlfloorfrac{operatorname{lim}}{x} ight floor}left(prod_{d|tx} [x|d]{left(frac{d}{x} ight)}^{left(muleft(frac{tx}{d} ight)varphi(x)leftlfloorfrac{c}{x} ight floor ight)} ight)^{leftlfloorfrac{a}{tx} ight floorleftlfloorfrac{b}{tx} ight floor}\ =& prod_{x=1} left(prod_{t=1}^{leftlfloorfrac{operatorname{lim}}{x} ight floor}left(prod_{d|tx} [x|d]{left(frac{d}{x} ight)}^{muleft(frac{tx}{d} ight)} ight)^{leftlfloorfrac{a}{tx} ight floorleftlfloorfrac{b}{tx} ight floor} ight)^{varphi(x)leftlfloorfrac{c}{x} ight floor} end{aligned}]

    考虑枚举 (dx),替换原来的 (d),注意一下这里的倍数关系。原式等于:

    [large prod_{x=1} left(prod_{t=1}^{leftlfloorfrac{operatorname{lim}}{x} ight floor}left(prod_{d|t}d^{muleft(frac{t}{d} ight)} ight)^{leftlfloorfrac{a}{tx} ight floorleftlfloorfrac{b}{tx} ight floor} ight)^{varphi(x)leftlfloorfrac{c}{x} ight floor} ]

    发现最内层的式子 (prod_{d|t}d^{muleft(frac{t}{d} ight)}),即为二面处理过的 (g_0(t))。直接代入,原式等于:

    [large prod_{x=1} left(prod_{t=1}^{leftlfloorfrac{operatorname{lim}}{x} ight floor}g_0(t)^{leftlfloorfrac{a}{tx} ight floorleftlfloorfrac{b}{tx} ight floor} ight)^{varphi(x)leftlfloorfrac{c}{x} ight floor} ]

    一个小结论,证明可以看 这里

    [forall a,b,cin mathbb{Z},leftlfloordfrac{a}{bc} ight floor = leftlfloor{dfrac{leftlfloordfrac{a}{b} ight floor}{c}} ight floor ]

    则原式等于:

    [large prod_{x=1} left(prod_{t=1}^{leftlfloorfrac{operatorname{lim}}{x} ight floor}g_0(t)^{leftlfloorfrac{leftlfloorfrac{a}{x} ight floor}{t} ight floorleftlfloorfrac{leftlfloorfrac{b}{x} ight floor}{t} ight floor} ight)^{varphi(x)leftlfloorfrac{c}{x} ight floor} ]

    于是可以先对外层整除分块,再对内层整除分块。

    然后就做完了,哈哈哈。


    一些实现上的小技巧:

    • 逆元能预处理就预处理。
    • 注意对指数取模,模数为 (p-1)
    //知识点:莫比乌斯反演 
    /*
    By:Luckyblock
    用了比较清晰易懂的变量名,阅读体验应该会比较好。  
    vsc 的自动补全真是太好啦!
    */
    #include <algorithm>
    #include <cctype>
    #include <cstdio>
    #include <cstring>
    using std::min;
    using std::max;
    #define LL long long
    const int Lim = 1e5;
    const int kN = 1e5 + 10;
    //=============================================================
    LL A, B, C, mod, ans;
    int T, p_num, p[kN];
    bool vis[kN];
    LL mu[kN], phi[kN], fac[kN], g[2][kN];
    LL sumphi[kN], prodt_phi[kN], prodi_i[kN], prodg[2][kN];
    LL inv[kN], inv_fac[kN], inv_prodt_phi[kN], inv_prodg[2][kN];
    //=============================================================
    inline int read() {
      int f = 1, w = 0;
      char ch = getchar();
      for (; !isdigit(ch); ch = getchar())
        if (ch == '-') f = -1;
      for (; isdigit(ch); ch = getchar()) {
        w = (w << 3) + (w << 1) + (ch ^ '0');
      }
      return f * w;
    }
    void Chkmax(int &fir_, int sec_) {
      if (sec_ > fir_) fir_ = sec_;
    }
    void Chkmin(int &fir_, int sec_) {
      if (sec_ < fir_) fir_ = sec_;
    }
    LL QPow(LL x_, LL y_) {
      x_ %= mod;
      y_ %= mod - 1;
      LL ret = 1;
      for (; y_; y_ >>= 1ll) {
        if (y_ & 1) ret = ret * x_ % mod;
        x_ = x_ * x_ % mod;
      }
      return ret;
    }
    LL Inv(LL x_) {
      return QPow(x_, mod - 2);
    }
    LL Sum(LL n_) {
      return ((n_ * (n_ + 1ll)) / 2ll) % (mod - 1);
    }
    void Euler() {
      vis[1] = true, mu[1] = phi[1] = 1; //初值
      for (int i = 2; i <= Lim; ++ i) {
        if (! vis[i]) {
          p[++ p_num] = i;
          mu[i] = -1;
          phi[i] = i - 1;
        }
        for (int j = 1; j <= p_num && i * p[j] <= Lim; ++ j) {
          vis[i * p[j]] = true;
          if (i % p[j] == 0) {
            mu[i * p[j]] = 0;
            phi[i * p[j]] = phi[i] * p[j];
            break;
          }
          mu[i * p[j]] = -mu[i];
          phi[i * p[j]] = phi[i] * (p[j] - 1);
        }
      }
    }
    void Prepare() {
      Euler();
      inv[1] = fac[0] = prodt_phi[0] = prodi_i[0] = 1;
      for (int i = 1; i <= Lim; ++ i) {
        g[0][i] = g[1][i] = 1;
        fac[i] = 1ll * fac[i - 1] * i % mod;
        sumphi[i] = (sumphi[i - 1] + phi[i]) % (mod - 1);
        prodi_i[i] = prodi_i[i - 1] * QPow(i, i) % mod;
        if (i > 1) inv[i] = (mod - mod / i) * inv[mod % i] % mod;
    
        prodt_phi[i] = prodt_phi[i - 1] * QPow(i, phi[i]) % mod;
        inv_prodt_phi[i] = Inv(prodt_phi[i]);
      }
    
      for (int d = 1; d <= Lim; ++ d) {
        for (int j = 1; d * j <= Lim; ++ j) {
          int t = d * j;
          if (mu[j] == 1) {
            g[0][t] = g[0][t] * d % mod;
            g[1][t] = g[1][t] * QPow(1ll * d, 1ll * t * t) % mod;
          } else if (mu[j] == -1) {
            g[0][t] = g[0][t] * inv[d] % mod;
            g[1][t] = g[1][t] * Inv(QPow(1ll * d, 1ll * t * t)) % mod;
          }
        }
      }
      inv_prodg[0][0] = prodg[0][0] = 1;
      inv_prodg[1][0] = prodg[1][0] = 1;
      inv_prodt_phi[0] = 1;
      for (int i = 1; i <= Lim; ++ i) {
        for (int j = 0; j <= 1; ++ j) {
          prodg[j][i] = prodg[j][i - 1] * g[j][i] % mod;
          inv_prodg[j][i] = Inv(prodg[j][i]);
        }
      }
    }
    LL f1(LL a_, LL b_, LL c_, int type) {
      if (! type) return QPow(fac[a_], b_ * c_);
      if (type == 1) return QPow(prodi_i[a_], Sum(b_) * Sum(c_));
      LL ret = 1, lim = min(min(a_, b_), c_);
      for (LL l = 1, r = 1; l <= lim; l = r + 1) {
        r = min(min(a_ / (a_ / l), b_ / (b_ / l)), c_ / (c_ / l));
        ret = ret * QPow(prodt_phi[r] * inv_prodt_phi[l - 1], (a_ / l) * (b_ / l) % (mod - 1) * (c_ / l)) % mod;
        ret = ret * QPow(fac[a_ / l], (sumphi[r] - sumphi[l - 1] + mod - 1) % (mod - 1) * (b_ / l) % (mod - 1) * (c_ / l)) % mod;
      }
      return ret;
    }
    LL f2_2(LL a_, LL b_) { 
      LL ret = 1;
      for (LL l = 1, r = 1; l <= min(a_, b_); l = r + 1) {
        r = min(a_ / (a_ / l), b_ / (b_ / l));
        ret = ret * QPow(prodg[0][r] * inv_prodg[0][l - 1], 1ll * (a_ / l) * (b_ / l)) % mod;
      }
      return ret;
    }
    LL f2(LL a_, LL b_, LL c_, int type) {
      LL ret = 1;
      if (! type) {
        for (LL l = 1, r = 1; l <= min(a_, b_); l = r + 1) {
          r = min(a_ / (a_ / l), b_ / (b_ / l));
          LL val = QPow(prodg[0][r] * inv_prodg[0][l - 1], 1ll * (a_ / l) * (b_ / l));
          ret = (ret * QPow(val, c_)) % mod;
        }
      } else if (type == 1) {
        for (LL l = 1, r = 1; l <= min(a_, b_); l = r + 1) {
          r = min(a_ / (a_ / l), b_ / (b_ / l));
          LL val = QPow(prodg[1][r] * inv_prodg[1][l - 1], Sum(a_ / l) * Sum(b_ / l));
          ret = (ret * QPow(val, Sum(c_))) % mod;
        }
      } else {
        LL lim = min(min(a_, b_), c_);
        for (LL l = 1, r = 1; l <= lim; l = r + 1) {
          r = min(min(a_ / (a_ / l), b_ / (b_ / l)), c_ / (c_ / l));
          ret = ret * QPow(f2_2(a_ / l, b_ / l), (sumphi[r] - sumphi[l - 1] + mod - 1) % (mod - 1) * (c_ / l)) % mod;
          ret = ret * QPow(prodt_phi[r] * inv_prodt_phi[l - 1], (a_ / l) * (b_ / l) % (mod - 1) * (c_ / l)) % mod;
        }
      }
      return ret;
    }
    //=============================================================
    int main() {
      T = read(), mod = read();
      Prepare();
      while (T -- ) {
        A = read(), B = read(), C = read();
        for (int i = 0; i <= 2; ++ i) {
          ans = f1(A, B, C, i) * f1(B, A, C, i) % mod;
          ans = ans * Inv(f2(A, B, C, i)) % mod * Inv(f2(A, C, B, i)) % mod;
          printf("%lld ", ans);  
        }
        printf("
    ");
      }
      return 0;
    }
    

    写在最后

    参考资料:
    Oi-Wiki-莫比乌斯反演
    算法学习笔记(35): 狄利克雷卷积 By: Pecco
    题解 SP5971 【LCMSUM - LCM Sum】 - BJpers2 的博客
    题解 SP5971 【LCMSUM - LCM Sum】 - Venus 的博客
    题解 P3327 【[SDOI2015]约数个数和】 - suncongbo 的博客

    把 Oi-Wiki 上的内容进行了复制 整理扩充。
    我是个没有感情的复读机(大雾)

  • 相关阅读:
    Django(十五)模板详解:模板标签、过滤器、模板注释、模板继承、html转义
    Python @函数装饰器及用法
    NPM概述及使用简介
    MVC、MVT简介
    十、Vue:Vuex实现data(){}内数据多个组件间共享
    九、响应式发:rem和less(适配移动端)
    八、Vue-lazyload
    Vue点到子路由,父级,无法高亮问题解决
    七、Vue组件库:Element、Swiper(轮播专用组件)
    六、Vue-Router:基础路由处理、路由提取成单独文件、路由嵌套、路由传参数、路由高亮、html5的history使用
  • 原文地址:https://www.cnblogs.com/luckyblock/p/12654760.html
Copyright © 2011-2022 走看看