fork download
copy 
  1. import numpy as np
  2. from typing import List, Set, Union
  3.  
  4.  
  5. def _make_rng(rng: Union[None, int, np.random.Generator]) -> np.random.Generator:
  6.     '''Create a deterministic RNG from None, an int seed or a Generator.'''
  7.     if rng is None:
  8.         return np.random.default_rng(0)
  9.     if isinstance(rng, np.random.Generator):
  10.         # copy state to avoid side‑effects
  11.         state = rng.bit_generator.state
  12.         return np.random.Generator(np.random.PCG64(state))
  13.     return np.random.default_rng(int(rng))
  14.  
  15.  
  16. def _conflicts(x: int, S: Set[int]) -> bool:
  17.     '''True iff adding x to S creates a 3‑AP with two elements already in S.'''
  18.     for y in S:
  19.         # x as middle term
  20.         if (x + y) % 2 == 0 and ((x + y) // 2) in S:
  21.             return True
  22.         # x as an endpoint
  23.         if (2 * x - y) in S or (2 * y - x) in S:
  24.             return True
  25.     return False
  26.  
  27.  
  28. def heuristic(n: int, rng: Union[None, int, np.random.Generator] = None) -> List[int]:
  29.     '''
  30.     Return a large 3‑AP‑free subset of {1,…,n}.
  31.     The algorithm performs a greedy randomized construction and a short
  32.     random‑swap local improvement phase.
  33.     '''
  34.     if n <= 0:
  35.         return []
  36.  
  37.     gen = _make_rng(rng)
  38.  
  39.     # ---------- 1. Greedy randomized construction ----------
  40.     all_numbers = np.arange(1, n + 1, dtype=int)
  41.     gen.shuffle(all_numbers)                     # random order, deterministic
  42.     S: Set[int] = set()
  43.     omitted: Set[int] = set()
  44.  
  45.     for x in all_numbers:
  46.         if not _conflicts(int(x), S):
  47.             S.add(int(x))
  48.         else:
  49.             omitted.add(int(x))
  50.  
  51.     # ---------- 2. Local improvement (random swaps) ----------
  52.     # Budget proportional to sqrt(n) * n to stay fast for a few thousand.
  53.     budget = int(3 * n * np.sqrt(n))
  54.     temperature_start = 1.0
  55.     temperature_end = 0.001
  56.  
  57.     omitted = set(range(1, n + 1)) - S  # recompute after greedy phase
  58.  
  59.     for step in range(budget):
  60.         # linear cooling
  61.         T = temperature_start + (temperature_end - temperature_start) * (step / budget)
  62.  
  63.         if not omitted:
  64.             break
  65.  
  66.         # pick a random element not in S
  67.         x = int(gen.choice(list(omitted)))
  68.  
  69.         # collect elements of S that would conflict with x
  70.         conflicting: List[int] = []
  71.         for y in S:
  72.             if (x + y) % 2 == 0 and ((x + y) // 2) in S:
  73.                 conflicting.append(y)
  74.             elif (2 * x - y) in S or (2 * y - x) in S:
  75.                 conflicting.append(y)
  76.  
  77.         # If no conflict we could simply add, but that was already checked.
  78.         # Try a replace‑one‑with‑up‑to‑two move.
  79.         if not conflicting:
  80.             continue
  81.  
  82.         # remove one conflicting element (randomly)
  83.         to_remove = int(gen.choice(conflicting))
  84.         S.remove(to_remove)
  85.         omitted.add(to_remove)
  86.  
  87.         # attempt to insert x and possibly another non‑conflicting omitted element
  88.         added: List[int] = []
  89.         candidates = [x] + list(omitted)
  90.         gen.shuffle(candidates)
  91.  
  92.         for y in candidates:
  93.             if y == to_remove:
  94.                 continue
  95.             if not _conflicts(y, S):
  96.                 S.add(y)
  97.                 omitted.remove(y)
  98.                 added.append(y)
  99.                 if len(added) == 2:
  100.                     break
  101.  
  102.         net_gain = len(added) - 1   # removed one, added len(added)
  103.  
  104.         # Acceptance test (accept if gain > 0, else Metropolis)
  105.         if net_gain > 0 or gen.random() < np.exp(net_gain / max(T, 1e-12)):
  106.             # keep the move (already applied)
  107.             pass
  108.         else:
  109.             # revert
  110.             for y in added:
  111.                 S.remove(y)
  112.                 omitted.add(y)
  113.             S.add(to_remove)
  114.             omitted.remove(to_remove)
  115.  
  116.     # Return as a sorted list for reproducibility
  117.     return S
  118. for i in range(10):
  119. 	s = heuristic(82, i)
  120. 	print(s, "|s|", len(s))
Success #stdin #stdout 3.89s 41248KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
{1, 3, 7, 6, 10, 16, 20, 21, 23, 27, 28, 52, 54, 58, 59, 61, 67, 71, 72, 74, 78, 79} |s| 22
{1, 2, 4, 8, 9, 11, 19, 26, 22, 23, 28, 31, 56, 57, 59, 66, 67, 71, 72, 80, 79, 74} |s| 22
{5, 1, 4, 8, 10, 18, 24, 23, 20, 27, 29, 56, 58, 59, 61, 67, 72, 71, 74, 80, 82, 79} |s| 22
{1, 3, 6, 7, 10, 15, 21, 22, 25, 31, 33, 48, 54, 53, 57, 67, 70, 69, 72, 79, 78, 82} |s| 22
{1, 4, 3, 11, 9, 20, 22, 23, 28, 27, 60, 46, 48, 49, 62, 59, 67, 66, 77, 79, 80, 82} |s| 22
{3, 56, 7, 10, 27, 18, 20, 21, 62, 6, 25, 28, 57, 1, 60, 73, 71, 66, 74, 79, 81, 78} |s| 22
{1, 2, 8, 10, 7, 11, 22, 23, 25, 26, 32, 31, 47, 58, 59, 64, 66, 67, 77, 79, 80, 82} |s| 22
{1, 2, 4, 8, 10, 13, 21, 23, 26, 27, 30, 49, 55, 57, 62, 63, 70, 73, 74, 79, 81, 82} |s| 22
{5, 4, 10, 8, 1, 21, 22, 27, 26, 29, 35, 56, 55, 59, 64, 61, 71, 74, 80, 76, 82} |s| 21
{1, 2, 4, 8, 9, 11, 19, 22, 23, 26, 28, 31, 46, 56, 57, 63, 65, 72, 71, 82, 76, 78} |s| 22
https://ideone.com/RkklDD
language:
Python 3 (python  3.12)
created:
1 day ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!