Skip to content

Instantly share code, notes, and snippets.

@Hermann-SW
Last active February 29, 2024 02:18
Show Gist options
  • Save Hermann-SW/f13b8adf7d7e3094f0b6db0bce29a7b8 to your computer and use it in GitHub Desktop.
Save Hermann-SW/f13b8adf7d7e3094f0b6db0bce29a7b8 to your computer and use it in GitHub Desktop.
Generate ternary quadratic form that represents n and has determinant 1, for determining sum of 3(4) squares for n
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"));
}
@Hermann-SW
Copy link
Author

Hermann-SW commented Dec 20, 2023

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:~ $ 

@Hermann-SW
Copy link
Author

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:~ $ 

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment