PE 216 【二次同余】

作为一个AFO的选手我又回到了机房（￣口￣）
能继续学OI当然很开心，但是看到机房里的人少了一半，连含爷也走了，外带坑了我妈两万软妹币,结果还是只能狗牌退出，还是好难过。。。各种纠结啊有木有
算了，做道数学题冷静一下（说好的数据结构呢...

———–正题—————–

题目大意：设 t(n)=2n2−1,求t(n)中素数的个数2<=n<=5∗107

这道题一看就觉得可以用Miller Rabin水过去，结果愉快地打完后一直WA。。。
没办法只有好好想。。。因为我弱嘛，所以想的办法当然也渣啦╮(╯▽╰)╭

考虑两个位置 i，j 和一个素数 p，其中 p|t(i),现在要找p|t(j),
因为 t(i)−t(j)=2(i−j)(i+j),p 是素数，所以只能是p|i−j或p|i+j,这样就能找到所有 p 的倍数
现在对于每个质数 p，要找到一个 i 使p|t(i)，即

i2=2−1(modp)

然后解个二次剩余。。。

神犇们肯定觉得很水对不对，但是我这种弱渣不会解二次剩余就GG啦
今天就搞了这一道题，果真太弱QAQ
话说今天整个机房都好浪→→（怎么又扯远了。。。不过反正也没人看。。。

二次剩余：

x2≡n(modp)(p为奇素数)

方程有解，iff np−12≡1(modp),n 叫做模 p 的平方剩余（此时有两解）
方程无解，iff np−12≡−1(modp),n 叫做模 p 的平方非剩余
在模 p 的简化系中，平方剩余和平方非剩余各有p−12个
设 w=a2−n,且x2≡w(modp)无解，那么x=(a+w−−√)p+12是x2≡n(modp)的解

写个题解都写了这么久。。。果真我太弱QAQ
废话好多QAQ

【答案】5437849

#include<iostream>#include<algorithm>#include<string>#include<cstring>#include<ctime>#include<cmath>#include<cstdlib>#include<cstdio>#define N 70710679#define LL long longusing namespace std;LL n,w,tot,cnt;LL a[N];int num[N],prime[N];bool v[N];struct T{LL p,d;};void get_prime(){    for (int i = 2;i < N;i ++)    {        if (!v[i]) prime[tot ++] = i;        for (int j = 0;j < tot && i * prime[j] < N;j ++)        {            v[i * prime[j]] = true;            if (i % prime[j] == 0) break;        }    }}LL ksm(LL a,LL b,LL c){    LL ret = 1;    for (;b;b >>= 1,a = a * a % c)        if (b & 1) ret = ret * a % c;    return ret;}T mul(T a,T b,LL c){    return (T){(a.p * b.p + a.d * b.d % c * w) % c,(a.p * b.d + b.p * a.d) % c};}T T_ksm(T a,LL b,LL c){    T ret = (T){1,0};    for (;b;b >>= 1,a = mul(a,a,c))        if (b & 1) ret = mul(ret,a,c);    return ret;}bool check(LL a,LL p){    return ksm(a,p-1 >> 1,p) == 1;}LL tr(LL n,LL p){    if (!check(n,p)) return -1;    LL a;    while (true)    {        a = rand() % p;        w = ((a * a - n) % p + p) % p;        if (!check(w,p)) break;    }    T ret = T_ksm((T){a,1},p+1 >> 1,p);    return ret.p;}void work(int p,int k){    int x = - p % k + k;    for (;x <= n;x += k)        for (;a[x] % k == 0;a[x] /= k) num[x] ++;    for (x = p % k;x <= n;x += k)        for (;a[x] % k == 0;a[x] /= k) num[x] ++;}int main(){    cin >> n;    get_prime();    for (LL i = 0;i <= n;i ++) a[i] = i * i * 2 - 1;    for (int i = 1;i < tot;i ++)    {        int x = tr(ksm(2,prime[i]-2,prime[i]),prime[i]);        if (~ x) work(x,prime[i]);    }    for (int i = 2;i <= n;i ++)        if (a[i] ^ 1) num[i] ++;    for (int i = 2;i <= n;i ++)        if (num[i] == 1) cnt ++;    cout << cnt << endl;    return 0;}