Project Euler 457 题解

初等数论小题目

求

配方，得到：

根据亨泽尔引理，只需得到 \((2n-3)^2\equiv 13 \pmod {p}\) 的解即可提升到 \(p^2\)。这是二次剩余，直接解。

单次求解 \(O(\log n)\)，时间复杂度 \(O(n)\)。

#include<bits/stdc++.h>
using namespace std;
#define int long long
namespace Cp{
	int mod,I;
	struct comp{
		int x,y;
		comp(int a=0,int b=0){
			x=a,y=b;
		}
	};
	bool operator ==(comp a,comp b){
		return a.x==b.x&&a.y==b.y;
	}
	comp operator *(comp a,comp b){
		return comp((a.x*b.x+I*a.y%mod*b.y)%mod,(a.x*b.y+a.y*b.x)%mod);
	}
	comp qp(comp a,int b){
		if(b==0)return comp(1);
		comp T=qp(a,b>>1);T=T*T;
		if(b&1)return T*a;
		return T;
	}
	bool ck(int x){
		return qp(x,(mod-1)/2)==comp(1,0);
	}
	mt19937 rng(time(0));
	void solve(int n,int p,int &x1,int &x2){
		int a=rng()%mod;
		while(!a||ck((a*a+mod-n)%mod))a=rng()%mod;
		I=(a*a+mod-n)%mod;
		x1=(qp(comp(a,1),(mod+1)>>1).x)%mod;
		x2=mod-x1;
	}
	pair<int,int> solve(int N,int P){
		mod=P;
		if(!ck(N))return {-1,-1};
		int x1,x2;solve(N,P,x1,x2);
		return {x1,x2};
	}
}
int baoli(int p){
	int mod=p*p;
	for(int i=0;i<p*p;i++)if((i*i-3*i+3*mod)%mod==1)return i;
	return 0;
}
int qp(int a,int b,int p){
	if(b==0)return 1;
	int T=qp(a,b>>1,p);T=T*T%p;
	if(b&1)return T*a%p;
	return T;
}
int hez(int x,int p){
	int M=p*p;
	if((8*x-12%p+p)%p!=0){
		int t=-qp((8*x-12%p+p)%p,p-2,p)*(((4*x*x-12*x-4))/p)%p;
		t=(t%p+p)%p;
		return t*p+x;
	}else if((4*x*x%M-12*x%M+9+M)%M==13)return x;
	return p*p;
}
int calc(int p){
	if(p==2)return 0;
	auto [a,b]=Cp::solve(13,p);
	if(a==-1)return 0;
	int I2=qp(2,p-2,p);a=(a+3)*I2%p,b=(b+3)*I2%p;
	int res=min(hez(a,p),hez(b,p));
	if(res>=p*p)return 0;
	return res;
}
const int maxn=1e7+5;
bool isp[maxn];
vector<signed> pr;int ans=0;
signed n;
signed main(){
	ios::sync_with_stdio(false);
	cin.tie(0);cout.tie(0);
	cin>>n;
	for(int i=2;i<=n;i++){
		if(!isp[i])ans+=calc(i),pr.push_back(i);
		for(auto u:pr){
			if(i*u>n)break;
			isp[i*u]=1;
			if(i%u==0)break;
		}
	}
	cout<<ans<<endl;
	
	return 0;
}