Representatins of n as sum of k squares (non-negative monotonic increasing, all, count)

SquaresRepresentations.gp

assert(b)=if(!(b),error());

SquaresRepresentations(n,k,a=0)={
  my(R=List(),m=sqrtint(n\k));
  if(k==1,return(if(m>=a&&n==m^2,[m],[])));
  for(b=a,m,foreach(SquaresRepresentations(n-b^2,k-1,b),s,listput(R,concat([b],s))));
  Vec(R);
}

signs(a)={
  if(#a==0,
    [[]],
    [concat(s,v)|s<-Set([-a[1],a[1]]);v<-signs(a[2..#a])]
  );
}

r3(n)=if(n%8==7,0,if(n%8==3,24*qfbclassno(-n),12*qfbclassno(-4*n)));

n=eval(getenv("n"));
k=eval(getenv("k"));
S=SquaresRepresentations(n,k);
{
  if(getenv("all"),
    L=List();
    foreach(S,s,
      forperm(s,p,
        foreach(signs(Vec(p)),v,listput(L,v))
      );
    );
    if(getenv("all")=="count",
      print("#(signed(permuted(elements)))=",#L);
      if(k==3,
        if(issquarefree(n),
          assert(#L==r3(n));
          print("is equal to r3(n)"),
          print("n not squarefree, so simple r3(n) cannot be used to compare")
        )
      ),
      print(Vec(L));
    ),
    print(S);
  );
}

pi@raspberrypi5:~/gp $ n=17 k=3 gp -q < SquaresRepresentations.gp 
[[0, 1, 4], [2, 2, 3]]
pi@raspberrypi5:~/gp $ n=17 k=3 all=1 gp -q < SquaresRepresentations.gp 
[[0, -1, -4], [0, -1, 4], [0, 1, -4], [0, 1, 4], [0, -4, -1], [0, -4, 1], [0, 4, -1], [0, 4, 1], [-1, 0, -4], [-1, 0, 4], [1, 0, -4], [1, 0, 4], [-1, -4, 0], [-1, 4, 0], [1, -4, 0], [1, 4, 0], [-4, 0, -1], [-4, 0, 1], [4, 0, -1], [4, 0, 1], [-4, -1, 0], [-4, 1, 0], [4, -1, 0], [4, 1, 0], [-2, -2, -3], [-2, -2, 3], [-2, 2, -3], [-2, 2, 3], [2, -2, -3], [2, -2, 3], [2, 2, -3], [2, 2, 3], [-2, -3, -2], [-2, -3, 2], [-2, 3, -2], [-2, 3, 2], [2, -3, -2], [2, -3, 2], [2, 3, -2], [2, 3, 2], [-3, -2, -2], [-3, -2, 2], [-3, 2, -2], [-3, 2, 2], [3, -2, -2], [3, -2, 2], [3, 2, -2], [3, 2, 2]]
pi@raspberrypi5:~/gp $ n=17 k=3 all=count gp -q < SquaresRepresentations.gp 
#(signed(permuted(elements)))=48
is equal to r3(n)
pi@raspberrypi5:~/gp $  

pi@raspberrypi5:~/gp $ n=10009 k=3 all=count gp -q < SquaresRepresentations.gp 
#(signed(permuted(elements)))=1152
is equal to r3(n)
pi@raspberrypi5:~/gp $ 
pi@raspberrypi5:~/gp $ n=10009 k=3 gp -q < SquaresRepresentations.gp 
[[0, 3, 100], [3, 28, 96], [3, 60, 80], [6, 57, 82], [8, 12, 99], [8, 27, 96], [8, 36, 93], [8, 69, 72], [9, 18, 98], [9, 62, 78], [12, 53, 84], [18, 42, 89], [18, 46, 87], [18, 66, 73], [24, 28, 93], [26, 42, 87], [26, 57, 78], [27, 64, 72], [28, 60, 75], [30, 55, 78], [35, 60, 72], [39, 42, 82], [42, 54, 73], [53, 60, 60], [54, 57, 62]]
pi@raspberrypi5:~/gp $ 

Revision 4 did not work for n=100009, k=3 as well as for n=16, k-3 because simple r3(n)
https://en.wikipedia.org/wiki/Sum_of_squares_function#k_=_3
applies to squarefree n only. Revision 5 takes care of that:

hermann@j4105:~$ n=10009 k=3 all=count gp -q < SquaresRepresentations.gp
#(signed(permuted(elements)))=1152
is equal to r3(n)
hermann@j4105:~$ n=100009 k=3 all=count gp -q < SquaresRepresentations.gp
#(signed(permuted(elements)))=3456
n not squarefree, so simple r3(n) cannot be used to compare
hermann@j4105:~$ 

Revision 6:
vecmax(factorint(n)[,2])==1 -> issquarefree(n)