Uniquely separated pixels
Python 3, 135 136 137
10 6830 20470 47750 370770 148190 306910 373250 267230 354030 30390 361470 118430 58910 197790 348450 381336 21710 183530 305050 2430 1810 365832 99038 381324 39598 262270 365886 341662 15478 9822 365950 44526 58862 24142 381150 31662 237614 118830 380846 7182 113598 306750 11950 373774 111326 272358 64310 43990 200278 381014 165310 254454 12394 382534 87894 6142 750 382478 15982 298326 70142 186478 152126 367166 1162 23426 341074 7306 76210 140770 163410 211106 207962 35282 165266 300178 120106 336110 30958 158 362758 382894 308754 88434 336918 244502 43502 54990 279910 175966 234054 196910 287284 288468 119040 275084 321268 17968 2332 86064 340044 244604 262436 111188 291868 367695 362739 370781 375723 360261 377565 383109 328689 347879 2415 319421 55707 352897 313831 302079 19051 346775 361293 328481 35445 113997 108547 309243 19439 199037 216463 62273 174471 207197 167695 296927
Found using a greedy algorithm which, at each stage, chooses the valid pixel whose set of distances to the chosen pixels overlaps the least with that of other pixels.
Specifically, the scoring is
score(P) = sum(number of pixels with D in its distance set
for each D in P's distance set)
and the pixel with the lowest score is chosen.
The search is kicked off with the point 10 (i.e. (0, 10)). This part is adjustable, so beginning with different pixels may lead to better or worse results.
It's quite a slow algorithm, so I'm trying to add optimisations/heuristics, and maybe some backtracking. PyPy is recommended for speed.
Anyone trying to come up with an algorithm should test on N = 10, for which I've got 9 (but this took a lot of tweaking and trying different initial points):
Code
from collections import Counter, defaultdict
import sys
import time
N = 619
start_time = time.time()
def norm(p1, p2):
return (p1//N - p2//N)**2 + (p1%N - p2%N)**2
selected = [10]
selected_dists = {norm(p1, p2) for p1 in selected for p2 in selected if p1 != p2}
pix2dist = {} # {candidate pixel: {distances to chosen}}
dist2pix = defaultdict(set)
for pixel in range(N*N):
if pixel in selected:
continue
dist_list = [norm(pixel, p) for p in selected]
dist_set = set(dist_list)
if len(dist_set) != len(dist_list) or dist_set & selected_dists:
continue
pix2dist[pixel] = dist_set
for dist in dist_set:
dist2pix[dist].add(pixel)
while pix2dist:
best_score = None
best_pixel = None
for pixel in sorted(pix2dist): # Sorting for determinism
score = sum(len(dist2pix[d]) for d in pix2dist[pixel])
if best_score is None or score < best_score:
best_score = score
best_pixel = pixel
added_dists = pix2dist[best_pixel]
selected_dists |= added_dists
del pix2dist[best_pixel]
selected.append(best_pixel)
for d in added_dists:
dist2pix[d].remove(best_pixel)
to_remove = set()
for pixel in pix2dist:
new_dist = norm(pixel, best_pixel)
if (new_dist in selected_dists or new_dist in pix2dist[pixel]
or added_dists & pix2dist[pixel]):
to_remove.add(pixel)
continue
pix2dist[pixel].add(new_dist)
dist2pix[new_dist].add(pixel)
for pixel in to_remove:
for d in pix2dist[pixel]:
dist2pix[d].remove(pixel)
del pix2dist[pixel]
print("Selected: {}, Remaining: {}, Chosen: ({}, {})".format(len(selected), len(pix2dist),
best_pixel//N, best_pixel%N))
sys.stdout.flush()
print(*selected)
print("Time taken:", time.time() - start_time)
SWI-Prolog, score 131
Barely better than the initial answer, but I guess this will get things started a bit more. The algorithm is the same as the Python answer except for the fact that it tries pixels in an alternate way, starting with the top left pixel (pixel 0), then the bottom right pixel (pixel 383160), then pixel 1, then pixel 383159, etc.
a(Z) :-
N = 619,
build_list(N,Z).
build_list(N,R) :-
M is N*N,
get_list([M,-1],[],L),
reverse(L,O),
build_list(N,O,[],[],R).
get_list([A,B|C],R,Z) :-
X is A - 1,
Y is B + 1,
(X =< Y,
Z = R
;
get_list([X,Y,A,B|C],[X,Y|R],Z)).
build_list(_,[],R,_,R) :- !.
build_list(N,[A|T],R,W,Z) :-
separated_pixel(N,A,R,W,S),
is_set(S),
flatten([W|S],V),!,
build_list(N,T,[A|R],V,Z)
;build_list(N,T,R,W,Z).
separated_pixel(N,A,L,W,R) :-
separated_pixel(N,A,L,[],W,R).
separated_pixel(N,A,[A|T],R,W,S) :-
separated_pixel(N,A,T,R,W,S).
separated_pixel(N,A,[B|T],R,W,S) :-
X is (A mod N - B mod N)*(A mod N - B mod N),
Y is (A//N - B//N)*(A//N - B//N),
Z is X + Y,
\+member(Z,W),
separated_pixel(N,A,T,[Z|R],W,S).
separated_pixel(_,_,[],R,_,R).
Input:
a(A).
Output:
Z = [202089, 180052, 170398, 166825, 235399, 138306, 126354, 261759, 119490, 117393, 281623, 95521, 290446, 299681, 304310, 78491, 314776, 63618, 321423, 60433, 323679, 52092, 331836, 335753, 46989, 40402, 343753, 345805, 36352, 350309, 32701, 32470, 352329, 30256, 28089, 357859, 23290, 360097, 22534, 362132, 20985, 364217, 365098, 17311, 365995, 15965, 15156, 368487, 370980, 371251, 11713, 372078, 372337, 10316, 373699, 8893, 374417, 8313, 7849, 7586, 7289, 6922, 376588, 6121, 5831, 377399, 377639, 4941, 378494, 4490, 379179, 3848, 379453, 3521, 3420, 379963, 380033, 3017, 380409, 2579, 380636, 2450, 2221, 2006, 381235, 1875, 381369, 381442, 381682, 1422, 381784, 1268, 381918, 1087, 382144, 382260, 833, 382399, 697, 382520, 622, 382584, 382647, 382772, 384, 382806, 319, 286, 382915, 382939, 190, 172, 383005, 128, 383050, 93, 383076, 68, 383099, 52, 40, 383131, 21, 383145, 10, 383153, 4, 383158, 1, 383160, 0]
Image from Stack Snippet
Haskell—115 130 131 135 136
My inspiration was the Sieve of Eratosthenes and in particular The Genuine Sieve of Eratosthenes, a paper by Melissa E. O'Neill of Harvey Mudd College. My original version (which considered points in index order) sieved points extremely quickly, for some reason I can't recall I decided to shuffle the points before “sieving” them in this version (I think solely to make generating different answers easier by using a new seed in the random generator). Because the points are no longer in any sort of order there's not really any sieving going on anymore, and as a result it takes a couple minutes just to produce this single 115 point answer. A knockout Vector would probably be a better choice now.
So with this version as a checkpoint I see two branches, returning to the “Genuine Sieve” algorithm and leveraging the List monad for choice, or swapping out the Set operations for equivalents on Vector.
Edit: So for working version two I veered back toward the sieve algorithm, improved the generation of “multiples” (knocking indices out by finding points at integer coordinates on circles with radius equal to the distance between any two points, akin to generating prime multiples) and making a few constant time improvements by avoiding some needless recalculation.
For some reason I can't recompile with profiling turned on, but I believe that the major bottleneck now is backtracking. I think exploring a bit of parallelism and concurrency will produce linear speedups, but memory exhaustion will probably cap me at a 2x improvement.
Edit: Version 3 meandered a bit, I first experimented with a heuristic in taking the next 𝐧 indices (after sieving from earlier choices) and choosing the one that produced the next minimum knockout set. This ended up being far too slow, so I went back to a whole-searchspace brute force method. An idea to order the points by distance from some origin came to me, and led to an improvement by one single point (in the time my patience lasted). This version picks index 0 as the origin, it may be worth trying the center point of the plane.
Edit: I picked up 4 points by re-ordering the search space to prioritize the most distant points from the center. If you're testing my code, 135 136 is actually the second third solution found. Fast edit: This version seems the most likely to keep being productive if left running. I suspect I may tie at 137, then run out of patience waiting for 138.
One thing I noticed (that may be of help to someone) is that if you set the point ordering from the center of the plane (i.e., remove (d*d -) from originDistance) the image formed looks a bit like a sparse prime spiral.
{-# LANGUAGE RecordWildCards #-}
{-# LANGUAGE BangPatterns #-}
module Main where
import Data.Function (on)
import Data.List (tails, sortBy)
import Data.Maybe (fromJust)
import Data.Ratio
import Data.Set (fromList, toList, union, difference, member)
import System.IO
sideLength :: Int
sideLength = 619
data Point = Point { x :: !Int, y :: !Int } deriving (Ord, Eq)
data Delta = Delta { da :: !Int, db :: !Int }
euclidean :: Delta -> Int
euclidean Delta{..} = da*da + db*db
instance Eq Delta where
(==) = (==) `on` euclidean
instance Ord Delta where
compare = compare `on` euclidean
delta :: Point -> Point -> Delta
delta a b = Delta (min dx dy) (max dx dy)
where
dx = abs (x a - x b)
dy = abs (y a - y b)
equidistant :: Dimension -> Point -> Point -> [Point]
equidistant d a b =
let
(dx, dy) = (x a - x b, y a - y b)
m = if dx == 0 then Nothing else Just (dy % dx) -- Slope
w = if dy == 0 then Nothing else Just $ maybe 0 (negate . recip) m -- Negative reciprocal
justW = fromJust w -- Moral bankruptcy
(px, py) = ((x a + x b) % 2, (y a + y b) % 2) -- Midpoint
b0 = py - (justW * px) -- Y-intercept
f q = justW * q + b0 -- Perpendicular bisector
in
maybe (if denominator px == 1 then map (Point (numerator px)) [0..d - 1] else [])
( map (\q -> Point q (numerator . f . fromIntegral $ q))
. filter ((== 1) . denominator . f . fromIntegral)
)
(w >> return [0..d - 1])
circle :: Dimension -> Point -> Delta -> [Point]
circle d p delta' =
let
square = (^(2 :: Int))
hypoteneuse = euclidean delta'
candidates = takeWhile ((<= hypoteneuse) . square) [0..d - 1]
candidatesSet = fromList $ map square [0..d - 1]
legs = filter ((`member` candidatesSet) . (hypoteneuse -) . square) candidates
pythagoreans = zipWith Delta legs
$ map (\l -> floor . sqrt . (fromIntegral :: Int -> Double) $ hypoteneuse - square l) legs
in
toList . fromList $ concatMap (knight p) pythagoreans
knight :: Point -> Delta -> [Point]
knight Point{..} Delta{..} =
[ Point (x + da) (y - db), Point (x + da) (y + db)
, Point (x + db) (y - da), Point (x + db) (y + da)
, Point (x - da) (y - db), Point (x - da) (y + db)
, Point (x - db) (y - da), Point (x - db) (y + da)
]
type Dimension = Int
type Index = Int
index :: Dimension -> Point -> Index
index d Point{..} = y * d + x
point :: Dimension -> Index -> Point
point d i = Point (i `rem` d) (i `div` d)
valid :: Dimension -> Point -> Bool
valid d Point{..} = 0 <= x && x < d
&& 0 <= y && y < d
isLT :: Ordering -> Bool
isLT LT = True
isLT _ = False
sieve :: Dimension -> [[Point]]
sieve d = [i0 : sieve' is0 [i0] [] | (i0:is0) <- tails . sortBy originDistance . map (point d) $ [0..d*d - 1]]
where
originDistance :: Point -> Point -> Ordering
originDistance = compare `on` ((d*d -) . euclidean . delta (point d (d*d `div` 2)))
sieve' :: [Point] -> [Point] -> [Delta] -> [Point]
sieve' [] _ _ = []
sieve' (i:is) ps ds = i : sieve' is' (i:ps) ds'
where
ds' = map (delta i) ps ++ ds
knockouts = fromList [k | d' <- ds
, k <- circle d i d'
, valid d k
, not . isLT $ k `originDistance` i
]
`union` fromList [k | q <- i : ps
, d' <- map (delta i) ps
, k <- circle d q d'
, valid d k
, not . isLT $ k `originDistance` i
]
`union` fromList [e | q <- ps
, e <- equidistant d i q
, valid d e
, not . isLT $ e `originDistance` i
]
is' = sortBy originDistance . toList $ fromList is `difference` knockouts
main :: IO ()
main = do let answers = strictlyIncreasingLength . map (map (index sideLength)) $ sieve sideLength
hSetBuffering stdout LineBuffering
mapM_ (putStrLn . unwords . map show) $ answers
where
strictlyIncreasingLength :: [[a]] -> [[a]]
strictlyIncreasingLength = go 0
where
go _ [] = []
go n (x:xs) = if n < length x then x : go (length x) xs else go n xs
Output
1237 381923 382543 382541 1238 1857 380066 5 380687 378828 611 5571 382553 377587 375113 3705 8664 376356 602 1253 381942 370161 12376 15475 7413 383131 367691 380092 376373 362114 36 4921 368291 19180 382503 26617 3052 359029 353451 29716 382596 372674 352203 8091 25395 12959 382479 381987 35894 346031 1166 371346 336118 48276 2555 332400 46433 29675 380597 13066 382019 1138 339859 368230 29142 58174 315070 326847 56345 337940 2590 382663 320627 70553 19278 7309 82942 84804 64399 5707 461 286598 363864 292161 89126 371267 377122 270502 109556 263694 43864 382957 824 303886 248218 18417 347372 282290 144227 354820 382909 380301 382808 334361 375341 2197 260623 222212 196214 231526 177637 29884 251280 366739 39442 143568 132420 334718 160894 353132 78125 306866 140600 297272 54150 240054 98840 219257 189278 94968 226987 265881 180959 142006 218763 214475