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

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