/* Lagrange's Four square theorem (+ optimality) in Picat. Inspired by CryptoMoonWatc1 https://x.com/CryptoMoonWatc1/status/1833558857718661354 """ An easy one: every natural number N is the sum of four natural numbers squared. Given N find the four naturals with minimum supremum minus infimum. """ This is Lagrange's Four square theorem (https://en.wikipedia.org/wiki/Lagrange%27s_four-square_theorem ) + requirement of optimality. Here are some random instances: [n = 524287,x = [346,349,373,379],d = 33] [n = 786431,x = [431,433,450,459],d = 28] [n = 262143,x = [243,246,267,267],d = 24] [n = 917503,x = [467,467,482,499],d = 32] [n = 393215,x = [301,301,318,333],d = 32] [n = 655359,x = [401,401,406,411],d = 10] [n = 131071,x = [169,173,190,191],d = 22] [n = 983039,x = [470,477,513,521],d = 51] [n = 458751,x = [321,323,350,359],d = 38] [n = 720895,x = [409,409,437,442],d = 33] [n = 196607,x = [201,215,234,235],d = 34] [n = 851967,x = [443,449,474,479],d = 36] [n = 327679,x = [275,277,277,314],d = 39] [n = 589823,x = [370,373,387,405],d = 35] [n = 65535,x = [119,121,133,138],d = 19] [n = 491519,x = [335,337,363,366],d = 31] This program was created by Hakan Kjellerstrand, hakank@gmail.com See also my Picat page: http://www.hakank.org/picat/ */ import sat. main => go. % Generate random instances go ?=> nolog, % N=random(1,1000000), % This does not work with fail/0 N = random_cp(1,1000000).first, four_squares(N, X,D), println([n=N,x=X,d=D]), fail, nl. go => true. /* Are there any number which has more than one optimal solution? Yes, quite many: ... [n = 12,d = 2,all = [[1,1,1,3],[0,2,2,2]],len = 2] [n = 19,d = 3,all = [[1,1,1,4],[0,1,3,3]],len = 2] [n = 22,d = 3,all = [[1,1,2,4],[0,2,3,3]],len = 2] [n = 27,d = 3,all = [[1,1,3,4],[0,3,3,3]],len = 2] [n = 28,d = 2,all = [[1,3,3,3],[2,2,2,4]],len = 2] [n = 37,d = 3,all = [[1,2,4,4],[2,2,2,5]],len = 2] [n = 42,d = 3,all = [[2,2,3,5],[1,3,4,4]],len = 2] [n = 48,d = 4,all = [[2,2,2,6],[0,4,4,4]],len = 2] [n = 49,d = 3,all = [[2,2,4,5],[1,4,4,4]],len = 2] [n = 52,d = 2,all = [[3,3,3,5],[2,4,4,4]],len = 2] [n = 60,d = 4,all = [[1,3,5,5],[2,2,4,6]],len = 2] [n = 63,d = 3,all = [[2,3,5,5],[3,3,3,6]],len = 2] [n = 70,d = 3,all = [[3,3,4,6],[2,4,5,5]],len = 2] [n = 71,d = 5,all = [[2,3,3,7],[1,3,5,6]],len = 2] [n = 76,d = 4,all = [[3,3,3,7],[1,5,5,5]],len = 2] [n = 78,d = 5,all = [[2,3,4,7],[1,4,5,6]],len = 2] [n = 79,d = 3,all = [[3,3,5,6],[2,5,5,5]],len = 2] [n = 84,d = 2,all = [[4,4,4,6],[3,5,5,5]],len = 2] [n = 87,d = 5,all = [[1,5,5,6],[2,3,5,7]],len = 2] [n = 88,d = 6,all = [[0,4,6,6],[2,2,4,8]],len = 2] [n = 92,d = 4,all = [[3,3,5,7],[2,4,6,6]],len = 2] [n = 97,d = 3,all = [[3,4,6,6],[4,4,4,7]],len = 2] [n = 98,d = 5,all = [[1,5,6,6],[2,3,6,7],[3,3,4,8]],len = 3] [n = 105,d = 5,all = [[3,4,4,8],[2,4,6,7]],len = 2] [n = 106,d = 3,all = [[4,4,5,7],[3,5,6,6]],len = 2] [n = 112,d = 4,all = [[2,6,6,6],[4,4,4,8]],len = 2] [n = 114,d = 5,all = [[2,5,6,7],[3,4,5,8]],len = 2] [n = 117,d = 3,all = [[4,4,6,7],[3,6,6,6]],len = 2] [n = 118,d = 5,all = [[3,3,6,8],[2,4,7,7]],len = 2] [n = 124,d = 2,all = [[5,5,5,7],[4,6,6,6]],len = 2] [n = 125,d = 5,all = [[2,6,6,7],[3,4,6,8]],len = 2] More than 3 optimal solutions? Again, yes: [n = 252,d = 6,all = [[5,5,9,11],[3,9,9,9],[4,6,10,10],[6,6,6,12]],len = 4] [n = 284,d = 10,all = [[1,9,9,11],[2,6,10,12],[4,6,6,14],[3,5,9,13]],len = 4] [n = 316,d = 6,all = [[5,7,11,11],[4,10,10,10],[6,6,10,12],[7,7,7,13]],len = 4] [n = 388,d = 6,all = [[6,8,12,12],[5,11,11,11],[8,8,8,14],[7,7,11,13]],len = 4] [n = 468,d = 6,all = [[8,8,12,14],[9,9,9,15],[6,12,12,12],[7,9,13,13]],len = 4] [n = 482,d = 8,all = [[5,12,12,13],[7,8,12,15],[6,9,13,14],[8,9,9,16]],len = 4] [n = 556,d = 6,all = [[9,9,13,15],[7,13,13,13],[8,10,14,14],[10,10,10,16]],len = 4] [n = 570,d = 8,all = [[9,10,10,17],[7,10,14,15],[8,9,13,16],[6,13,13,14]],len = 4] [n = 652,d = 6,all = [[9,11,15,15],[8,14,14,14],[11,11,11,17],[10,10,14,16]],len = 4] [n = 666,d = 8,all = [[7,14,14,15],[9,10,14,17],[8,11,15,16],[10,11,11,18]],len = 4] ... More than 4 optimal solutions? Yes, here are some with 5 and even 6 optimal solutions: [n = 6082,d = 15,all = [[34,34,37,49],[28,40,43,43],[33,33,40,48],[29,37,44,44],[31,34,43,46]],len = 5] [n = 6394,d = 15,all = [[32,35,44,47],[30,38,45,45],[35,35,38,50],[29,41,44,44],[34,34,41,49]],len = 5] [n = 6714,d = 15,all = [[31,39,46,46],[35,35,42,50],[36,36,39,51],[30,42,45,45],[33,36,45,48]],len = 5] [n = 7042,d = 15,all = [[31,43,46,46],[36,36,43,51],[37,37,40,52],[34,37,46,49],[32,40,47,47]],len = 5] [n = 7212,d = 22,all = [[26,46,46,48],[32,34,46,54],[34,34,42,56],[27,41,49,49],[29,37,49,51],[35,37,37,57]],len = 6] [n = 7378,d = 15,all = [[33,41,48,48],[32,44,47,47],[38,38,41,53],[35,38,47,50],[37,37,44,52]],len = 5] [n = 7722,d = 15,all = [[38,38,45,53],[39,39,42,54],[33,45,48,48],[34,42,49,49],[36,39,48,51]],len = 5] [n = 8074,d = 15,all = [[34,46,49,49],[37,40,49,52],[40,40,43,55],[35,43,50,50],[39,39,46,54]],len = 5] [n = 8434,d = 15,all = [[40,40,47,55],[41,41,44,56],[35,47,50,50],[38,41,50,53],[36,44,51,51]],len = 5] [n = 8802,d = 15,all = [[37,45,52,52],[41,41,48,56],[39,42,51,54],[42,42,45,57],[36,48,51,51]],len = 5] [n = 8834,d = 17,all = [[41,42,45,58],[38,42,51,55],[35,48,51,52],[36,45,52,53],[40,41,48,57]],len = 5] [n = 9178,d = 15,all = [[40,43,52,55],[37,49,52,52],[38,46,53,53],[42,42,49,57],[43,43,46,58]],len = 5] [n = 9210,d = 17,all = [[41,42,49,58],[37,46,53,54],[39,43,52,56],[42,43,46,59],[36,49,52,53]],len = 5] [n = 9562,d = 15,all = [[44,44,47,59],[39,47,54,54],[41,44,53,56],[43,43,50,58],[38,50,53,53]],len = 5] [n = 9594,d = 17,all = [[38,47,54,55],[43,44,47,60],[42,43,50,59],[40,44,53,57],[37,50,53,54]],len = 5] [n = 9798,d = 17,all = [[43,43,50,60],[42,43,52,59],[37,52,53,54],[44,45,46,61],[38,48,55,55],[39,46,55,56]],len = 6] [n = 9954,d = 15,all = [[40,48,55,55],[44,44,51,59],[42,45,54,57],[39,51,54,54],[45,45,48,60]],len = 5] [n = 9986,d = 17,all = [[43,44,51,60],[41,45,54,58],[44,45,48,61],[39,48,55,56],[38,51,54,55]],len = 5] [n = 10124,d = 22,all = [[37,45,57,59],[40,42,54,62],[35,49,57,57],[42,42,50,64],[43,45,45,65],[34,54,54,56]],len = 6] [n = 10188,d = 24,all = [[38,42,56,62],[37,43,57,61],[41,41,51,65],[42,42,48,66],[34,48,58,58],[33,51,57,57]],len = 6] [n = 10194,d = 17,all = [[44,44,51,61],[43,44,53,60],[45,46,47,62],[39,49,56,56],[40,47,56,57],[38,53,54,55]],len = 6] I haven't yet found any with 7 optimal solutions but there are some with 8 solutions: [n = 81794,d = 29,all = [[127,135,152,156],[125,138,153,154],[122,147,150,151],[129,133,150,158],[124,140,153,153],[134,135,138,163],[132,132,145,161],[131,132,147,160]],len = 8] [n = 82938,d = 29,all = [[126,139,154,155],[130,134,151,159],[125,141,154,154],[135,136,139,164],[132,133,148,161],[123,148,151,152],[128,136,153,157],[133,133,146,162]],len = 8] [n = 84090,d = 29,all = [[133,134,149,162],[136,137,140,165],[124,149,152,153],[131,135,152,160],[127,140,155,156],[134,134,147,163],[126,142,155,155],[129,137,154,158]],len = 8] [n = 85250,d = 29,all = [[125,150,153,154],[130,138,155,159],[128,141,156,157],[127,143,156,156],[135,135,148,164],[132,136,153,161],[134,135,150,163],[137,138,141,166]],len = 8] [n = 86418,d = 29,all = [[138,139,142,167],[136,136,149,165],[129,142,157,158],[131,139,156,160],[126,151,154,155],[133,137,154,162],[128,144,157,157],[135,136,151,164]],len = 8] [n = 87594,d = 29,all = [[130,143,158,159],[129,145,158,158],[139,140,143,168],[134,138,155,163],[127,152,155,156],[136,137,152,165],[132,140,157,161],[137,137,150,166]],len = 8] See go3/0 for more on this... */ go2 ?=> nolog, % member(N,1..1000000), member(N,0..1000000), % N = random_cp(1,1000000).first, four_squares(N, _X,D), % find the optimal solution for this N % Find all optimal solutions All = findall(X2,four_squares(N, X2,D)), % More than one optimal solution? All.len > 1, % All.len > 3, % more than 3 optimal solutions? All.len > 4, % more than 4 optimal solutions? println([n=N,d=D,all=All,len=All.len]), fail, nl. /* Number of optimal solutions for N=1..1000: 1 = 605 2 = 348 3 = 29 4 = 18 Here's the distribution of the number of optimal solutions printed each 1000 N: 1000 = (map)[1 = 605,2 = 348,3 = 29,4 = 18] 2000 = (map)[1 = 1150,2 = 721,3 = 62,4 = 67] 3000 = (map)[1 = 1693,2 = 1076,3 = 111,4 = 120] 4000 = (map)[1 = 2226,2 = 1423,3 = 166,4 = 185] 5000 = (map)[1 = 2757,2 = 1773,3 = 214,4 = 256] 6000 = (map)[1 = 3280,2 = 2126,3 = 262,4 = 332] 7000 = (map)[1 = 3796,2 = 2474,3 = 318,4 = 409,5 = 3] 8000 = (map)[1 = 4313,2 = 2822,3 = 372,4 = 486,5 = 6,6 = 1] 9000 = (map)[1 = 4819,2 = 3183,3 = 425,4 = 562,5 = 10,6 = 1] 10000 = (map)[1 = 5327,2 = 3540,3 = 478,4 = 637,5 = 16,6 = 2] 11000 = (map)[1 = 5831,2 = 3904,3 = 530,4 = 709,5 = 20,6 = 6] 12000 = (map)[1 = 6335,2 = 4262,3 = 586,4 = 783,5 = 24,6 = 10] 13000 = (map)[1 = 6832,2 = 4625,3 = 636,4 = 865,5 = 30,6 = 12] 14000 = (map)[1 = 7338,2 = 4988,3 = 685,4 = 938,5 = 34,6 = 17] 15000 = (map)[1 = 7829,2 = 5363,3 = 735,4 = 1016,5 = 38,6 = 19] 16000 = (map)[1 = 8326,2 = 5732,3 = 783,4 = 1096,5 = 42,6 = 21] 17000 = (map)[1 = 8821,2 = 6097,3 = 831,4 = 1182,5 = 46,6 = 23] 18000 = (map)[1 = 9315,2 = 6466,3 = 878,4 = 1263,5 = 50,6 = 28] 19000 = (map)[1 = 9807,2 = 6837,3 = 924,4 = 1347,5 = 54,6 = 31] 20000 = (map)[1 = 10302,2 = 7208,3 = 975,4 = 1423,5 = 56,6 = 36] 21000 = (map)[1 = 10795,2 = 7578,3 = 1020,4 = 1508,5 = 60,6 = 39] 22000 = (map)[1 = 11284,2 = 7949,3 = 1067,4 = 1592,5 = 64,6 = 44] 23000 = (map)[1 = 11769,2 = 8321,3 = 1110,4 = 1683,5 = 68,6 = 49] 24000 = (map)[1 = 12250,2 = 8705,3 = 1157,4 = 1764,5 = 70,6 = 54] 25000 = (map)[1 = 12739,2 = 9074,3 = 1205,4 = 1849,5 = 75,6 = 58] 26000 = (map)[1 = 13216,2 = 9454,3 = 1253,4 = 1934,5 = 78,6 = 65] 27000 = (map)[1 = 13701,2 = 9832,3 = 1295,4 = 2018,5 = 83,6 = 71] 28000 = (map)[1 = 14182,2 = 10211,3 = 1340,4 = 2102,5 = 86,6 = 79] */ go3 => nolog, Map = new_map(), foreach(N in 1..100000) if four_squares(N, _X,D) then % find the optimal solution for this N All = findall(X2,four_squares(N, X2,D)), Len = All.len, Map.put(Len,Map.get(Len,0)+1) else println(N=no_solution) end, if N > 0, N mod 1000 == 0 then println(N=Map) end end, println(Map), nl. /* What is the distribution of the optimal value (D)? For N=1..1000 0 = 15 1 = 48 2 = 76 3 = 102 4 = 144 5 = 148 6 = 137 7 = 135 8 = 90 9 = 54 10 = 25 12 = 18 16 = 6 20 = 1 24 = 1 */ go4 => nolog, Map = new_map(), foreach(N in 1..1000) four_squares(N, _X,D), Map.put(D,Map.get(D,0)+1) end, println(Map), nl. four_squares(N, X,D) => X=new_list(4), X :: 0..ceiling(sqrt(N)), increasing(X), N #= sum([X[I]**2: I in 1..4]), D #= max(X)-min(X), if var(D) then solve($[min(D),degree,updown],X) else % Generate all optimal solutions solve(X) end. random_cp(N,MaxVal) = X => X = new_array(N), X :: 1..MaxVal, solve($[rand],X).