【51nod】1602 矩阵方程的解
【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;
}
