Generate ternary quadratic form that represents n and has determinant 1, for determining sum of 3(4) squares for n

tqf.gp

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

dbg(a,b,c="")=if(getenv("dbg"),print(a,b,c));

get_tqf(n,vstart)={
    assert(n%4!=0);
    assert(n%8!=7);
    my(Q,v,b,p,a12);
    if(n%4==2,
        for(vv=vstart,oo,
            if(ispseudoprime((4*vv+1)*n-1),
                v=vv; break()));
        write("/dev/stderr",v+1);
        b=4*v+1;
        p=b*n-1;
        assert(kronecker(-b,p)==1);
        a12=lift(sqrt(Mod(-b,p)));
        ,
        my(c=4-(n%4));
        assert(((c*n-1)/2)%2==1);
        for(vv=vstart,oo,
            if(ispseudoprime(((8*vv+c)*n-1)/2),
                v=vv; break()));
        write("/dev/stderr",v+1);
        b=8*v+c;
        p=(b*n-1)/2;
        assert(kronecker(-b,p)==1);
        my(r0=lift(sqrt(Mod(-b,p))));
        a12=abs(r0-(1-(r0%2))*p));
    [(b+a12^2)/(b*n-1),a12,1;a12,b*n-1,0;1,0,n];
}

squares34(n)={
    my(Q,G,V=valuation(n,4),n1=n>>(2*V),res=[],x=[0,0,1]~);
    if(n1%8==7,n=sqrtint(n1);if(n%4==0,n-=1);res=[n*2^V];n1-=n^2;);
    vstart=if(getenv("vstart"),eval(getenv("vstart")),0);

    Q=get_tqf(n1,vstart);
    if(type(Q)!="t_MAT",return([]));
    assert(qfeval(Q,x)==n1);
    assert(matdet(Q)==1);

    dbg(" Q=",Q);
    G=qflllgram(Q,flag=1);
    dbg(" G=",G);
    assert(G~*Q*G==matdiagonal([1,1,1]));

    res=concat([2^V*i|i<-G^-1*x],res);
    abs(res);
}


{
    if(getenv("n")!=0,
        n=eval(getenv("n"));
    
        sq=squares34(n);
        print(getenv("n"),"=norml2(",sq,")");
        assert(norml2(sq)==n);

        print("all asserts OK"));
}

tqf.gp can be used to easily implement Gauss's Eureka theorem:
https://warwick.ac.uk/fac/sci/maths/people/staff/michaud/thirdyearessay.pdf#page=27

$ N=42 gp -q <eureka.gp 2>/dev/null
Δ(2)+Δ(8)+Δ(2)
=3+36+3
=42
$ N=1234567891 gp -q <eureka.gp 2>/dev/null
Δ(39136)+Δ(879)+Δ(30605)
=765832816+386760+468348315
=1234567891
$

$ cat eureka.gp
read("tqf.gp");

delta(n)=n*(n+1)/2;

assert(getenv("N")!=0);
n=eval(getenv("N"));

sq3=squares34(8*n+3)\2;
print("Δ(",sq3[1],")+Δ(",sq3[2],")+Δ(",sq3[3],")");
print("=",delta(sq3[1]),"+",delta(sq3[2]),"+",delta(sq3[3]));
print("=",delta(sq3[1])+delta(sq3[2])+delta(sq3[3]));
$