Project Euler 457 题解
初等数论小题目
求
\[n^2-3n-1\equiv 0\pmod {p^2}
\]
配方,得到:
\[(2n-3)^2\equiv 13 \pmod {p^2}
\]
根据亨泽尔引理,只需得到 \((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;
}
YJX AK IOI