描述

求解關於xx的方程：
$x^2\equiv a(\mathrm{m}\mathrm{o}$

輸入

第一行一個數，表示有T個方程
以下T行，每行兩個數表示a,p

輸出

T行，每行一個數，表示方程的解。無解輸出-1

樣例輸入

2
1 2
13 27

樣例輸出

1
16

提示

保證所有p≤1e15，T≤5000，不保證a<p

解析：

去另外一篇部落格上看：https://blog.csdn.net/zxyoi_dreamer/article/details/85195819

程式碼：

#include<bits/stdc++.h>
using namespace std;
#define ll long long 
#define re register
#define gc getchar
#define pc putchar
#define puts put_s
#define cs const

namespace IO{
	namespace
 
 READONLY{
		cs int Rlen=1<<18|1;
		char buf[Rlen],*p1,*p2;
		char obuf[Rlen],*p3=obuf;
		char ch[23];
	}
	inline char get_char(){
		using namespace READONLY;
		return (p1==p2)&&(p2=(p1=buf)+fread(buf,1,Rlen,stdin),p1==p2)?EOF:*p1++;
	}
	inline void put_char(cs char &c){
		using namespace READONLY;
		*p3++=c;
		if(p3==obuf+Rlen)fwrite(obuf,1,Rlen,stdout),p3=obuf;
	}
	inline void put_s(cs char *s){
		for(;*s;++s)pc(*s);
		pc('\n');
	}
	inline void FLUSH(){
		using namespace READONLY;
		fwrite(obuf,1,p3-obuf,stdout);
		p3=obuf;
	}
	
	inline ll getint(){
		re ll num;
		re char c;
		while(!isdigit(c=gc()));num=c^48;
		while(isdigit(c=gc()))num=(num+(num<<2)<<1)+(c^48);
		return num;
	}
	inline void outint(ll a){
		using namespace READONLY;
		if(a==0)pc('0');
		if(a<0)pc('-'),a=-a;
		while(a)ch[++ch[0]]=a-a/10*10,a/=10;
		while(ch[0])pc(ch[ch[0]--]^48);
	}
}
using namespace IO;

namespace Linear_sieves{
	cs int P=300005;
	int prime[P],pcnt;
	bool mark[P];
	
	inline void init(int len=P-5){
		mark[1]=true;
		for(int re i=2;i<=len;++i){
			if(!mark[i])prime[++pcnt]=i;
			for(int re j=1;j<=pcnt&&i*prime[j]<=len;++j){
				mark[i*prime[j]]=true;
				if(i%prime[j]==0)break;
			}
		}
	}
}

namespace Find_root{
	#define complex COMPLEX
	using namespace Linear_sieves;
	
	inline ll mul(cs ll &a,cs ll &b,cs ll &mod){
		return (a*b-(ll)((long double)a/mod*b)*mod+mod)%mod;
	}
	inline ll quickpow(ll a,ll b,cs ll &mod,ll res=1){
		for(;b;b>>=1,a=mul(a,a,mod))if(b&1)res=mul(res,a,mod);
		return res;
	}
	
	inline ll ex_gcd(cs ll &a,cs ll &b,ll &x,ll &y){
		if(!b){
			y=0;
			x=1;
			return a;
		}
		ll t=ex_gcd(b,a-a/b*b,y,x);
		y-=(a/b)*x;
		return t;
	}
	inline ll inv(cs ll a,cs ll mod){
		ll x,y;
		ll t=ex_gcd(a,mod,x,y);
		return (x%mod+mod)%mod;
	}
	
	ll W,Mod;
	class complex{
		public:
			ll x,y;
			complex(cs ll &_x=0,cs ll &_y=0):x(_x),y(_y){}
			inline friend complex operator*(cs complex &a,cs complex &b){
				return complex(
					(mul(a.x,b.x,Mod)+mul(mul(a.y,b.y,Mod),W,Mod))%Mod,
					(mul(a.x,b.y,Mod)+mul(a.y,b.x,Mod))%Mod);
			}
	};
	
	complex quickpow(complex a,ll b){
		complex res(1,0);
		for(;b;b>>=1,a=a*a)if(b&1)res=res*a;
		return res;
	}
	
	inline bool isprime(ll x){
		if(x<=P-5)return !mark[x];
		if(x%2==0||x%3==0||x%5==0||x%7==0||x%31==0||x%24251==0)return false;
		re ll t=x-1,s;
		t>>=(s=__builtin_ctzll(t));
		for(int re i=1;i<=5;++i){
			re ll p=prime[rand()%pcnt+1];
			re ll num=quickpow(p,t,x),pre=num;
			for(int re j=0;j<s;++j){
				num=mul(num,num,x);
				if(num==1&&pre!=x-1&&pre!=1)return false;
				pre=num;
				if(num==1)break;
			}
			if(num^1)return false;
		}
		return true;
	}
	
	inline ll Pollard_rho(ll x){
		if(x%2==0)return 2;
		if(x%3==0)return 3;
		if(x%5==0)return 5;
		if(x%7==0)return 7;
		if(x%31==0)return 31;
		if(x%24251==0)return 24251;
		re ll n=0,m=0,t=1,q=1,c=rand()%(x-2)+2;
		for(int re k=2;;k<<=1,m=n,q=1){
			for(int re i=1;i<=k;++i){
				n=(mul(n,n,x)+c)%x;
				q=mul(q,abs(n-m),x);
			}
			if((t=__gcd(q,x))>1)return t;
		}
	}
	
	ll fact[60],cntf;
	inline void sieves(ll x){
		if(x==1)return ;
		if(isprime(x)){fact[++cntf]=x;return;}
		re ll p=x;
		while(p==x)p=Pollard_rho(p);
		sieves(p);
		while(x%p==0)x/=p;
		sieves(x);
	}
	
	inline ll solve_2k(ll a,ll k){
		if(a%(1<<k)==0)return 0;
		a%=(1<<k);
		re ll t=0,res=1;
		a>>=(t=__builtin_ctzll(a));
		if((a&7)^1)return -1;
		if(t&1)return -1;
		k-=t;
		for(int re i=4;i<=k;++i){
			res=(res+(a%(1<<i)-res*res)/2)%(1<<k);
		}
		res%=1<<k;
		if(res<0)res+=1<<k;
		return res<<(t>>1); 
	}
	
	inline ll solve_p(ll a,ll p){
		a%=p;
		if(quickpow(a,(p-1)>>1,p)==p-1)return -1;
		re ll b;
		Mod=p;
		while(true){
			b=rand()%p;
			W=(mul(b,b,p)-a+p)%p;
			if(quickpow(W,(p-1)>>1,p)==p-1)break;
		}
		re ll ans=quickpow(complex(b,1),(p+1)>>1).x;
		return min(ans,p-ans);
	}
	
	inline ll solve_pk(ll a,ll k,ll p,ll mod){
		if(a%mod==0)return 0;
		re ll t=0,hmod=1;
		while(a%p==0)a/=p,++t,hmod*=(t&1)?p:1;
		if(t&1)return -1;
		k-=t;
		mod/=hmod*hmod;
		re ll res=solve_p(a,p);
		if(res==-1)return -1;
		complex tmp(res,1);
		W=a;
		Mod=mod;
		tmp=quickpow(tmp,k);
		res=mul(tmp.x,inv(tmp.y,Mod),Mod);
		return res*hmod;
	}
	
	ll remain[20],mod[20],p;
	inline ll CRT(){
		re ll ans=0;
		for(int re i=1;i<=cntf;++i){
			ans=(ans+mul(mul(p/mod[i],inv(p/mod[i],mod[i]),p),remain[i],p))%p;
		}
		return ans;
	}
	
	inline ll solve(ll a,ll pmod){
		a%=pmod;
		cntf=0;
		p=pmod;
		sieves(pmod);
		if(cntf>1)sort(fact+1,fact+cntf+1);
		if(cntf>1)cntf=unique(fact+1,fact+cntf+1)-fact-1;
		for(int re i=1;i<=cntf;++i){
			re ll now=0,rmod=1;
			while(pmod%fact[i]==0)pmod/=fact[i],++now,rmod*=fact[i];
			mod[i]=rmod;
			if(fact[i]==2)remain[i]=solve_2k(a,now);
			else remain[i 
 

            
  
           

              
              
            
            
相關推薦
			   
            
            
            
 

    

    
    
2018.12.30【NOIP訓練】任意模數二次剩餘（高階數論大雜燴）
     
  
 
  
  
 描述 
 求解關於xx的方程： 
      
       
        
         
          
           x
          
          
           2
          
         
         
    

  
 

    

    
    
2018.12.30【NOIP訓練】【SCOI2018】Numazu 的蜜柑（二次剩餘）
     
  
 
  
  
 題面傳送門 
  
 解析： 
 直接解方程可以得到
     
      
       
        
         
          a
         
         
          u
         
        
        
      

  
 

    

    
    
2018.12.31【NOIP訓練】三七二十一（生成函式）
     
  
 
  
  
 傳送門 
  
 解析： 
 設
     
      
       
        
         n
        
       
       
        n
       
      
     n位數的答案為
     
      
       
   

  
 

    

    
    
2018.12.31【NOIP訓練】偶數個5（生成函式）（快速冪）
     
  
 
  
  
 傳送門 
  
 解析： 
 設
     
      
       
        
         
          a
         
         
          n
         
        
       
       
         

  
 

    

    
    
2018.12.30【NOIP提高組】模擬A組 JZOJ 5354 導彈攔截
     
  
 
  
  
 大意 
 三維導彈攔截 
  
 思路 
 對於最多幹掉的導彈數，由於導彈同時發射的，所以我們可以排序，然後
     
      
       
        
         O
        
        
         (
        
        
  

  
 

    

    
    
2018.12.30【NOIP提高組】模擬賽C組總結
     
  
 
  
  
 2018.12.30【NOIP提高組】模擬賽C組總結 
 今天成功迴歸開始做比賽 感覺十分良（zhōng）好（chà）。 
  
   統計數字(count.pas/c/cpp) 
   字串的展開(expand.pas/c/cpp) 
   矩陣取數遊戲(game.pas/c/cpp) 

  
 

    

    
    
2018.11.07【NOIP訓練】lzy的遊戲（01揹包）
     
  
 
  
  
 傳送門 
  
 解析： 
 一個月前口胡了一下這道題，然後現在才在OJ上找到。。。 
 其實最困擾的是後效性的處理，我第一次口胡的時候總是覺得這個後效性環形怎麼處理都不太對，要麼不對要麼複雜度會爆炸，但是其實我們只需要知道一件事情，就是我們選擇的卡牌不超過總的張數，我們就總有辦法構造出 

  
 

    

    
    
2018.11.01【NOIP訓練】某種密碼（折半搜尋）
     
 
								
								            
							
							
							傳送門

解析：
這道題主要是折半的思想，搜尋倒沒有那麼重要。
考慮直接搜尋列舉選或不選2402^{40}240是不現實的，
那直接把數列分成兩半，搜尋兩個2202^{20}220，這個複雜度是可以接受 

  
 

    

    
    
2018.12.30【國家集訓隊】【洛谷P1903】數顏色 / 維護佇列（帶修莫隊）
     
  
 
  
  
 傳送門 
  
 解析： 
 這道題好像以前在BZOJ上做過。 
 但是因為BZOJ資料較水，所以被我複雜度不對的程式碼搞過去了。。 
 真正的排序策略應該是這樣的： 
 塊大小設定成
     
      
       
        
         
          n 

  
 

    

    
    
2018.12.08【NOIP提高組】模擬B組 總結
     
  
 
  
  
 與紀中的一群大佬比賽，竟然水到了
     
      
       
        
         R
        
        
         a
        
        
         k
        
        
         1 

  
 

    

    
    
2018.12.08【NOIP提高組】模擬B組 JZOJ 5223 B
     
  
 
  
  
 描述 
 給定一個
     
      
       
        
         3
        
        
         ×
        
        
         3
        
       
       
        3\ 

  
 

    

    
    
2018.12.08【NOIP提高組】模擬B組 JZOJ 3518 進化序列
     
  
 
  
  
 描述 
 有
     
      
       
        
         
          A
         
         
          1
         
        
        
         ,
        
      

  
 

    

    
    
2018.12.08【NOIP提高組】模擬B組 JZOJ 100042 保留道路
     
  
 
  
  
 描述 
 很久很久以前有一個國家，這個國家有N個城市，城市由1,2,3,…,
     
      
       
        
         N
        
       
       
        N
       
      
     N標號，城市間有M 

  
 

    

    
    
2018.12.08【NOIP提高組】模擬B組 JZOJ 5123 diyiti
     
  
 
  
  
 描述 
 給定
     
      
       
        
         n
        
       
       
        n
       
      
     n根長度為
     
      
       
        
     

  
 

    

    
    
2018.12.15【NOIP提高組】模擬B組 JZOJ 100047 基因變異
     
  
 
  
  
 題目 
 JZOJ 100047 基因變異 
  
 思路 
 設
     
      
       
        
         x
        
        
          
        
        
         x
       

  
 

    

    
    
2018.12.15【NOIP提高組】模擬B組 JZOJ 100046 收集卡片
     
  
 
  
  
 題目 
 JZOJ 100046 收集卡片 
  
 思路 
 維護一個區間，列舉結束訂閱的時間判斷能否晚點訂閱，開一個變數模擬指標維護即可。 ——鳴謝
      
       
        
         
          w
         
         
  

  
 

    

    
    
2018.12.15【NOIP提高組】模擬B組 題目
     
  
 
  
  
 T1 收集卡片 
 一個長度為
     
      
       
        
         n
        
       
       
        n
       
      
     n的字串，每次可以刪除一段連續的子串，求最小刪除子串長度和使得刪 

  
 

    

    
    
2018.12.15【NOIP提高組】模擬B組
     
  
 
  
  
 
 
  解題報告
  
   JZOJ 100046 收集卡片
   
    題目
    分析
    程式碼
   
   JZOJ 100047 基因變異
   
    題目
    分析
    程式碼
   
   JZOJ 100044 abcd
   
    題目
 

  
 

    

    
    
2018.12.08【NOIP提高組】模擬B組
     
  
 
  
  
 
 
  遲到的解題報告
  
   JZOJ 5123 diyiti
   
    分析
    程式碼
   
   JZOJ 100042 保留道路
   
    題目大意
    分析
    程式碼
   
   JZOJ 3518 進化序列
   
    題目大意
   

  
 

    

    
    
2018.12.22【NOIP提高組】模擬B組 總結
     
  
 
  
  
 Part 1(8:20~10:00) 
 剛看道題覺得
     
      
       
        
         T
        
        
         1
        
       
       
        T1