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

Literature:

get_tqf() sets values as on slide 23:
https://warwick.ac.uk/fac/sci/maths/people/staff/michaud/threesquarestalk.pdf#page=23

Same author article with much more detailed information:
https://warwick.ac.uk/fac/sci/maths/people/staff/michaud/thirdyearessay.pdf

However, finding such a matrix M is not a simple task, we refer to [7, pp. 49–51] for some detailed examples.

Example 1 from [7] E. Grosswald. Representations of Integers as Sums of Squares. SpringerVerlag, New York, 1985:

Found interesting paper, using ternary quadratic forms as well, but different intial matrix M:

Dirichlet's proof of the three-square theorem: an algorithmic perspective
PAUL POLLACK AND PETER SCHORN
Math. Comp. 88 (2019), 1007--1019

Available on Pollack's website:
https://pollack.uga.edu/finding3squares-6.pdf

tqf.gp determines sum of three squares if and only if $n$ is not of the form $n = 4^a(8b + 7)$ for nonnegative integers $a$ and $b$.

pi@raspberrypi400:~ $ n=255 gp -q < tqf.gp
255=[2, 13, 9, 1]
all asserts OK
pi@raspberrypi400:~ $ dbg=1 n=255 gp -q < tqf.gp
 Q=[797, 1622, 1; 1622, 3301, 0; 1, 0, 254]
 G=[-2, 521, -753; 1, -256, 370; 0, -2, 3]
255=[2, 13, 9, 1]
all asserts OK
pi@raspberrypi400:~ $ 

Determines sum of 3[4] squares for n=1..127:

Handling of $n$ of the form $4^a*(8*b+7)$ with nonnegative integers $a$ and $b$ has been tuned.
Mersenne primes $M_p=2^p-1$ are of that form and don't have sum of three squares representation.

Under "Contributed GP scripts" this script can be found, to compute sums of 2, 3 or 4 squares:
https://pari.math.u-bordeaux.fr/Scripts/foursquares.gp
Running with gp 2.15.3 was really bad for "foursquares()".
But gp 2.16.1 on Ubuntu 22.04 is much better.
All measurements below were done on AMD 7600X CPU, with stable gp 2.15.4.

n	p	squares34(M_p) [s]	foursquares(M_p) [s]
22	9941	9.664	0.599
23	11213	15.080	0.798
24	19937	37.610	3.371
25	21701	463.653	3.924
26	23209	95.273	4.591
27	44497	4227.301	22.315
28	86243	47109.029	112.653
29	110503	-	206.783
30	132049	-	320.044
31	216091	-	Killed

$ gp -q tqf.gp
? sq=squares34(2^23209-1);
? ##
  ***   last result computed in 1min, 35,273 ms.
? #sq
4
? norml2(sq)==2^23209-1
1
?

$ gp -q foursquares.gp
? sq=foursquares(2^23209-1);
? ##
  ***   last result computed in 4,642 ms.
? #sq
4
? norml2(sq)==2^23209-1
1
?

Next step is tuning of this loop in "get_tqf(n)", because it takes nearly all the runtime:

        for(vv=0,oo,
            if(ispseudoprime((4*vv+1)*n-1),
                v=vv; break()));

Last update was cosmetic, to make printed output correct.

Before:

pi@raspberrypi400:~ $ n=255 gp -q < tqf.gp
255=[2, 13, 9, 1]
all asserts OK
pi@raspberrypi400:~ $

Now:

pi@raspberrypi5:~ $ n=255 gp -q < tqf.gp 
255=norml2([2, 5, 1, 15])
all asserts OK
pi@raspberrypi5:~ $ 

I learned on pari-users list from Bill how to determine all sums of three squares using qfminim():
https://pari.math.u-bordeaux.fr/archives/pari-users-2312/msg00074.html
In S2.b.gp from that posting qflllgram() is not used.
And number of determined vectors is identical to $r_3(n) = 12*h(-4*n)$:
https://en.wikipedia.org/wiki/Sum_of_squares_function#k_=_3

pi@raspberrypi5:~ $ n=101 gp -q < S2.b.gp 
101=norml2([1, 6, 8])
all asserts OK
528979 [711, -153, -7]~ [2, 9, 4]~
528979 [-711, 153, 7]~ [-2, -9, -4]~
10232019 [3127, -673, -31]~ [2, 9, -4]~
10232019 [-3127, 673, 31]~ [-2, -9, 4]~
13150419 [3545, -763, -35]~ [0, -1, 10]~
13150419 [-3545, 763, 35]~ [0, 1, -10]~
43831329 [6472, -1393, -64]~ [-1, -6, 8]~
43831329 [-6472, 1393, 64]~ [1, 6, -8]~
96137219 [9585, -2063, -95]~ [0, -1, -10]~
96137219 [-9585, 2063, 95]~ [0, 1, 10]~
133712449 [11304, -2433, -112]~ [-1, -6, -8]~
133712449 [-11304, 2433, 112]~ [1, 6, 8]~
427236237 [20206, -4349, -200]~ [-1, -8, 6]~
427236237 [-20206, 4349, 200]~ [1, 8, -6]~
594231237 [23830, -5129, -236]~ [-1, -8, -6]~
594231237 [-23830, 5129, 236]~ [1, 8, 6]~
790321533 [27482, -5915, -272]~ [1, 0, 10]~
790321533 [-27482, 5915, 272]~ [-1, 0, -10]~
1073278614 [32026, -6893, -317]~ [2, 4, 9]~
1073278614 [-32026, 6893, 317]~ [-2, -4, -9]~
1175890933 [33522, -7215, -332]~ [1, 0, -10]~
1175890933 [-33522, 7215, 332]~ [-1, 0, 10]~
1292732233 [35148, -7565, -348]~ [-1, -10, 0]~
1292732233 [-35148, 7565, 348]~ [1, 10, 0]~
1468551054 [37462, -8063, -371]~ [2, 4, -9]~
1468551054 [-37462, 8063, 371]~ [-2, -4, 9]~
4063947019 [62319, -13413, -617]~ [4, 9, 2]~
4063947019 [-62319, 13413, 617]~ [-4, -9, -2]~
4223026299 [63527, -13673, -629]~ [4, 9, -2]~
4223026299 [-63527, 13673, 629]~ [-4, -9, 2]~
4468603938 [65348, -14065, -647]~ [0, -10, 1]~
4468603938 [-65348, 14065, 647]~ [0, 10, -1]~
4551590738 [65952, -14195, -653]~ [0, -10, -1]~
4551590738 [-65952, 14195, 653]~ [0, 10, 1]~
4764374329 [67476, -14523, -668]~ [1, -6, 8]~
4764374329 [-67476, 14523, 668]~ [-1, 6, -8]~
5471166489 [72308, -15563, -716]~ [1, -6, -8]~
5471166489 [-72308, 15563, 716]~ [-1, 6, 8]~
5767594747 [74241, -15979, -735]~ [4, 7, 6]~
5767594747 [-74241, 15979, 735]~ [-4, -7, -6]~
6344416747 [77865, -16759, -771]~ [4, 7, -6]~
6344416747 [-77865, 16759, 771]~ [-4, -7, 6]~
6781754154 [80504, -17327, -797]~ [4, 6, 7]~
6781754154 [-80504, 17327, 797]~ [-4, -6, -7]~
6901225957 [81210, -17479, -804]~ [1, -8, 6]~
6901225957 [-81210, 17479, 804]~ [-1, 8, -6]~
7479857494 [84546, -18197, -837]~ [2, -4, 9]~
7479857494 [-84546, 18197, 837]~ [-2, 4, -9]~
7512803914 [84732, -18237, -839]~ [4, 6, -7]~
7512803914 [-84732, 18237, 839]~ [-4, -6, 7]~
7530904237 [84834, -18259, -840]~ [1, -8, -6]~
7530904237 [-84834, 18259, 840]~ [-1, 8, 6]~
8472634894 [89982, -19367, -891]~ [2, -4, -9]~
8472634894 [-89982, 19367, 891]~ [-2, 4, 9]~
9674396433 [96152, -20695, -952]~ [1, -10, 0]~
9674396433 [-96152, 20695, 952]~ [-1, 10, 0]~
11793127002 [106160, -22849, -1051]~ [4, 2, 9]~
11793127002 [-106160, 22849, 1051]~ [-4, -2, -9]~
13031800602 [111596, -24019, -1105]~ [4, 2, -9]~
13031800602 [-111596, 24019, 1105]~ [-4, -2, 9]~
14788772059 [118881, -25587, -1177]~ [2, -9, 4]~
14788772059 [-118881, 25587, 1177]~ [-2, 9, -4]~
15395980059 [121297, -26107, -1201]~ [2, -9, -4]~
15395980059 [-121297, 26107, 1201]~ [-2, 9, 4]~
17736274062 [130190, -28021, -1289]~ [6, 8, 1]~
17736274062 [-130190, 28021, 1289]~ [-6, -8, -1]~
17901226262 [130794, -28151, -1295]~ [6, 8, -1]~
17901226262 [-130794, 28151, 1295]~ [-6, -8, 1]~
18349082122 [132420, -28501, -1311]~ [4, -2, 9]~
18349082122 [-132420, 28501, 1311]~ [-4, 2, -9]~
19311678947 [135849, -29239, -1345]~ [6, 7, 4]~
19311678947 [-135849, 29239, 1345]~ [-6, -7, -4]~
19886508202 [137856, -29671, -1365]~ [4, -2, -9]~
19886508202 [-137856, 29671, 1365]~ [-4, 2, 9]~
20004682467 [138265, -29759, -1369]~ [6, 7, -4]~
20004682467 [-138265, 29759, 1369]~ [-6, -7, 4]~
25023013094 [154638, -33283, -1531]~ [6, 4, 7]~
25023013094 [-154638, 33283, 1531]~ [-6, -4, -7]~
26410041534 [158866, -34193, -1573]~ [6, 4, -7]~
26410041534 [-158866, 34193, 1573]~ [-6, -4, 7]~
26549203674 [159284, -34283, -1577]~ [4, -6, 7]~
26549203674 [-159284, 34283, 1577]~ [-4, 6, -7]~
27977342554 [163512, -35193, -1619]~ [4, -6, -7]~
27977342554 [-163512, 35193, 1619]~ [-4, 6, 7]~
28887709947 [166151, -35761, -1645]~ [4, -7, 6]~
28887709947 [-166151, 35761, 1645]~ [-4, 7, -6]~
30161621067 [169775, -36541, -1681]~ [4, -7, -6]~
30161621067 [-169775, 36541, 1681]~ [-4, 7, 6]~
31287979089 [172916, -37217, -1712]~ [7, 6, 4]~
31287979089 [-172916, 37217, 1712]~ [-7, -6, -4]~
31692784539 [174031, -37457, -1723]~ [6, 1, 8]~
31692784539 [-174031, 37457, 1723]~ [-6, -1, -8]~
32168405089 [175332, -37737, -1736]~ [7, 6, -4]~
32168405089 [-175332, 37737, 1736]~ [-7, -6, 4]~
33477128219 [178863, -38497, -1771]~ [6, 1, -8]~
33477128219 [-178863, 38497, 1771]~ [-6, -1, 8]~
34088561899 [180489, -38847, -1787]~ [4, -9, 2]~
34088561899 [-180489, 38847, 1787]~ [-4, 9, -2]~
34546393659 [181697, -39107, -1799]~ [4, -9, -2]~
34546393659 [-181697, 39107, 1799]~ [-4, 9, 2]~
35985153829 [185442, -39913, -1836]~ [7, 4, 6]~
35985153829 [-185442, 39913, 1836]~ [-7, -4, -6]~
36655393619 [187161, -40283, -1853]~ [6, -1, 8]~
36655393619 [-187161, 40283, 1853]~ [-6, 1, -8]~
37405376989 [189066, -40693, -1872]~ [7, 4, -6]~
37405376989 [-189066, 40693, 1872]~ [-7, -4, 6]~
38572516179 [191993, -41323, -1901]~ [6, -1, -8]~
38572516179 [-191993, 41323, 1901]~ [-6, 1, 8]~
43686366034 [204324, -43977, -2023]~ [8, 6, 1]~
43686366034 [-204324, 43977, 2023]~ [-8, -6, -1]~
43945029474 [204928, -44107, -2029]~ [8, 6, -1]~
43945029474 [-204928, 44107, 2029]~ [-8, -6, 1]~
44906644134 [207158, -44587, -2051]~ [6, -4, 7]~
44906644134 [-207158, 44587, 2051]~ [-6, 4, -7]~
46758398654 [211386, -45497, -2093]~ [6, -4, -7]~
46758398654 [-211386, 45497, 2093]~ [-6, 4, 7]~
54282305547 [227759, -49021, -2255]~ [6, -7, 4]~
54282305547 [-227759, 49021, 2255]~ [-6, 7, -4]~
55440035147 [230175, -49541, -2279]~ [6, -7, -4]~
55440035147 [-230175, 49541, 2279]~ [-6, 7, 4]~
57901872782 [235230, -50629, -2329]~ [6, -8, 1]~
57901872782 [-235230, 50629, 2329]~ [-6, 8, -1]~
58103395299 [235639, -50717, -2333]~ [8, 1, 6]~
58103395299 [-235639, 50717, 2333]~ [-8, -1, -6]~
58199603862 [235834, -50759, -2335]~ [6, -8, -1]~
58199603862 [-235834, 50759, 2335]~ [-6, 8, 1]~
59254645269 [237962, -51217, -2356]~ [7, -4, 6]~
59254645269 [-237962, 51217, 2356]~ [-7, 4, -6]~
59904336339 [239263, -51497, -2369]~ [8, 1, -6]~
59904336339 [-239263, 51497, 2369]~ [-8, -1, 6]~
61073205069 [241586, -51997, -2392]~ [7, -4, -6]~
61073205069 [-241586, 51997, 2392]~ [-7, 4, 6]~
64179725829 [247654, -53303, -2452]~ [9, 4, 2]~
64179725829 [-247654, 53303, 2452]~ [-9, -4, -2]~
64758934579 [248769, -53543, -2463]~ [8, -1, 6]~
64758934579 [-248769, 53543, 2463]~ [-8, 1, -6]~
64807361309 [248862, -53563, -2464]~ [9, 4, -2]~
64807361309 [-248862, 53563, 2464]~ [-9, -4, 2]~
66291800409 [251696, -54173, -2492]~ [7, -6, 4]~
66291800409 [-251696, 54173, 2492]~ [-7, 6, -4]~
66659459779 [252393, -54323, -2499]~ [8, -1, -6]~
66659459779 [-252393, 54323, 2499]~ [-8, 1, 6]~
67570563049 [254112, -54693, -2516]~ [7, -6, -4]~
67570563049 [-254112, 54693, 2516]~ [-7, 6, 4]~
70836156177 [260180, -55999, -2576]~ [9, 2, 4]~
70836156177 [-260180, 55999, 2576]~ [-9, -2, -4]~
72157816577 [262596, -56519, -2600]~ [9, 2, -4]~
72157816577 [-262596, 56519, 2600]~ [-9, -2, 4]~
83868562114 [283104, -60933, -2803]~ [8, -6, 1]~
83868562114 [-283104, 60933, 2803]~ [-8, 6, -1]~
84226809714 [283708, -61063, -2809]~ [8, -6, -1]~
84226809714 [-283708, 61063, 2809]~ [-8, 6, 1]~
85856762297 [286440, -61651, -2836]~ [9, -2, 4]~
85856762297 [-286440, 61651, 2836]~ [-9, 2, -4]~
87311201577 [288856, -62171, -2860]~ [9, -2, -4]~
87311201577 [-288856, 62171, 2860]~ [-9, 2, 4]~
93210511219 [298455, -64237, -2955]~ [10, 1, 0]~
93210511219 [-298455, 64237, 2955]~ [-10, -1, 0]~
94287327509 [300174, -64607, -2972]~ [9, -4, 2]~
94287327509 [-300174, 64607, 2972]~ [-9, 4, -2]~
95047741869 [301382, -64867, -2984]~ [9, -4, -2]~
95047741869 [-301382, 64867, 2984]~ [-9, 4, 2]~
97163554038 [304718, -65585, -3017]~ [10, 0, 1]~
97163554038 [-304718, 65585, 3017]~ [-10, 0, -1]~
97549123438 [305322, -65715, -3023]~ [10, 0, -1]~
97549123438 [-305322, 65715, 3023]~ [-10, 0, 1]~
101592175419 [311585, -67063, -3085]~ [10, -1, 0]~
101592175419 [-311585, 67063, 3085]~ [-10, 1, 0]~
#S2=168
12*h(-4*n)=168
101592175419 [311585, -67063, -3085]~ [10, -1, 0]~
pi@raspberrypi5:~ $ 

With this little diff using "qflllgram()" ...

pi@raspberrypi5:~ $ diff S2.b.gp S2.c.gp 
3c3
< Q=get_tqf(n);
---
> Q=qflllgram(get_tqf(n))^-1;
pi@raspberrypi5:~ $ 

... the vector $[0,0,1]$~ from Dirichlet construction ("get_tqf(n)") is first sorted vector:

pi@raspberrypi5:~ $ n=101 gp -q < S2.c.gp 
101=norml2([1, 6, 8])
all asserts OK
1 [0, 0, 1]~ [-1, -6, -8]~
1 [0, 0, -1]~ [1, 6, 8]~
651 [25, -5, -1]~ [6, 1, 8]~
651 [-25, 5, 1]~ [-6, -1, -8]~
6886 [81, -18, -1]~ [-8, 6, -1]~
6886 [-81, 18, 1]~ [8, -6, 1]~
23826 [151, -32, -1]~ [6, -8, -1]~
23826 [-151, 32, 1]~ [-6, 8, 1]~
23891 [151, -33, -1]~ [-8, 1, -6]~
23891 [-151, 33, 1]~ [8, -1, 6]~
40181 [196, -42, -1]~ [1, -8, -6]~
40181 [-196, 42, 1]~ [-1, 8, 6]~
43561 [204, -44, -3]~ [-1, 6, 8]~
43561 [-204, 44, 3]~ [1, -6, -8]~
43970 [205, -44, -3]~ [2, 4, 9]~
43970 [-205, 44, 3]~ [-2, -4, -9]~
48350 [215, -46, -3]~ [4, 2, 9]~
48350 [-215, 46, 3]~ [-4, -2, -9]~
49307 [217, -47, -3]~ [-4, 7, 6]~
49307 [-217, 47, 3]~ [4, -7, -6]~
57835 [235, -51, -3]~ [-6, 7, 4]~
57835 [-235, 51, 3]~ [6, -7, -4]~
59731 [239, -51, -3]~ [6, -1, 8]~
59731 [-239, 51, 3]~ [-6, 1, -8]~
74662 [267, -58, -3]~ [-8, 6, 1]~
74662 [-267, 58, 3]~ [8, -6, -1]~
77357 [272, -58, -3]~ [7, -4, 6]~
77357 [-272, 58, 3]~ [-7, 4, -6]~
94105 [300, -64, -3]~ [7, -6, 4]~
94105 [-300, 64, 3]~ [-7, 6, -4]~
96781 [304, -66, -3]~ [-9, 4, -2]~
96781 [-304, 66, 3]~ [9, -4, 2]~
115417 [332, -72, -3]~ [-9, 2, -4]~
115417 [-332, 72, 3]~ [9, -2, 4]~
118762 [337, -72, -3]~ [6, -8, 1]~
118762 [-337, 72, 3]~ [-6, 8, -1]~
139475 [365, -79, -3]~ [-8, -1, -6]~
139475 [-365, 79, 3]~ [8, 1, 6]~
142411 [369, -79, -3]~ [4, -9, -2]~
142411 [-369, 79, 3]~ [-4, 9, 2]~
156667 [387, -83, -3]~ [2, -9, -4]~
156667 [-387, 83, 3]~ [-2, 9, 4]~
158386 [389, -84, -3]~ [-6, -4, -7]~
158386 [-389, 84, 3]~ [6, 4, 7]~
166606 [399, -86, -3]~ [-4, -6, -7]~
166606 [-399, 86, 3]~ [4, 6, 7]~
167405 [400, -86, -3]~ [-1, -8, -6]~
167405 [-400, 86, 3]~ [1, 8, 6]~
182014 [417, -90, -5]~ [-4, 6, 7]~
182014 [-417, 90, 5]~ [4, -6, -7]~
206377 [444, -96, -5]~ [-7, 6, 4]~
206377 [-444, 96, 5]~ [7, -6, -4]~
228114 [467, -100, -5]~ [6, -4, 7]~
228114 [-467, 100, 5]~ [-6, 4, -7]~
270987 [509, -109, -5]~ [6, -7, 4]~
270987 [-509, 109, 5]~ [-6, 7, -4]~
356957 [584, -126, -5]~ [-7, -4, -6]~
356957 [-584, 126, 5]~ [7, 4, 6]~
375467 [599, -129, -5]~ [-4, -7, -6]~
375467 [-599, 129, 5]~ [4, 7, 6]~
393242 [613, -132, -7]~ [-2, 4, 9]~
393242 [-613, 132, 7]~ [2, -4, -9]~
432542 [643, -138, -7]~ [4, -2, 9]~
432542 [-643, 138, 7]~ [-4, 2, -9]~
478341 [676, -146, -7]~ [-9, 4, 2]~
478341 [-676, 146, 7]~ [9, -4, -2]~
574411 [741, -159, -7]~ [4, -9, 2]~
574411 [-741, 159, 7]~ [-4, 9, -2]~
604545 [760, -164, -7]~ [-9, -2, -4]~
604545 [-760, 164, 7]~ [9, 2, 4]~
661315 [795, -171, -7]~ [-2, -9, -4]~
661315 [-795, 171, 7]~ [2, 9, 4]~
708739 [823, -177, -9]~ [0, 1, 10]~
708739 [-823, 177, 9]~ [0, -1, -10]~
717349 [828, -178, -9]~ [1, 0, 10]~
717349 [-828, 178, 9]~ [-1, 0, -10]~
729706 [835, -180, -9]~ [-6, 4, 7]~
729706 [-835, 180, 9]~ [6, -4, -7]~
745541 [844, -182, -9]~ [-7, 4, 6]~
745541 [-844, 182, 9]~ [7, -4, -6]~
819406 [885, -190, -9]~ [4, -6, 7]~
819406 [-885, 190, 9]~ [-4, 6, -7]~
845531 [899, -193, -9]~ [4, -7, 6]~
845531 [-899, 193, 9]~ [-4, 7, -6]~
872459 [913, -197, -9]~ [-10, 1, 0]~
872459 [-913, 197, 9]~ [10, -1, 0]~
899410 [927, -200, -9]~ [-10, 0, -1]~
899410 [-927, 200, 9]~ [10, 0, 1]~
980369 [968, -208, -9]~ [1, -10, 0]~
980369 [-968, 208, 9]~ [-1, 10, 0]~
998710 [977, -210, -9]~ [0, -10, -1]~
998710 [-977, 210, 9]~ [0, 10, 1]~
1013281 [984, -212, -9]~ [-7, -6, -4]~
1013281 [-984, 212, 9]~ [7, 6, 4]~
1023571 [989, -213, -9]~ [-6, -7, -4]~
1023571 [-989, 213, 9]~ [6, 7, 4]~
1112366 [1031, -222, -11]~ [-4, 2, 9]~
1112366 [-1031, 222, 11]~ [4, -2, -9]~
1114429 [1032, -222, -11]~ [-1, 0, 10]~
1114429 [-1032, 222, 11]~ [1, 0, -10]~
1125219 [1037, -223, -11]~ [0, -1, 10]~
1125219 [-1037, 223, 11]~ [0, 1, -10]~
1177826 [1061, -228, -11]~ [2, -4, 9]~
1177826 [-1061, 228, 11]~ [-2, 4, -9]~
1211721 [1076, -232, -11]~ [-9, 2, 4]~
1211721 [-1076, 232, 11]~ [9, -2, -4]~
1296490 [1113, -240, -11]~ [-10, 0, 1]~
1296490 [-1113, 240, 11]~ [10, 0, -1]~
1329299 [1127, -243, -11]~ [-10, -1, 0]~
1329299 [-1127, 243, 11]~ [10, 1, 0]~
1338331 [1131, -243, -11]~ [2, -9, 4]~
1338331 [-1131, 243, 11]~ [-2, 9, -4]~
1408221 [1160, -250, -11]~ [-9, -4, -2]~
1408221 [-1160, 250, 11]~ [9, 4, 2]~
1415190 [1163, -250, -11]~ [0, -10, 1]~
1415190 [-1163, 250, 11]~ [0, 10, -1]~
1437209 [1172, -252, -11]~ [-1, -10, 0]~
1437209 [-1172, 252, 11]~ [1, 10, 0]~
1469371 [1185, -255, -11]~ [-4, -9, -2]~
1469371 [-1185, 255, 11]~ [4, 9, 2]~
1632531 [1249, -269, -13]~ [-6, 1, 8]~
1632531 [-1249, 269, 13]~ [6, -1, -8]~
1679987 [1267, -273, -13]~ [-8, 1, 6]~
1679987 [-1267, 273, 13]~ [8, -1, -6]~
1725001 [1284, -276, -13]~ [1, -6, 8]~
1725001 [-1284, 276, 13]~ [-1, 6, -8]~
1801037 [1312, -282, -13]~ [1, -8, 6]~
1801037 [-1312, 282, 13]~ [-1, 8, -6]~
1949830 [1365, -294, -13]~ [-8, -6, -1]~
1949830 [-1365, 294, 13]~ [8, 6, 1]~
1978410 [1375, -296, -13]~ [-6, -8, -1]~
1978410 [-1375, 296, 13]~ [6, 8, 1]~
2227502 [1459, -314, -15]~ [-4, -2, 9]~
2227502 [-1459, 314, 15]~ [4, 2, -9]~
2239819 [1463, -315, -15]~ [-6, -1, 8]~
2239819 [-1463, 315, 15]~ [6, 1, -8]~
2258042 [1469, -316, -15]~ [-2, -4, 9]~
2258042 [-1469, 316, 15]~ [2, 4, -9]~
2295347 [1481, -319, -15]~ [-8, -1, 6]~
2295347 [-1481, 319, 15]~ [8, 1, -6]~
2316769 [1488, -320, -15]~ [-1, -6, 8]~
2316769 [-1488, 320, 15]~ [1, 6, -8]~
2367217 [1504, -324, -15]~ [-9, -2, 4]~
2367217 [-1504, 324, 15]~ [9, 2, -4]~
2404757 [1516, -326, -15]~ [-1, -8, 6]~
2404757 [-1516, 326, 15]~ [1, 8, -6]~
2456149 [1532, -330, -15]~ [-9, -4, 2]~
2456149 [-1532, 330, 15]~ [9, 4, -2]~
2478307 [1539, -331, -15]~ [-2, -9, 4]~
2478307 [-1539, 331, 15]~ [2, 9, -4]~
2517382 [1551, -334, -15]~ [-8, -6, 1]~
2517382 [-1551, 334, 15]~ [8, 6, -1]~
2536699 [1557, -335, -15]~ [-4, -9, 2]~
2536699 [-1557, 335, 15]~ [4, 9, -2]~
2549842 [1561, -336, -15]~ [-6, -8, 1]~
2549842 [-1561, 336, 15]~ [6, 8, -1]~
2992266 [1691, -364, -17]~ [-6, -4, 7]~
2992266 [-1691, 364, 17]~ [6, 4, -7]~
3024245 [1700, -366, -17]~ [-7, -4, 6]~
3024245 [-1700, 366, 17]~ [7, 4, -6]~
3027646 [1701, -366, -17]~ [-4, -6, 7]~
3027646 [-1701, 366, 17]~ [4, 6, -7]~
3077675 [1715, -369, -17]~ [-4, -7, 6]~
3077675 [-1715, 369, 17]~ [4, 7, -6]~
3124657 [1728, -372, -17]~ [-7, -6, 4]~
3124657 [-1728, 372, 17]~ [7, 6, -4]~
3142707 [1733, -373, -17]~ [-6, -7, 4]~
3142707 [-1733, 373, 17]~ [6, 7, -4]~
#S2=168
12*h(-4*n)=168
708739 [823, -177, -9]~ [0, 1, 10]~
pi@raspberrypi5:~ $ 

Before the get_tqf() loops always started at 0, now they start at vstart environment variable.
This allows to get more than 1 solution for a given n.
Next vstart start value for (different) output is printed to /dev/stderr:

pi@raspberrypi5:~ $ n=101 gp -q < tqf.gp
1
101=norml2([1, 6, 8])
all asserts OK
pi@raspberrypi5:~ $ vstart=1 n=101 gp -q < tqf.gp
13
101=norml2([6, 4, 7])
all asserts OK
pi@raspberrypi5:~ $ vstart=13 n=101 gp -q < tqf.gp
15
101=norml2([9, 4, 2])
all asserts OK
pi@raspberrypi5:~ $ vstart=15 n=101 gp -q < tqf.gp
16
101=norml2([0, 10, 1])
all asserts OK
pi@raspberrypi5:~ $ 

Redirecting 3 to standard output allows to always use the same command for iteration over possible sum of squares:

pi@raspberrypi5:~ $ exec 3>&1
pi@raspberrypi5:~ $ vstart=
pi@raspberrypi5:~ $ vstart=$(vstart=$vstart n=101 gp -q < tqf.gp 2>&1 1>&3)
101=norml2([1, 6, 8])
all asserts OK
pi@raspberrypi5:~ $ vstart=$(vstart=$vstart n=101 gp -q < tqf.gp 2>&1 1>&3)
101=norml2([6, 4, 7])
all asserts OK
pi@raspberrypi5:~ $ vstart=$(vstart=$vstart n=101 gp -q < tqf.gp 2>&1 1>&3)
101=norml2([9, 4, 2])
all asserts OK
pi@raspberrypi5:~ $ vstart=$(vstart=$vstart n=101 gp -q < tqf.gp 2>&1 1>&3)
101=norml2([0, 10, 1])
all asserts OK
pi@raspberrypi5:~ $