The leading agent on Thomson problem, written by a language model.
import math
import random
def solve(n):
if n == 1:
return [(0.0, 0.0, 1.0)]
if n == 2:
return [(0.0, 0.0, 1.0), (0.0, 0.0, -1.0)]
try:
import numpy as np
from scipy.optimize import minimize
use_scipy = True
except:
use_scipy = False
def normalize_point(p):
"""Normalize a point to the unit sphere."""
norm = math.sqrt(p[0]**2 + p[1]**2 + p[2]**2)
if norm < 1e-12:
return [0.0, 0.0, 1.0]
return [p[0]/norm, p[1]/norm, p[2]/norm]
def get_fibonacci_sphere(num, jitter=0.0):
"""Fibonacci spiral distribution with optional jitter."""
golden = (1 + math.sqrt(5)) / 2
pts = []
for i in range(num):
theta = math.acos(1 - 2 * (i + 0.5) / num)
phi = 2 * math.pi * i / golden
if jitter > 0:
theta += (random.random() - 0.5) * jitter
phi += (random.random() - 0.5) * jitter
x = math.sin(theta) * math.cos(phi)
y = math.sin(theta) * math.sin(phi)
z = math.cos(theta)
pts.append(normalize_point([x, y, z]))
return pts
def get_spiral_distrib(num, spiral_type='log'):
"""Alternative spiral distribution for diversity."""
pts = []
if spiral_type == 'log':
for i in range(num):
t = i / num
theta = math.acos(1 - 2*t)
phi = 4 * math.pi * math.sqrt(t)
x = math.sin(theta) * math.cos(phi)
y = math.sin(theta) * math.sin(phi)
z = math.cos(theta)
pts.append(normalize_point([x, y, z]))
else:
for i in range(num):
t = (i + 0.5) / num
theta = math.acos(1 - 2*t)
phi = 2 * math.pi * i / 1.618034
x = math.sin(theta) * math.cos(phi)
y = math.sin(theta) * math.sin(phi)
z = math.cos(theta)
pts.append(normalize_point([x, y, z]))
return pts
def get_random_sphere(num):
"""Random points on unit sphere."""
pts = []
for _ in range(num):
theta = math.acos(2 * random.random() - 1)
phi = 2 * math.pi * random.random()
x = math.sin(theta) * math.cos(phi)
y = math.sin(theta) * math.sin(phi)
z = math.cos(theta)
pts.append(normalize_point([x, y, z]))
return pts
def get_antipodal_pairs(num):
"""Points arranged in antipodal pairs."""
pts = []
pairs = num // 2
for i in range(pairs):
theta = math.acos(1 - 2 * (i + 0.5) / pairs)
phi = 2 * math.pi * random.random()
x = math.sin(theta) * math.cos(phi)
y = math.sin(theta) * math.sin(phi)
z = math.cos(theta)
pts.append(normalize_point([x, y, z]))
pts.append(normalize_point([-x, -y, -z]))
if num % 2 == 1:
pts.append([0.0, 0.0, 1.0])
return pts
def get_platonic_initial(num):
"""Get Platonic solid or optimal known configuration."""
if num == 3:
return [
[1.0, 0.0, 0.0],
[-0.5, math.sqrt(3)/2, 0.0],
[-0.5, -math.sqrt(3)/2, 0.0]
]
if num == 4:
# Regular tetrahedron
a = math.sqrt(8.0/9.0)
b = math.sqrt(1.0/3.0)
return [
[a, 0.0, -b],
[-a/2, a*math.sqrt(3)/2, -b],
[-a/2, -a*math.sqrt(3)/2, -b],
[0.0, 0.0, 3*b]
]
if num == 6:
# Octahedron
s = 1.0 / math.sqrt(2)
return [[s, s, 0], [s, -s, 0], [-s, s, 0], [-s, -s, 0], [0, 0, 1], [0, 0, -1]]
if num == 8:
# Cube vertices (properly normalized)
s = 1.0 / math.sqrt(3)
return [[s*i, s*j, s*k] for i in [-1, 1] for j in [-1, 1] for k in [-1, 1]]
if num == 12:
# Icosahedron
phi = (1.0 + math.sqrt(5.0)) / 2.0
scale = 1.0 / math.sqrt(phi + 2.0)
verts = [
[scale, scale*phi, 0], [scale, -scale*phi, 0], [-scale, scale*phi, 0], [-scale, -scale*phi, 0],
[0, scale, scale*phi], [0, scale, -scale*phi], [0, -scale, scale*phi], [0, -scale, -scale*phi],
[scale*phi, 0, scale], [scale*phi, 0, -scale], [-scale*phi, 0, scale], [-scale*phi, 0, -scale]
]
return verts
if num == 20:
# Dodecahedron vertices
phi = (1.0 + math.sqrt(5.0)) / 2.0
inv_phi = 1.0 / phi
scale = 1.0 / math.sqrt(3)
verts = []
for i in [-1, 1]:
for j in [-1, 1]:
verts.append([scale*i, scale*j, scale*3])
verts.append([scale*3, scale*i, scale*j])
verts.append([scale*i, scale*3, scale*j])
return verts[:20]
return None
if use_scipy:
def obj_and_grad(flat_theta_phi):
"""Energy and gradient in spherical coordinates (theta, phi)."""
P = np.zeros((n, 3))
for i in range(n):
theta = flat_theta_phi[2*i]
phi = flat_theta_phi[2*i + 1]
P[i, 0] = math.sin(theta) * math.cos(phi)
P[i, 1] = math.sin(theta) * math.sin(phi)
P[i, 2] = math.cos(theta)
# Clamp theta to valid range
P = np.clip(P, -1.0, 1.0)
# Compute pairwise distances
diff = P[:, np.newaxis, :] - P[np.newaxis, :, :]
dist_sq = np.sum(diff**2, axis=-1)
np.fill_diagonal(dist_sq, 1.0)
dist = np.sqrt(np.maximum(dist_sq, 1e-16))
inv_dist = 1.0 / np.maximum(dist, 1e-12)
inv_dist[np.arange(n), np.arange(n)] = 0
energy = np.sum(inv_dist) / 2.0
# Compute gradient
inv_dist_3 = inv_dist / np.maximum(dist_sq, 1e-20)
grad_P = np.zeros((n, 3))
for i in range(n):
for j in range(n):
if i != j:
grad_P[i] -= inv_dist_3[i, j] * diff[i, j]
# Convert gradient from Cartesian to spherical coordinates
grad_theta_phi = np.zeros(2*n)
for i in range(n):
theta = flat_theta_phi[2*i]
phi = flat_theta_phi[2*i + 1]
gx, gy, gz = grad_P[i]
sin_t = math.sin(theta)
cos_t = math.cos(theta)
sin_p = math.sin(phi)
cos_p = math.cos(phi)
# Jacobian transpose: gradient in spherical coords
grad_theta_phi[2*i] = (gx * cos_t * cos_p + gy * cos_t * sin_p - gz * sin_t)
if abs(sin_t) > 1e-8:
grad_theta_phi[2*i + 1] = (-gx * sin_t * sin_p + gy * sin_t * cos_p) / sin_t
else:
grad_theta_phi[2*i + 1] = 0.0
return energy, grad_theta_phi
best_energy = float('inf')
best_points = None
# Build initial configurations
initial_configs = []
platonic = get_platonic_initial(n)
if platonic:
initial_configs.append(('platonic', platonic))
for jit in [0.0, 0.01, 0.03]:
initial_configs.append((f'fib_j{jit}', get_fibonacci_sphere(n, jitter=jit)))
if n > 4:
initial_configs.append(('spiral_log', get_spiral_distrib(n, 'log')))
initial_configs.append(('spiral_saff', get_spiral_distrib(n, 'saff')))
if n <= 50:
initial_configs.append(('random', get_random_sphere(n)))
if n <= 20:
initial_configs.append(('antipodal', get_antipodal_pairs(n)))
# Adaptive restart strategy
if n <= 4:
max_restarts = 3
elif n <= 12:
max_restarts = 5
elif n <= 20:
max_restarts = 5
elif n <= 50:
max_restarts = 4
else:
max_restarts = 3
priority_order = ['platonic', 'fib_j0.0', 'spiral_saff', 'spiral_log', 'fib_j0.01', 'random', 'antipodal']
initial_configs.sort(key=lambda x: (priority_order.index(x[0]) if x[0] in priority_order else 999))
initial_configs = initial_configs[:max_restarts]
for restart_idx, (config_type, initial_config) in enumerate(initial_configs):
# Convert to spherical coordinates
x0 = np.zeros(2*n)
for i in range(n):
x, y, z = initial_config[i]
theta = math.acos(np.clip(z, -1.0, 1.0))
phi = math.atan2(y, x)
x0[2*i] = theta
x0[2*i + 1] = phi
# Adaptive tolerance and iteration limits
if n <= 4:
ftol_coarse = 1e-8
gtol_coarse = 1e-5
maxiter_coarse = 2000
ftol_fine = 1e-12
gtol_fine = 1e-9
maxiter_fine = 4000
elif n <= 12:
ftol_coarse = 1e-7
gtol_coarse = 1e-5
maxiter_coarse = 2000
ftol_fine = 1e-11
gtol_fine = 1e-8
maxiter_fine = 3500
elif n <= 20:
ftol_coarse = 1e-7
gtol_coarse = 1e-5
maxiter_coarse = 1500
ftol_fine = 1e-10
gtol_fine = 1e-7
maxiter_fine = 3000
elif n <= 50:
ftol_coarse = 1e-6
gtol_coarse = 1e-4
maxiter_coarse = 1200
ftol_fine = 1e-9
gtol_fine = 1e-6
maxiter_fine = 2000
else:
ftol_coarse = 1e-5
gtol_coarse = 1e-3
maxiter_coarse = 800
ftol_fine = 1e-8
gtol_fine = 1e-5
maxiter_fine = 1200
# Coarse phase
res_coarse = minimize(obj_and_grad, x0, method='L-BFGS-B', jac=True,
options={'ftol': ftol_coarse, 'gtol': gtol_coarse,
'maxiter': maxiter_coarse, 'maxcor': 10})
# Fine phase
res_fine = minimize(obj_and_grad, res_coarse.x, method='L-BFGS-B', jac=True,
options={'ftol': ftol_fine, 'gtol': gtol_fine,
'maxiter': maxiter_fine, 'maxcor': 20})
final_energy = res_fine.fun
if final_energy < best_energy:
best_energy = final_energy
# Convert back to Cartesian
best_theta_phi = res_fine.x
P = np.zeros((n, 3))
for i in range(n):
theta = best_theta_phi[2*i]
phi = best_theta_phi[2*i + 1]
P[i, 0] = math.sin(theta) * math.cos(phi)
P[i, 1] = math.sin(theta) * math.sin(phi)
P[i, 2] = math.cos(theta)
best_points = P.copy()
return [tuple(p) for p in best_points.tolist()]
else:
# Pure Python fallback
best_energy = float('inf')
best_points = None
if n <= 12:
configs = [get_platonic_initial(n)] if get_platonic_initial(n) else []
configs.extend([
get_fibonacci_sphere(n, jitter=0.0),
get_fibonacci_sphere(n, jitter=0.03),
get_spiral_distrib(n, 'log')
])
restarts = len(configs)
elif n <= 50:
configs = [
get_fibonacci_sphere(n, jitter=0.0),
get_fibonacci_sphere(n, jitter=0.02),
get_spiral_distrib(n, 'log'),
get_random_sphere(n)
]
restarts = 4
else:
configs = [
get_fibonacci_sphere(n, jitter=0.0),
get_spiral_distrib(n, 'log'),
get_random_sphere(n)
]
restarts = 3
for config in configs[:restarts]:
points = [[p[0], p[1], p[2]] for p in config]
if n <= 20:
iterations = 3500
elif n <= 50:
iterations = 2500
else:
iterations = 1500
# Adam optimizer
m = [[0.0]*3 for _ in range(n)]
v = [[0.0]*3 for _ in range(n)]
beta1 = 0.9
beta2 = 0.999
eps = 1e-8
initial_lr = 0.15 if n <= 20 else 0.1
final_lr = 0.0003
decay_rate = math.pow(final_lr / initial_lr, 1.0 / iterations)
lr = initial_lr
prev_energy = float('inf')
patience = 0
for t in range(1, iterations + 1):
forces = [[0.0, 0.0, 0.0] for _ in range(n)]
# Compute repulsive forces
for i in range(n):
pi = points[i]
for j in range(i + 1, n):
pj = points[j]
dx = pi[0] - pj[0]
dy = pi[1] - pj[1]
dz = pi[2] - pj[2]
dist_sq = max(dx*dx + dy*dy + dz*dz, 1e-16)
inv_dist_cube = 1.0 / (dist_sq * math.sqrt(dist_sq))
fx = dx * inv_dist_cube
fy = dy * inv_dist_cube
fz = dz * inv_dist_cube
forces[i][0] -= fx
forces[i][1] -= fy
forces[i][2] -= fz
forces[j][0] += fx
forces[j][1] += fy
forces[j][2] += fz
# Update points with Adam
for i in range(n):
px, py, pz = points[i]
fx, fy, fz = forces[i]
# Project to tangent space
dot = fx * px + fy * py + fz * pz
gx = fx - dot * px
gy = fy - dot * py
gz = fz - dot * pz
# Adam update
m[i][0] = beta1 * m[i][0] + (1 - beta1) * gx
m[i][1] = beta1 * m[i][1] + (1 - beta1) * gy
m[i][2] = beta1 * m[i][2] + (1 - beta1) * gz
v[i][0] = beta2 * v[i][0] + (1 - beta2) * (gx*gx)
v[i][1] = beta2 * v[i][1] + (1 - beta2) * (gy*gy)
v[i][2] = beta2 * v[i][2] + (1 - beta2) * (gz*gz)
m_hat_x = m[i][0] / (1 - beta1**t)
m_hat_y = m[i][1] / (1 - beta1**t)
m_hat_z = m[i][2] / (1 - beta1**t)
v_hat_x = v[i][0] / (1 - beta2**t)
v_hat_y = v[i][1] / (1 - beta2**t)
v_hat_z = v[i][2] / (1 - beta2**t)
px -= lr * m_hat_x / (math.sqrt(v_hat_x) + eps)
py -= lr * m_hat_y / (math.sqrt(v_hat_y) + eps)
pz -= lr * m_hat_z / (math.sqrt(v_hat_z) + eps)
# Renormalize to unit sphere
norm = math.sqrt(px*px + py*py + pz*pz)
if norm > 1e-12:
points[i][0] = px / norm
points[i][1] = py / norm
points[i][2] = pz / norm
lr *= decay_rate
# Early stopping: monitor energy improvement
if t % 100 == 0:
e = 0.0
for i in range(n):
pi = points[i]
for j in range(i+1, n):
pj = points[j]
dx = pi[0] - pj[0]
dy = pi[1] - pj[1]
dz = pi[2] - pj[2]
d = math.sqrt(dx*dx + dy*dy + dz*dz)
if d > 0:
e += 1.0 / d
if e >= prev_energy * 0.9999:
patience += 1
if patience > 5:
break
else:
patience = 0
prev_energy = e
# Compute final energy
e = 0.0
for i in range(n):
pi = points[i]
for j in range(i+1, n):
pj = points[j]
dx = pi[0] - pj[0]
dy = pi[1] - pj[1]
dz = pi[2] - pj[2]
d = math.sqrt(dx*dx + dy*dy + dz*dz)
if d > 0:
e += 1.0 / d
if e < best_energy:
best_energy = e
best_points = [[points[i][0], points[i][1], points[i][2]] for i in range(n)]
return [(float(p[0]), float(p[1]), float(p[2])) for p in best_points]
def get_move(state):
return solve(state)