[51node 1965] 奇怪的式子

Posted 2021-01-27 guessycb

[51node 1965] 奇怪的式子

题目大意

求(prod_{i=1}^n sigma(i)^{mu(i) + i})，数据范围(nleq 10^{11})，答案对(10^{12}+39)取模。

题解

这不显然要拆幂吗？
(Ans = (prod_{i=1}^n sigma(i)^{mu(i)})(prod_{i=1}^n sigma(i)^i))。
然后咋办呢.....
= =
(F_1 = prod_{i=1}^n sigma(i)^{mu(i)})，唔[托腮]......
考虑一下(mu(i) eq 0)的条件，每个质因子只出现一次。
所以(F_1 = prod_{i=1}^n 2^{mu(i)d(i)} = 2^{sum_{i=1}^n mu(i)d(i)})，其中(d(i))为质因子个数。
所以要求(F_1' = sum_{i=1}^n mu(i)d(i))，嗯[托腮$ imes$2]......得出结论，不会......
然后就get到了这个新的套路。
我们考虑把min-25筛的第一部分的那个过程给逆过来。
设(f(n,j) = sum_{i=1}^n [iin prime or p_{min}(i) geq p_j]mu(i)d(i))。
注意到(d(x))的变化是加一变化，所以再设(g(n,j) = sum_{i=1}^n [iin prime or p_{min}(i) geq p_j]mu(i))。
考虑把min-25筛的那个过程倒过来，我们逆推。
初始状态(f(n,|P|+1))、(g(n,|P|+1))我们拿真正的min-25筛给弄出来。
考虑每次加入只包含一个(p_j)的那些数的贡献，有：
(g(n,j) = g(n,j+1) + mu(p_j)(g(lfloorfrac{n}{p_j} floor,j+1) - gpre_{j})) ；
(f(n,j) = f(n,j+1) + mu(p_j)[(f(lfloor frac{n}{p_j} floor,j+1)-fpre_{j}) + (g(lfloor frac{n}{p_j} floor,j+1) -gpre_{j})])
其中(fpre_j = sum_{k=1}^j mu(p_k)d(p_k))，(gpre_j = sum_{k=1}^j mu(p_k))。
注意(gpre)、(fpre)每次减的都是(gpre_{j})而不是(gpre_{j-1})，因为若(p_j)多次出现则(mu=0)不满足积性。
这做法简直NB的不行。
= =
(F_2 = prod_{i=1}^n sigma(i)^{i} = prod_{pin prime} prod_k (k+1)^{delta(n,p,k)})
定义(Sum(n) = sum_{i=1}^n i)
容斥有(delta(n,p,k) = sum_{i=1}^n [p^k | i and i\%p^{k+1} eq 0]i = p^k Sum(lfloor frac{n}{p^k} floor) -p^{k+1} Sum(lfloor frac{n}{p^{k+1}} floor))。
考虑分块。
对于(pleq sqrt{n})，暴力计算即可，复杂度(O(sqrt{n}logn))。
对于(p > sqrt{n})，：
(F_2' = prod_{pin prime , p > sqrt{n}}2^{delta(n,p,1)} = 2^{sum_{pin prime,p>sqrt{n}} delta(n,p,1)})
然后(delta(n,p,1) = pSum(lfloor frac{n}{p} floor))。
注意到(lfloor frac{n}{p} floor)的取值只有(sqrt{n})个，所以考虑先抛弃(p > sqrt{n})这个条件，我们来算：
(F_2'' = sum_{pin prime} pSum(lfloor frac{n}{p} floor) = sum_{d=lfloor frac{n}{x} floor} Sum(d) sum_{pin prime} p[lfloor frac{n}{p} floor = d])
对于(sum_{pin prime} p[lfloor frac{n}{p} floor = d])，显然对应一段区间，而我们知道数论分块后区间右端点封闭。
所以只需要求：(sum_{pin prime} p)，这个就很舒服了，直接min-25筛。
得到(F_2'')后，暴力把(pin prime,pleq sqrt{n})的贡献全部减掉，稍加处理即可得到(F_2)。
= =
最后(Ans = F_1F_2)，注意不要爆(long long)，记得时刻注意模数是(mod)还是(varphi(mod))。

实现代码

卡常警告......
由于卡常数这份代码已经丑陋不堪了，但是这道15sec的题跑了13.4sec还是跑过去了。

#include<bits/stdc++.h>
#define IL inline
#define _ 1000005
#define ll long long
using namespace std ;

ll Num[_],ID1[_],ID2[_],fpre[_],gpre[_],ppre[_],pcnt[_],pval[_],n,m,Base,tot,prm[_],mod ; 
bool Is[_] ; ll f[_],g[_] ;

IL ll mul(ll x , ll y , ll Mod) {
    ll ret = x * y - (long long)((long double)x / Mod * y + 0.5) * Mod ;
    return ret < 0 ? ret + Mod : ret ;
}
IL ll Pow(ll ts , ll js , ll Mod) {
    ll al = 1 ;
    while(js) {
        if(js & 1) al = mul(al , ts , Mod) ;
        ts = mul(ts , ts , Mod) ; js >>= 1 ; 
    }return al ;
}
IL ll ID(ll x) {return x > Base ? ID2[n / x] : ID1[x] ; }
IL ll Sum(ll x , ll Mod) {return (x & 1) ? mul(x , (x + 1) / 2 , Mod) : mul(x / 2 , x + 1 , Mod) ; }

IL void Sieve() {
    Base = sqrt(n) ; tot = 0 ; memset(Is , false , sizeof(Is)) ;
    while(1ll * (Base + 1) * (Base + 1) <= n) ++ Base ; 
    for(int i = 2; i <= Base; i ++) {
        if(!Is[i]) prm[++tot] = i ;
        for(int j = 1; j <= tot && 1ll * prm[j] * i <= Base; j ++) {
            Is[i * prm[j]] = true ; if(i % prm[j] == 0) break ; 
        }
    }
    for(int i = 1; i <= tot; i ++) fpre[i] = gpre[i] = (mod - 1) - i ;
    for(int i = 1; i <= tot; i ++) {
        ppre[i] = ppre[i - 1] + prm[i] ; if(ppre[i] >= (mod - 1)) ppre[i] -= (mod - 1) ;
    }
}
IL void Pre() {
    m = 0 ;
    for(ll l = 1,r; l <= n; l = n / (n / l) + 1) {
        (n / l) <= Base ? ID1[n / l] = ++ m : ID2[n / (n / l)] = ++ m ;
        Num[m] = n / l ; 
    }
    for(int i = 1; i <= m; i ++) pval[i] = (Sum(Num[i] , mod - 1) + mod - 2) % (mod - 1) , pcnt[i] = Num[i] - 1 ;
    for(int i = 1; i <= tot; i ++)
        for(int j = 1; j <= m && 1ll * prm[i] * prm[i] <= Num[j] ; j ++) {
            ll v = ID(Num[j] / prm[i]) , w ;
            w = pcnt[v] - (i-1) + (mod - 1) ; if(w >= (mod - 1)) w -= (mod - 1) ; 
            pcnt[j] = pcnt[j] - w + (mod - 1) ;
            if(pcnt[j] >= mod - 1) pcnt[j] -= mod - 1 ;
            w = pval[v] - ppre[i - 1] + (mod - 1) ; if(w >= (mod - 1)) w -= (mod - 1) ; 
            pval[j] = pval[j] - mul(prm[i] , w , mod - 1) + (mod - 1);
            if(pval[j] >= mod - 1) pval[j] -= mod - 1 ; 
            
        }
    return ; 
}
IL ll Calc1() {
    for(int j = 1; j <= m; j ++) f[j] = g[j] = (mod - 1) - pcnt[j] ;
    for(int j = tot; j >= 1; j --)
        for(int i = 1; i <= m && 1ll * prm[j] * prm[j] <= Num[i] ; i ++) {
            ll v = ID(Num[i] / prm[j]) , w ;
            w = g[v] - gpre[j] + (mod - 1) ; if(w >= (mod - 1)) w -= (mod - 1) ;
            g[i] = g[i] - w + (mod - 1) ; if(g[i] >= mod - 1) g[i] -= mod - 1 ;
            w = f[v] - fpre[j] + (mod - 1) ; if(w >= (mod - 1)) w -= (mod - 1) ;  
            f[i] = f[i] - w + (mod - 1) ; if(f[i] >= mod - 1) f[i] -= mod - 1 ;
            w = g[v] - gpre[j] + (mod - 1) ; if(w >= (mod - 1)) w -= (mod - 1) ;
            f[i] = f[i] - w + (mod - 1) ; if(f[i] >= mod - 1) f[i] -= mod - 1 ;
        }
    return Pow(2ll , f[1] , mod) ;
}

IL ll slv(ll p1,ll p2){
    ll w = (mul(p1 , Sum(n/p1,mod-1) , mod-1) - mul(p2 , Sum(n/p2,mod-1) , mod-1) + (mod-1)) + (mod - 1) ;
    return w >= mod - 1 ? w - (mod - 1) : w ; 
}
IL ll Calc2() {
    ll ans = 1 ;
    for(ll j = 1; j <= tot; j ++)
        for(ll pk = prm[j],k = 1 ; pk <= n; pk = 1ll * pk * prm[j] , ++ k) 
            ans = mul(ans , Pow(k + 1 , slv(pk , pk * prm[j]) , mod) , mod) ;
    ll ret = 0 ;
    for(ll i = m,w; i >= 1; i --) {
        w = pval[i] - pval[i+1] + (mod-1) ; if(w >= (mod - 1)) w -= (mod - 1) ; 
        ret = (ret + mul(Sum(n / Num[i] , mod - 1) , w , mod - 1)) % (mod - 1) ;
    }
    for(ll i = 1; i <= tot; i ++) {
        ret = ret - mul(Sum(n / prm[i] , mod - 1) , prm[i] , mod - 1) + (mod - 1) ;
        if(ret >= (mod - 1)) ret -= (mod - 1) ; 
    }
    return mul(ans , Pow(2ll , ret , mod) , mod) ; 
}

IL void Work() {Sieve() ; Pre() ; cout << mul(Calc1() , Calc2() , mod) << endl ; }
int main() {
    //freopen("testdata.in","r",stdin) ;
    int Case ; cin >> Case ;
    mod = 1e12 ; mod += 39 ; 
    while(Case --) cin >> n , Work() ; return 0 ;
}