Last active
February 29, 2024 02:18
-
-
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
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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")); | |
} |
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
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:~ $
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
I learned on pari-users list from Bill how to determine all sums of three squares using qfminim():$r_3(n) = 12*h(-4*n)$ :
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
https://en.wikipedia.org/wiki/Sum_of_squares_function#k_=_3