bzoj3309 DZY Loves Math

bzoj3309

---

前戏

$\sum_{i=1}^{a}\sum_{j=1}^{b}f(gcd(i,j))$
$= \sum_{d=1}^{a}\sum_{i=1}^{a}\sum_{j=1}^{b}e(gcd(i,j)==d)f(d)$
$= \sum_{d=1}^{a}\sum_{i=1}^{a/d}\sum_{j=1}^{b/d}e(gcd(i,j)==1)f(d)$
$= \sum_{d=1}^{a}\sum_{i=1}^{a/d}\sum_{j=1}^{b/d}\sum_{d'|i且d'|j}\mu(d')f(d)$
$= \sum_{d=1}^{a}f(d)\sum_{d'=1}^{a/d}\mu(d')\sum_{i=1}^{a/dd'}\sum_{j=1}^{b/dd'}1$
$= \sum_{d=1}^{a}f(d)\sum_{d'=1}^{a/d}\mu(d') \left \lfloor \frac {a} {dd'} \right \rfloor\left \lfloor \frac {b} {dd'} \right \rfloor$
考虑枚举dd’
$= \sum_{i=1}^{a} \left \lfloor \frac {a} {i} \right \rfloor\left \lfloor \frac {b} {i} \right \rfloor \sum_{d|i} f(d)\mu(i/d)$
很显然f虽然不是积性函数但是可以被筛出来，
对于后面一部分$g(i)=\sum_{d|i} f(d)\mu(i/d)$我一开始是准备埃氏筛法，但是时限太紧了TAT
然而这玩意儿一副欧拉筛不可做的样子。

解决

然后我们就可以发现我们是不用考虑mu值等于0的那一部分。
考虑$i=q_{1}^{r_{1}}q_{2}^{r_{2}} \dots q_{t}^{r_{t}}$假设其中$r_{w}$最大且一共有k个r是与$r_{w}$相等的.
很显然，对于我们考虑的$\mu(d)$，d是$\{ q_{1},q_{2} \dots , q_{n} \}$这个集合的一个子集内所有元素相乘的一个值.
我们可以将这些子集按包不包含元素${r_{w}}$来分成两两一组（也就是说每一组除去$r_{w}$都应该是一样的），假设除去$r_{w}$后集合为$T$如果$\exists j, r_{j} \notin T$那么显然他们对应的f值是相等的，也就是说这一对对g(i)的贡献相加=0.
所以我们应该考虑的是所有包含了这k个元素的集合，那么最后的答案应该是
$\sum_{i=k-1}^{t-1}C_{t-1-(k-1)}^{i-(k-1)}(-1)^{i+1}*(ans+1)- C_{t-1-(k-1)}^{i-(k-1)}*(-1)^i*ans$
=$\sum_{i=k-1}^{t-1}C_{t-1-(k-1)}^{i-(k-1)}(-1)^{i}=\sum_{i=0}^{t-k}C_{t-k}^{i}(-1)^{i+k-1}$
当$t=k$时，原式=$(-1)^{k+1}$
当$t \neq 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;
}