【UOJ #221】【NOI 2016】循環之美

http://uoj.ac/problem/221
由於\(a\)和\(b\)不互質時，\(\frac ab=\frac{\frac a{(a,b)}}{\frac b{(a,b)}}\)，因此只用求\(a\)和\(b\)互質時的知足條件的個數。
\(\frac ab\)在\(k\)進制下是純循環小數，咱們先假設循環節長度爲\(l\)，這樣\(\frac ab\times k^l-\frac ab\)的小數部分就是0，也就是這是個整數。
\(\frac{a\left(k^l-1\right)}b\)是個整數，就是說\(b|a\left(k^l-1\right)\)。又由於\(a\bot b\)，因此\(b|\left(k^l-1\right)\)。
判斷\(\frac ab\)是不是\(k\)進制下的純循環小數就轉變成了判斷是否存在一個\(l,l\geq 0\)，知足\(\left(k^l-1\right)\mod b=0\)。
當\(k\)和\(b\)不互質時，對於全部的\(l\)，\(b\)和\(k^l\)有共同的質因子，又由於\(k^l\bot\left(k^l-1\right)\)，因此\(b\)和\(k^l\)的共同的質因子是\(k^l-1\)沒有的，因此不存在\(l\)知足\(b|\left(k^l-1\right)\)。
當\(k\)和\(b\)互質時，由歐拉定理：\(k^{\varphi(b)}\mod b=1\)，存在\(l=\varphi(b)\)。
這樣答案就變成了：\[\sum_{a=1}^n\sum_{b=1}^m[b\bot k][a\bot b]\]
化一波式子：
\[ \begin{aligned} &\sum_{a=1}^n\sum_{b=1}^m[b\bot k][a\bot b]\\ =&\sum_{d=1}^n\mu(d)\sum_{a=1}^{\left\lfloor\frac nd\right\rfloor}\sum_{b=1}^{\left\lfloor\frac md\right\rfloor}[bd\bot k]\\ =&\sum_{d=1}^n\mu(d)[d\bot k]\left\lfloor\frac nd\right\rfloor\sum_{b=1}^{\left\lfloor\frac md\right\rfloor}[b\bot k] \end{aligned} \]
\(\sum\limits_{i=1}^n[i\bot k]\)很好預處理後\(O(1)\)計算。
重點是怎麼算\(\sum\limits_{d=1}^n\mu(d)[d\bot k]\)？
先放寬限制，把\([d\bot k]\)的限制去掉，這樣就是對\(O\left(\sqrt n+\sqrt m\right)\)個下取整取值求\(\mu\)的前綴和，能夠先\(O\left(n^{\frac 23}\right)\)大力杜教篩一波。
而後考慮\(S(i,n)\)表示\(1\sim n\)中與\(k\)的前\(i\)個質因子互質的數的\(\mu\)值和，這樣\(S(i,n)=S(i-1,n)-\mu\left(p_i\right)S\left(i,\left\lfloor\frac n{p_i}\right\rfloor\right)\)遞推求出，時間複雜度\(O\left(\left(\sqrt n+\sqrt m\right)\log k\right)\)（假設\(k\)的不一樣質因子有\(\log k\)個）。
最後對\(O\left(\sqrt n+\sqrt m\right)\)個不一樣的\(\left\lfloor\frac nd\right\rfloor\)和\(\left\lfloor\frac md\right\rfloor\)的取值進行分段求和便可。
時間複雜度\(O\left(n^{\frac 23}+\left(\sqrt n+\sqrt m\right)\log k\right)\)。

#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
typedef long long ll;

const int N = 1000003;

bool notp[N];
int mu[N], sum_mu[N], prime[N], Num = 0;

void Euler_shai() {
    mu[1] = sum_mu[1] = 1;
    for (int i = 2; i <= 1000000; ++i) {
        if (!notp[i]) prime[++Num] = i, mu[i] = -1;
        for (int j = 1; j <= Num && prime[j] * i <= 1000000; ++j) {
            notp[prime[j] * i] = true;
            if (i % prime[j] == 0) break;
            mu[prime[j] * i] = -mu[i];
        }
        sum_mu[i] = sum_mu[i - 1] + mu[i];
    }
}

int R[N], cnt = 0, sum_mu2[N], n, m, k, S[N], pr[N], prnum = 0;

const int mo = 2333333;

struct HashTable {
    int pos[mo], num[mo];
    
    void ins(int nu, int po) {
        int tmp = nu % mo;
        while (num[tmp]) {++tmp; if (tmp == mo) tmp = 0;}
        num[tmp] = nu;
        pos[tmp] = po;
    }
    
    int query(int nu) {
        int tmp = nu % mo;
        while (num[tmp] != nu) {++tmp; if (tmp == mo) tmp = 0;}
        return pos[tmp];
    }
} HT;

int sum[N], id[N];

int gcd(int a, int b) {return b ? gcd(b, a % b) : a;}

void pre_sum() {
    for (int i = 1; i <= k; ++i) {
        sum[i] = sum[i - 1];
        if (gcd(i, k) == 1)
            ++sum[i];
    }
}

int cal(int num) {return sum[k] * (num / k) + sum[num % k];}

int Sum(int num) {return num <= 1000000 ? S[id[num]] : S[HT.query(num)];}

int main() {
    scanf("%d%d%d", &n, &m, &k);
    Euler_shai();
    
    for (int i = 1, j = 1; i <= n && i <= m; ++i) {
        j = m / (m / i);
        i = n / (n / i);
        if (i > j) i = j;
        R[++cnt] = i;
        if (i > 1000000) HT.ins(i, cnt);
    }
    
    for (int i = 1; i <= cnt; ++i) {
        int num = R[i], &ret = S[i];
        if (num <= 1000000) {ret = sum_mu[num]; id[num] = i; continue;}
        ret = 1;
        for (int j = 2, pre = 1; j <= num; pre = j, ++j) {
            j = num / (num / j);
            ret -= 1ll * Sum(num / j) * (j - pre);
        }
    }
    
    for (int i = 1; i <= Num && prime[i] <= k; ++i)
        if (k % prime[i] == 0)
            pr[++prnum] = prime[i];
    
    for (int i = 1, pi = pr[1]; i <= prnum; pi = pr[++i])
        for (int j = 1; j <= cnt; ++j)
            S[j] += Sum(R[j] / pi);
    
    pre_sum();
    
    ll ans = 0;
    for (int tmp = 1, i = 1, j = 1; i <= n && i <= m; ++tmp, ++i) {
        j = m / (m / i);
        i = n / (n / i);
        if (i > j) i = j;
        ans += 1ll * (S[tmp] - S[tmp - 1]) * (n / i) * cal(m / i);
    }
    
    printf("%lld\n", ans);
    return 0;
}