Submission #3027726

---

Source Code Expand

#include <cstdio>
#include <cstring>
#include <iostream>
#include <string>
#include <cmath>
#include <bitset>
#include <vector>
#include <map>
#include <set>
#include <queue>
#include <deque>
#include <algorithm>
#include <complex>
#include <unordered_map>

using namespace std;
typedef long long int ll;
typedef pair<int, int> P;

const int MAX=5000051;

int l, s;

vector<int> prime;
bool isprime[MAX];

ll gcd(ll a, ll b){
	if(b==0) return a;
	return gcd(b, a%b);
}

ll powmod(ll a, ll k, ll m){
    ll ap=a, ans=1;
    while(k>0){
        if(k%2==1){
            ans*=ap;
            ans%=m;
        }
        ap=ap*ap;
        ap%=m;
        k/=2;
    }
    return ans;
}

ll modsqrt(int p){
	int q=p-1, s=0;
	while(q%2==0){
		q/=2;
		s++;
	}
	int z=2;
	while(1){
		if(powmod(z, (p-1)/2, p)==p-1) break;
		z++;
	}
	ll c=powmod(z, q, p);
	int r;
	if((q+1)/2%2==0) r=1;
	else r=p-1;
	for(int i=0; i<s-2; i++){
		c=c*c%(ll)p;
	}
	return (ll)r*c%(ll)p;
}

P factor2(int p){
	if(p==2) return P(1, 1);
	if(p%4==3) return P(-1, -1);
	ll x=modsqrt(p);
	if(x>p/2) x=p-x;
	ll m=(x*x+1)/(ll)p;
	ll x1=x, y1=1;
	while(m>1){
		ll x2=x1%m, y2=y1%m;
		if(x2>m/2) x2-=m;
		if(y2>m/2) y2-=m;
		ll r=(x2*x2+y2*y2)/m;
		ll x3=(x1*x2+y1*y2)/m, y3=(x1*y2-x2*y1)/m;
		m=r;
		if(x3<0) x1=-x3;
		else x1=x3;
		if(y3<0) y1=-y3;
		else y1=y3;
	}
	return P(x1, y1);
}

P f2[MAX];

void sieve(){
	for(ll i=3; i<MAX; i+=2){
		isprime[i]=1;
	}
	isprime[2]=1;
	prime.push_back(2);
	for(ll i=3; i<MAX; i++){
		if(isprime[i]){
			prime.push_back(i);
          if(i%4==1){
            f2[i]=factor2(i);
          }
			for(ll j=2*i; j<MAX; j+=i){
				isprime[j]=0;
			}
		}
	}
	return;
}

int ans0;

vector<P> factors[10001], factor[MAX/5];

ll x3[MAX/5], y3[MAX/5], x3s[10001], y3s[10001];

bool muri[MAX];

int ansm[MAX];

void fact(){
	for(auto p:prime){
      if(p==2 && s%2==0) continue;
      int c=min(s%p, (p-s%p)%p), c1=max(s%p, (p-s%p)%p);
      bool reigai=0;
      if(c==c1) reigai=1;
      for(; c<=l/2; ){
        if(muri[c] || c<=s || c%2==0){
          if(reigai){
            c+=p;
          }else{
            int c2=c1;
            c1=c+p;
            c=c2;
          }
          continue;
        }
        int e=0, c0;
        if((c-s)%p==0){
          c0=c-s;
          while(c0%p==0){
            c0/=p; e++;
          }
        }
        if((c+s)%p==0){
          c0=c+s;
          while(c0%p==0){
            c0/=p; e++;
          }
        }
        if(p%4==3 && e%2==1){
          muri[c]=1;
          if(reigai){
            c+=p;
          }else{
            int c2=c1;
            c1=c+p;
            c=c2;
          }
          continue;
        }
		if(p%4==3){
			if(c<=10000){
				for(int i=0; i<e/2; i++) x3s[c]*=p, y3s[c]*=p;
			}else if(c>=2*l/5){
				for(int i=0; i<e/2; i++) x3[c-2*l/5]*=p, y3[c-2*l/5]*=p;
			}
		}else if(p==2){
			if(c<=10000){
				for(int i=0; i<e; i++){
					ll a1=x3s[c]-y3s[c], b1=x3s[c]+y3s[c];
					x3s[c]=a1, y3s[c]=b1;
				}
			}else if(c>=2*l/5){
				for(int i=0; i<e; i++){
					ll a1=x3[c-2*l/5]-y3[c-2*l/5], b1=x3[c-2*l/5]+y3[c-2*l/5];
					x3[c-2*l/5]=a1, y3[c-2*l/5]=b1;
				}
			}
		}else{
          if(c<=10000){
				factors[c].push_back(P(p, e));
			}else if(c>=2*l/5){
				factor[c-2*l/5].push_back(P(p, e));
			}else{
				if(s%p!=0){
					ansm[c]*=(e+1);
				}
			}
		}
        if(reigai){
          c+=p;
        }else{
          int c2=c1;
          c1=c+p;
          c=c2;
        }
      }
    }
}

ll a, b;
int ans;

void solve(int i, ll c){
	if((c<=10000 && i==factors[c].size()) || (c>=2*l/5 && i==factor[c-2*l/5].size())){
      if(gcd(gcd(abs(a), abs(b)), c)>1) return;
		if(a>0 && b>0 && a<=b){
			if(a+b+c<=l && a+b>c && b<=c) ans++;
		}
		if(a>0 && b<0 && -b<=a){
			if(a-b+c<=l && a-b>c && a<=c) ans++;
		}
		if(a<0 && b>0 && b<=-a){
			if(-a+b+c<=l && -a+b>c && -a<=c) ans++;
		}
		if(a<0 && b<0 && -a<=-b){
			if(-a-b+c<=l && -a-b>c && -b<=c) ans++;
		}
		return;
	}
	int p, e;
	if(c<=10000){
		p=factors[c][i].first, e=factors[c][i].second;
	}else if(c>=2*l/5){
		p=factor[c-2*l/5][i].first, e=factor[c-2*l/5][i].second;
	}
	ll x1=1, y1=0;
	ll a1=a, b1=b;
	for(int j=0; j<e/2; j++){
		x1*=p;
	}
	if(e%2){
		ll x11=x1*f2[p].first-y1*f2[p].second, y11=y1*f2[p].first+x1*f2[p].second;
		x1=x11, y1=y11;
	}
	a=a1*x1-b1*y1, b=a1*y1+b1*x1;
	solve(i+1, c);
  if(y1!=0){
	a=a1*x1+b1*y1, b=-a1*y1+b1*x1;
	solve(i+1, c);
  }
	int j=2+e%2;
	for(; j<=e; j+=2){
		x1/=p, y1/=p;
		ll x11=x1*f2[p].first-y1*f2[p].second, y11=y1*f2[p].first+x1*f2[p].second;
		x1=x11, y1=y11;
		x11=x1*f2[p].first-y1*f2[p].second, y11=y1*f2[p].first+x1*f2[p].second;
		x1=x11, y1=y11;
		a=a1*x1-b1*y1, b=a1*y1+b1*x1;
		solve(i+1, c);
		a=a1*x1+b1*y1, b=-a1*y1+b1*x1;
		solve(i+1, c);
	}
}

int main()
{
	cin>>l>>s;
  if(l>1e7) return 0;
  sieve();
	for(int i=1; i<=10000; i++) x3s[i]=1;
	for(int i=2*l/5; i<=l/2; i++) x3[i-2*l/5]=1;
	for(int i=10000; i<=2*l/5; i++) ansm[i]=1;
  fact();
	for(ll c=s+3-(s%2); c<=l/2; c+=2){
      if(muri[c]) continue;
		if(c<=10000){
			a=x3s[c], b=y3s[c];
			ans=0;
			solve(0, c);
			ans0+=ans;
		}else if(c>=2*l/5){
			a=x3[c-2*l/5], b=y3[c-2*l/5];
			ans=0;
			solve(0, c);
			ans0+=ans;
		}else{
			ans0+=((ansm[c]+ansm[c]%2)/2);
		}
	}
  cout<<ans0<<endl;
	return 0;
}

Submission Info

Judge Result