【51nod】1602 矩阵方程的解

这个行向量显然就是莫比乌斯函数啦，好蠢的隐藏方法= =

然后我们尝试二分，二分的话要求一个这个东西

\(H(n) = \sum_{i = 1}^{n} \mu(i) == d\)

当然\(\mu(x)\)由于一些很好的性质，这个东西可以用分类讨论做出来

众所周知，求\(\mu\)不为0的数的方法就是容斥求无平方因子数

\(G(n) = \sum_{i = 1}^{\sqrt{N}} \mu(i) \lfloor \frac{N}{i^{2}} \rfloor\)

这样我们已经可以快速求得\(\mu(i) = 0\)的位置有多少了，就是\(N - G(n)\)

再求一个\(F(n)\)作为莫比乌斯函数的前缀和，也就是1和-1的差值

可以用\(F(n)\)和\(G(n)\)和一点特判快速求出来1和-1的个数，显然\(F(n)\)需要的复杂度远大于\(G(n)\)，是\(O(n^{2/3})\)的

梦想选手快乐的写了一发，只过了d = 0的点，惨兮兮

然后经过努力的尝试，发现如果\(R - L < 10^{6}\)的时候停止二分，转而改为通过类似埃氏筛的方法暴力求这\(10^{6}\)个位置的mu函数，然后挨个枚举看看哪里是答案

这样的话可以在3s左右通过了

#include <bits/stdc++.h>

#define fi first

#define se second

#define pii pair<int,int>

#define mp make_pair

#define pb push_back

#define space putchar(' ')

#define enter putchar('\n')

#define eps 1e-10

#define ba 47

#define MAXN 1000005

//#define ivorysi

using namespace std;

typedef long long int64;

typedef unsigned int u32;

typedef double db;

template<class T>

    void read(T &res) {

    res = 0;T f = 1;char c = getchar();

    while(c < '0' || c > '9') {

        if(c == '-') f = -1;

        c = getchar();

    }

    while(c >= '0' && c <= '9') {

        res = res * 10 +c - '0';

        c = getchar();

    }

    res *= f;

}

template<class T>

    void out(T x) {

    if(x < 0) {x = -x;putchar('-');}

    if(x >= 10) {

        out(x / 10);

    }

    putchar('0' + x % 10);

}

int prime[10000005],tot,mu[10000005],M[10000005];

int c[10000005];

int64 val[10000005];

bool nonprime[10000005];

struct hash {

    struct node {

        int64 x,val;int nxt;

    }E[10000005];

    int sumE,head[1000000],mo = 974711;

    void add(int64 x,int64 v) {

        int u = x % mo;

        E[++sumE].x = x;

        E[sumE].val = v;

        E[sumE].nxt = head[u];

        head[u] = sumE;

    }

    int64 Query(int64 x) {

        int u = x % mo;

        for(int i = head[u] ; i ; i = E[i].nxt) {

            if(E[i].x == x) return E[i].val;

        }

        return -1e18;

    }

}hsh;

map<int64,int64> zz;

int64 F(int64 n) {

    if(n <= 10000000) return M[n];

    //if(zz.count(n)) return zz[n];

    int64 t = hsh.Query(n);

    if(t > -1e18) return t;

    int64 res = 1;

    for(int64 i = 2 ; i <= n ; ++i) {

        int64 r = n / (n / i);

        res -= F(n / i) * (r - i + 1);

        i = r;

    }

    //zz[n] =res;

    hsh.add(n,res);

    return res;

}

int64 G(int64 n) {

    int64 ans = 0;

    for(int64 i = 1 ; i <= n / i ; ++i) {

        ans += 1LL * mu[i] * (n / (i * i));

    }

    return ans;

}

int getmu(int64 x) {

    int r = 1;

    for(int i = 1 ; prime[i] <= x / prime[i] ; ++i) {

        if(x % prime[i] == 0) {

            if(x % (prime[i] * prime[i]) == 0) return 0;

            r = -r;

            x /= prime[i];

        }

    }

    if(x != 1) r = -r;

    return r;

}

int64 check(int64 n,int d) {

    int64 g = G(n);

    if(d == 0) return n - g;

    int64 f = F(n);

    int64 t = (g - abs(f)) / 2;

    if(d == -1 && f < 0) t += abs(f);

    if(d == 1 && f > 0) t += f;

    return t;

}

int main(){

    #ifdef ivorysi

    freopen("f1.in","r",stdin);

    #endif

    mu[1] = 1;M[1] = 1;

    for(int i = 2 ; i <= 10000000 ; ++i) {

        if(!nonprime[i]) {

            prime[++tot] = i;

            mu[i] = -1;

        }

        for(int j = 1 ; j <= tot ; ++j) {

            if(prime[j] > 10000000 / i) break;

            nonprime[i * prime[j]] = 1;

            if(i % prime[j] == 0) break;

            else mu[i * prime[j]] = -mu[i];

        }

        M[i] = M[i - 1] + mu[i];

    }

    int64 d,k;

    read(d);read(k);

    int64 L = 1,R = 1e10;

    while(R - L > 1e6) {

        int64 mid = (L + R) >> 1;

        int64 x = check(mid,d);

        if(x >= k) R = mid - 2 * (x - k);

        else L = mid + 2 * (k - x);

    }

    for(int i = 0 ; i <= R - L ; ++i) {c[i] = 1;val[i] = i + L;}

    for(int i = 1 ; i <= tot ; ++i) {

        if(1LL * prime[i] * prime[i] > R) break;

        int64 l = ((L - 1)/ prime[i] + 1) * prime[i];

        int64 r = (R / prime[i]) * prime[i];

        for(int64 j = l - L ; j <= r - L ; j += prime[i]) {

            if(!c[j]) continue;

            if(val[j] % (1LL * prime[i] * prime[i]) == 0) c[j] = 0;

            else {c[j] = -c[j];val[j] /= prime[i];}

        }

    }

    int64 x = check(L - 1,d);

    for(int i = 0 ; i <= R - L ; ++i) {

        if(c[i] && val[i] != 1) c[i] = -c[i];

        if(c[i] == d) ++x;

        if(x == k) {out(i + L);enter;break;}

    }

    //printf("%.3lf\n",1.0 * clock() / CLOCKS_PER_SEC);

    return 0;

}

巴特西

【51nod】1602 矩阵方程的解

【51nod】1602 矩阵方程的解

最新文章

热门文章