Re: foursquares.gp

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

On Sun, Nov 19, 2023 at 07:47:05PM +0100, Bill Allombert wrote:
> > ? sq=foursquares_fact(2^756839-1);
> >   ***   at top-level: sq=foursquares_fact(2^756839-1)
> >   ***                    ^----------------------------
> >   ***   in function foursquares_fact: abs(
> >   ***   qfsolve(matdiagonal([1,1,1,1,-n]))[1..4])~
> >   ***   ^------------------------------------------
> >   *** qfsolve: the PARI stack overflows !
> >   current stack size: 15000002560 (14305.117 Mbytes)
> >   [hint] you can increase 'parisizemax' using default()
> > 
> >   ***   Break loop: type 'break' to go back to GP prompt
> > break>
> > 
> > 
> > How can it be that computing "sq=foursquares(2^756839-1);" worked (which
> > calls foursquares_fact()), and calling
> > it directly fails?

? setrand(1); sq=foursquares_fact(2^756839-1);
  *** qfsolve: Warning: increasing stack size to 64000000.
 ?  ##
  ***   last result computed in 3h, 34min, 7,871 ms.
? \v
    GP/PARI CALCULATOR Version 2.16.1 (development 28847-9ca7583835)

So indeed if it needs more than 15G, there is certainly a bug.

In any case, I have commited a change 6957a927 which can be used to
speed up a bit foursquares by avoiding factoring the same number twice.
I join a suitable version of foursquares.gp to match.

old:
? foursquares(2^9689-1);
? ##
  ***   last result computed in 675 ms.
new:
? foursquares(2^9689-1);
? ##
  ***   last result computed in 485 ms.

Cheers,
Bill.

/*
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.
*/

/* this version requires PARI 2.16.1 and up */

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

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));
}