Re: efficient foursquare() and/or threesquare1m4() functions

[Bill Allombert <Bill.Allombert@math.u-bordeaux.fr>]

On Fri, Oct 06, 2023 at 01:06:31AM +0200, hermann@stamm-wilbrandt.de wrote:
> In this posting Bill provided function "foursquare(n)":
> https://pari.math.u-bordeaux.fr/archives/pari-users-2310/msg00004.html
> 
> foursquare(n) = abs(qfsolve(matdiagonal([1,1,1,1,-n]))[1..4]);

I join a script that should handle all cases.
(two/three/four squares with or without factorization, except
for two squares without factoring, which I do not know how to do.) 

Cheers,
Bill

twosquares(n)=abs(qfbsolve(Qfb(1,0,1),n))
threesquares_fact(n)=abs(qfsolve(matdiagonal([1,1,1,-n]))[1..3])~
fouresquares_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(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(isfact(n),return(foursquares_fact(n)));
  forstep(i=sqrtint(n),1,-1,
    my(P=n-i^2, v = valuation(P,2)\2);
    if (P/4^v%8!=7,
      my(F=factor(P,2^20)[,1]);
      if(ispseudoprime(F[#F]),
        return(concat(i,threesquares_fact(P))))));
}

test()=
{
  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));
}