bzoj3309 DZY Loves Math

bzoj3309

---

前戏

∑ai=1∑bj=1f(gcd(i,j))
=∑ad=1∑ai=1∑bj=1e(gcd(i,j)==d)f(d)
=∑ad=1∑a/di=1∑b/dj=1e(gcd(i,j)==1)f(d)
=∑ad=1∑a/di=1∑b/dj=1∑d′|i且d′|jμ(d′)f(d)
=∑ad=1f(d)∑a/dd′=1μ(d′)∑a/dd′i=1∑b/dd′j=11
=∑ad=1f(d)∑a/dd′=1μ(d′)⌊add′⌋⌊bdd′⌋
考虑枚举dd’
=∑ai=1⌊ai⌋⌊bi⌋∑d|if(d)μ(i/d)
很显然f虽然不是积性函数但是可以被筛出来，
对于后面一部分g(i)=∑d|if(d)μ(i/d)我一开始是准备埃氏筛法，但是时限太紧了TAT
然而这玩意儿一副欧拉筛不可做的样子。

解决

然后我们就可以发现我们是不用考虑mu值等于0的那一部分。
考虑i=qr11qr22…qrtt假设其中rw最大且一共有k个r是与rw相等的.
很显然，对于我们考虑的μ(d)，d是{q1,q2…,qn}这个集合的一个子集内所有元素相乘的一个值.
我们可以将这些子集按包不包含元素rw来分成两两一组（也就是说每一组除去rw都应该是一样的），假设除去rw后集合为T如果∃j,rj∉T那么显然他们对应的f值是相等的，也就是说这一对对g(i)的贡献相加=0.
所以我们应该考虑的是所有包含了这k个元素的集合，那么最后的答案应该是
∑t−1i=k−1Ci−(k−1)t−1−(k−1)(−1)i+1∗(ans+1)−Ci−(k−1)t−1−(k−1)∗(−1)i∗ans
=∑t−1i=k−1Ci−(k−1)t−1−(k−1)(−1)i=∑t−ki=0Cit−k(−1)i+k−1
当t=k时，原式=(−1)k+1
当t≠k时，根据二项式定理，原式=(−1+1)t−k=0
所以我们如果能统计出每个数对应的t和k就搞定辣，这个很好筛出来的.
发现标程有辣么短，自己好长QAQ

#include<cstdio>
#include<cmath>
#include<algorithm>
#include<ctime>
const int N = 1e7 + 7;

typedef long long ll;
bool vis[N];
int prime[N], tot, wor[N], wor2[N], sl[N], num[N], jx[N];

void G (int &num) {
    static char a;
    for (a = getchar (); a > '9' || a < '0'; a = getchar ()) ;
    for (num = 0; a >= '0' && a <= '9'; a = getchar ()) num = (num << 3) + (num << 1) + a - '0';
}

void Pre () {
    vis[1] = true;
    for (int i = 2; i <= 10000000; ++i) {
        if (!vis[i]) prime[++tot] = i, wor[i] = wor2[i] = sl[i] = num[i] = 1;
        for (int j = 1; j <= tot && 1ll * i * prime[j] <= 10000000; ++j) {
            vis[i * prime[j]] = true;
            if (i % prime[j] == 0) {
                sl[i * prime[j]] = sl[i] + 1;
                wor[i * prime[j]] = wor[i];
                if (sl[i * prime[j]] == wor2[i]) num[i * prime[j]] = num[i] + 1, wor2[i * prime[j]] = wor2[i]; 
                else if (sl[i * prime[j]] > wor2[i]) num[i * prime[j]] = 1, wor2[i * prime[j]] = sl[i * prime[j]];
                else num[i * prime[j]] = num[i], wor2[i * prime[j]] = wor2[i];
                break;
            }
            sl[i * prime[j]] = 1;
            wor[i * prime[j]] = wor[i] + 1;
            wor2[i * prime[j]] = wor2[i];
            num[i * prime[j]] = num[i];
            if (wor2[i] == 1) ++num[i * prime[j]];
        }
        jx[i] = jx[i - 1];
        if (wor[i] == num[i]) jx[i] += num[i] & 1 ? 1 : -1;
    }
}

ll calc (int a, int b) {
    ll ret = 0;
    for (int i = 1, lo = 1; i <= a; i = lo + 1) {
        lo = std :: min (a / (a / i), b / (b / i));
        ret += 1ll * (jx[lo] - jx[i - 1]) * (a / i) * (b / i);
    }
    return ret;
}

void Solve () {
    int T, a, b;
    G (T);
    while (T--) {
        G (a), G (b);
        if (a > b) a ^= b ^= a ^= b;
        printf ("%lld\n", calc (a, b));
    }
}

void IO () {
    freopen ("3309.in", "r", stdin);
    freopen ("3309.ans", "w", stdout);
}

int main () {
    IO ();
    Pre ();
    Solve ();
    return 0;
}

update :
看了别人的代码，仔细研究了一番，发现我的信息记得过于全面，还是有很多玩意儿是可以省略的，现在附上简洁版的代码

/**************
+1
written by DraZxlNDdt
Accepted
***************/ 
#include<cstdio>
#include<cmath>
#include<algorithm>
#include<ctime>
const int N = 1e7 + 7;

typedef long long ll;
bool vis[N];
int prime[N], tot, sl[N], jx[N], zx[N], g[N];

void G (int &num) {
    static char a;
    for (a = getchar (); a > '9' || a < '0'; a = getchar ()) ;
    for (num = 0; a >= '0' && a <= '9'; a = getchar ()) num = (num << 3) + (num << 1) + a - '0';
}

void Pre () {
    vis[1] = true;
    for (int i = 2; i <= 10000000; ++i) {
        if (!vis[i]) prime[++tot] = i, zx[i] = i, sl[i] = g[i] = 1;
        for (int j = 1; j <= tot && 1ll * i * prime[j] <= 10000000; ++j) {
            vis[i * prime[j]] = true;
            if (i % prime[j] == 0) {
                sl[i * prime[j]] = sl[i] + 1;
                zx[i * prime[j]] = zx[i] * prime[j];
                if (i == zx[i]) g[i * prime[j]] = 1;
                else g[i * prime[j]] = sl[i * prime[j]] == sl[i / zx[i]] ? -g[i / zx[i]] : 0;
                break;
            }
            sl[i * prime[j]] = 1;
            zx[i * prime[j]] = prime[j];
            g[i * prime[j]] = sl[i] == 1 ? -g[i] : 0;
        }
        jx[i] = jx[i - 1];
        if (g[i]) jx[i] += g[i];
    }
}

ll calc (int a, int b) {
    ll ret = 0;
    for (int i = 1, lo = 1; i <= a; i = lo + 1) {
        lo = std :: min (a / (a / i), b / (b / i));
        ret += 1ll * (jx[lo] - jx[i - 1]) * (a / i) * (b / i);
    }
    return ret;
}

void Solve () {
    int T, a, b;
    G (T);
    while (T--) {
        G (a), G (b);
        if (a > b) a ^= b ^= a ^= b;
        printf ("%lld\n", calc (a, b));
    }
}

void IO () {
    freopen ("3309.in", "r", stdin);
    freopen ("3309.ans", "w", stdout);
}

int main () {
//  IO ();
    Pre ();
    Solve ();
    return 0;
}