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

首先是一个结论，若 \(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;
}