首页 > 其他分享 >把一个整数拆分成两个平方的和 $x^2+y^2=n$

把一个整数拆分成两个平方的和 $x^2+y^2=n$

时间:2023-06-11 10:11:42浏览次数:36  
标签:分成 node 平方 return int res ll 整数 ans

首先是一个结论,若 \(n=2^\alpha\prod p_i^{\beta_i}\prod q_i^{2\gamma_i}\) ,其中 \(p_i\equiv1 \pmod 4,q_i\equiv 3 \pmod 4\) ,那么有解,否则无解。

然后考虑如果 \(n_1=x_1^2+y_1^2,n_2=x_2^2+y_2^2\) ,那么 \(n_1n_2\) 也是可以被表示的,即 \(n_1n_2=(x_1y_1+x_2y_2)^2+(x_1y_2-x_2y_1)^2\) 。

注意到 \(2^\alpha\) 是好分解的,按照 \(\alpha\) 的奇偶性讨论即可,对于 \(\beta_i,2\gamma_i\) 其实可以直接把偶数次方拆出来最后再乘,那么现在就是要分解 \(p\equiv 1\pmod 4\) 的质数了。

考虑如果构造了 \(x^2+y^2=mp\) ,那么根据这个构造一组 \(a^2+b^2=rp\) ,其中 \(r<m\) ,不断重复就可以得到一组解。

考虑构造一组 \((u,v)\) 满足 \(u\equiv x\pmod m,v\equiv y\pmod m\) ,那么有 \(x^2+y^2\equiv u^2+v^2\pmod m\) ,不妨直接令 \(u^2+v^2=mr\) 。

把两个平方和的式子相加可以得到 \((xu+yv)^2+(xv-yu)^2=m^2rp\) ,那么根据同余的性质,可以得到 \((\frac{xu+yv}{m})^2+(\frac{xy-yu}{m})^2=rp\) ,取 \(\mid u,v\mid\leq \frac{m}{2}\) ,可以使得 \(r\leq \frac{m}{2}\) ,这样就能 \(O(\log V)\) 的构造了。

注意到一开始要构造一组 \(x^2+y^2=mp\) ,可以想到令 \(x^2\equiv -1 \pmod p,y=1\) ,那么直接随机 \(a\) ,判断 \(a^{\frac{p-1}{2}}\) 是否为 \(-1\) ,根据二次剩余的结论,这个随机几乎是 \(O(1)\) 的,得到 \(a\) 之后只需要令 \(x\equiv a^{\frac{p-1}{4}}\pmod p\) 即可。

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef __int128 I;
random_device RD;
mt19937 rnd(RD());
int n;
inline ll add(ll a,ll b,ll p){a+=b;return a>=p?a-p:a;}
inline ll ksc(ll a,ll b,ll p){ll x=(long double)a/p*b;return ((I)a*b-(I)x*p+p)%p;}
inline ll ksm(ll a,ll b,ll p){
    ll res=1;
    a%=p;
    while (b){
        if (b&1) res=ksc(res,a,p);
        a=ksc(a,a,p);
        b>>=1;
    }
    return res;
}
namespace MillerRabin{
	const int P=15;
	const int pr[P]={2,3,5,7,11,13,17,19,23,29,31,37,41,43,19260817};
	inline bool check(ll x,ll p){
		if (x%p==0 || ksm(p%x,x-1,x)!=1) return 0;
		ll a=x-1;
		while (!(a&1)){
			ll t=ksm(p%x,a>>=1,x);
			if (t!=1 && t!=x-1) return 0;
			if (t==x-1) return 1;
		}
		return 1;
	}
	inline bool Isprime(ll x){
		if (x>2 && x%2==0) return 0;
		if (x<=2) return x==2;
		for (int i=0;i<P;++i){
			if (x==pr[i]) return 1;
			if (!check(x,pr[i])) return 0;
		}
		return 1;
	}
}
namespace PollardRho{
    vector<ll> fac;
    ll gcd(ll x,ll y){return !y?x:gcd(y,x%y);}
    void clear(){fac.clear();}
    void unique(){
        sort(fac.begin(),fac.end());
        auto it=unique(fac.begin(),fac.end());
        fac.erase(it,fac.end());
    }
    ll c;
    inline ll f(ll x,ll p){return ((I)x*x+c)%p;}
    ll pollard_rho(ll x){
    	if (x==4) return 2;
    	while (1){
	    	c=uniform_int_distribution<ll>(1,x-1)(rnd);
	    	ll s=f(f(0,x),x),t=f(0,x);
	    	while (s!=t){
	    		ll val=1;
	    		for (int i=0;i<128 && s!=t;++i){
	    			ll tval=(I)val*abs(s-t)%x;
	    			if (!tval) break;
	    			val=tval;
	    			s=f(f(s,x),x),t=f(t,x);
				}
				if (val){
					ll d=gcd(val,x);
					if (d>1) return d;
				}
			}
		}
	}
    void factor(ll x){
        if (x<2) return;
        if (MillerRabin::Isprime(x)){
            fac.emplace_back(x);
            return;
        }
        ll p=pollard_rho(x);
        while ((x%p)==0) x/=p;
        factor(x);factor(p);
    }
    vector<pair<ll,int>> get_factor(ll n){
        clear(),factor(n),unique();
        vector<pair<ll,int>> res(0);
        for (ll p:fac){
            int c=0;
            while (n%p==0) ++c,n/=p;
            res.emplace_back(p,c);
        }
        return res;
    }
}
struct node{ll x,y;node(ll a=0,ll b=0){x=a,y=b;}};
inline node operator + (node a,node b){return node(a.x*b.x+a.y*b.y,abs(a.x*b.y-a.y*b.x));}
ll pow(ll x,ll y){
    ll res=1;
    while (y){
        if (y&1) res*=x;
        x*=x;
        y>>=1;
    }
    return res;
}
ll find(ll p){// 4k+1
    ll x=uniform_int_distribution<ll>(1,p-1)(rnd);
    while (ksm(x,(p-1)/2,p)==1) x=uniform_int_distribution<ll>(1,p-1)(rnd);
    return ksm(x,(p-1)/4,p);
}
I myabs(I x){return x<0?-x:x;}
void work(){
    ll n;scanf("%lld",&n);
    vector<pair<ll,int>> pr=PollardRho::get_factor(n);
    node ans(1,0);
    for (auto x:pr){
        ll p=x.first;int c=x.second;
        if (p==2){
            if (c&1) ans=ans+node(pow(2,c/2),pow(2,c/2));
            else ans=ans+node(pow(2,c/2),0);
        }
        else if (p%4==3){
            if (c&1) return puts("-1"),void();
            ans.x*=pow(p,c/2),ans.y*=pow(p,c/2);
        }
        else{
            ans.x*=pow(p,c/2),ans.y*=pow(p,c/2);
            if (c&1){
                I x=find(p),y=1;
				while (x*x+y*y!=p){
				    I k=(x*x+y*y)/p;
				    I u=x%k,v=y%k;
				    u-=(2*u>k)*k,v-=(2*v>k)*k;
				    tie(x,y)=make_pair((u*y-v*x)/k,(u*x+v*y)/k);
				    x=myabs(x),y=myabs(y);
				}
                ans=ans+node(x,y);
            }
        }
    }
    printf("%lld %lld\n",ans.x,ans.y);
}
int main(){
    int T;scanf("%d",&T);
    while (T--) work();
    return 0;
}

标签:分成,node,平方,return,int,res,ll,整数,ans
From: https://www.cnblogs.com/jerryjiang/p/17472551.html

相关文章

  • Python快速判断若干整数是否互不相同
    让我们先来看一个简单的问题:给定两个整数x和y,如果这两个数不相等就输出Yes,否则输出No。遇到这样的问题,一般都会毫不犹豫地给出类似于下面的代码:如果问题性质不变,简单地增加一下问题规模:给定三个整数x、y和z,如果这三个数互不相等就输出Yes,否则输出No。估计很多人会在上面代码的基础......
  • Python把列表中的数字尽量等分成n份
    问题描述:假设一个列表中含有若干整数,现在要求将其分成n个子列表,并使得各个子列表中的整数之和尽可能接近。下面的代码并没有使用算法,而是直接将原始列表分成n个子列表,然后再不断地调整各个子列表中的数字,从元素之和最大的子列表中拿出最小的元素放到元素之核最小的子列表中,重复这个......
  • Python组合列表中多个整数得到最小整数(一个算法的巧妙实现)
    '''程序功能:  给定一个含有多个整数的列表,将这些整数任意组合和连接,  返回能得到的最小值。  代码思路:  将这些整数变为相同长度(按最大的进行统一),短的右侧使用个位数补齐  然后将这些新的数字升序排列,将低位补齐的数字删掉,  把剩下的数字连接起来,即可得到满足要......
  • Python+tensorflow计算整数阶乘的方法与局限性
    本文代码主要演示tensorflow的基本用法。importtensorflowas#创建变量,保存计算结果start=tf.Variable(1,dtype=tf.int64)#初始化变量的opinit_op=tf.global_variables_initializer()#启用默认图withtf.Session()assess:#初始化变量sess.run(ini......
  • 北京划定63处禁止开发区域 总面积逾3千平方公里 (zz)
    本报北京10月14日电备受关注的《北京市主体功能区规划》近日浮出水面。其中除了人们熟悉的首都功能核心区、城市功能拓展区、城市发展新区、生态涵养发展区四类功能区域外,首次设立“禁止开发区域”,并明确在该区域内,除必要的交通、保护、修复、监测及科学实验设施外,禁止......
  • LeetCode> 69. 求x的平方根
    目录题目题目描述解题思路参考题目地址:LeetCode69.x的平方根题目描述给你一个非负整数x,计算并返回x的算术平方根。由于返回类型是整数,结果只保留整数部分,小数部分将被舍去。注意:不允许使用任何内置指数函数和算符,例如pow(x,0.5)或者x**0.5。示例1:输入:x......
  • 输入正整数N,检查它是否可以被其数字之和整除
    题目:*输入正整数N,检查它是否可以被其数字之和整除,*输出YES或者NO。不考虑不合理的输入等特殊情况。eg:*例如:78的各位数字之和是:7+8=15,则78是一个各位数字之和能被15整除的整数。classTest53{publicstaticvoidmain(String[]args){Scannerinpu......
  • P4451 [国家集训队]整数的lqp拆分
    Description求\[\begin{aligned}&\sum\prod_{i=1}^mF_{a_i}\\&m>0\\&a_1,a_2\ldotsa_m>0\\&a_1+a_2+\ldots+a_m=n\end{aligned}\]由于答案可能非常大,所以要对\(10^9+7\)取模。Solution题目中有整数拆分的部分,可以想到用生成函数的相关知识。设斐波那契数......
  • 一个整数数组里面,除了两个数之外,其他的数字都出现了两次,写一个程序找出这两个数...
    一个整数数组里面,除了两个数之外,其他的数字都出现了两次,写一个程序找出这两个数,要求算法的时间复杂度为O(n). n为数组的长度。  程序代码如下: //取二进制中首个为1的位置intfindFirstOne(intvalue){intpos=0;while((value&1)!=1){value=value>>1;pos++;......
  • LuoguP4318 完全平方数
    标签:莫比乌斯函数,容斥完全平方数题目描述小X自幼就很喜欢数。但奇怪的是,他十分讨厌完全平方数。他觉得这些数看起来很令人难受。由此,他也讨厌所有是完全平方数的正整数倍的数。然而这丝毫不影响他对其他数的热爱。这天是小X的生日,小W想送一个数给他作为生日礼物。当然他......