/* Copyright (C) 2023 The PARI group PARI/GP is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version. It is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY WHATSOEVER. Check the License for details. You should have received a copy of it, along with the package; see the file 'COPYING'. If not, write to the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */ twosquares(n)=abs(qfbsolve(Qfb(1,0,1),n,2)); threesquares_fact(n)=abs(qfsolve(matdiagonal([1,1,1,-n]))[1..3])~; isthreesquares(n)=n/4^(valuation(n,2)\2)%8!=7; foursquares_fact(n)=abs(qfsolve(matdiagonal([1,1,1,1,-n]))[1..4])~; isfact(n)=my(F=factor(n,2^20)[,1]);ispseudoprime(F[#F]); threesquares(n)= { if(n<0,error("n not a sum of three squares")); if(n<=1,return([n,0,0])); if(!isthreesquares(n),error("n not a sum of three squares")); if(isfact(n),return(threesquares_fact(n))); forstep(i=sqrtint(n),1,-1, my(P=n-i^2,v = valuation(P,2), Q = P/2^v); if(Q%4==1 && ispseudoprime(Q),return(concat(i,twosquares(P))))); } foursquares(n)= { if(n<0,error("n not a sum of four squares")); if(n<=1,return([n,0,0,0])); if(isfact(n),return(foursquares_fact(n))); forstep(i=sqrtint(n),1,-1, my(P=n-i^2); if (isthreesquares(P) && isfact(P), return(concat(i,threesquares_fact(P))))); } foursquares_test()= { my(P=randomprime(2^500,Mod(7,8))); my(F=foursquares(P)); if (#F!=4 || norml2(F)!=P,error("foursquares",P)); my(P=randomprime(2^500,Mod(1,8))*randomprime(2^500,Mod(7,8))); my(F=foursquares(P)); if (#F!=4 || norml2(F)!=P,error("foursquares",P)); my(P=randomprime(2^500,Mod(1,8))*randomprime(2^500,Mod(3,8))); my(F=threesquares(P)); if (#F!=3 || norml2(F)!=P,error("threesquares",P)); my(P=randomprime(2^500,Mod(1,8))*randomprime(2^500,Mod(5,8))); my(F=threesquares(P)); if (#F!=3 || norml2(F)!=P,error("threesquares",P)); my(P=4*randomprime(2^500,Mod(1,4))); my(F=twosquares(P)); if (#F!=2 || norml2(F)!=P,error("twosquares",P)); }