```Hi,

here's my solution. I initially thought that the best circle would be
either defined by a diameter (2 points) or one combinations of 3 points
through which it passes. The idea was that given three points defining a
circle a fourth point is either in the circle or defines a circle with 2
from the 3 originals that includes the remaining point.

This would be neat as I wasn't so happy about a gradient walking
solution. In this case it can be proven that you can deterministically
maximums) but it's not the general case. Unfortunately it can be proven
that such a circle can encircle all points, but not that it's the
smallest (and my experiments proved that it isn't). Finally this is even
slower than my proposed solution below, ...

So here's a more brute force method : minimizing the maximum distance
with some optimizations (heuristic for initial values and a not so naive
2D-space walker). There's probably some more optimizations for walking
the gradient but I'm happy with the current speed (~10s for 10000 points
on my 1.06GHz Core Duo).

This seems similar in spirit to Frank's solution, but less clean from
the Math PoV and more optimized for speed.

It outputs the solution and creates a gnuplot plotting file for you to
verify (you have the path used by the algorithm to find the center). I'm
no Gnuplot expert so unfortunately fancy colors are out and crappy
legends in...

Lionel

=============================================
include Math

# Uses 10 points by default
POINT_COUNT = (ARGV[0].to_i == 0) ? 10 : ARGV[0].to_i

class Point < Struct.new(:x, :y)
def self.random
Point.new(rand, rand)
end

def self.generate_samples(n)
(1..n).map { self.random }
end

def to_s
"(#{x}, #{y})"
end
end

class PotentialCenter < Struct.new(:x, :y, :points)
# Square distance with a point
def square_distance(point)
((point.x - x) ** 2) + ((point.y - y) ** 2)
end

# Maximum square distance with my points
# It's cached because it's called several times and can be very costly
def max_square_distance
@sq_dist ||= (points.map { |point| square_distance(point) }.max)
end

# Compute the best move with given distance (ro)
# and angles (thetas)
# returns nil if none is better than the current position
def best_move(ro, thetas)
potential_moves = thetas.map { |theta|
new_x = x + (ro * cos(theta))
new_y = y + (ro * sin(theta))
PotentialCenter.new(new_x,new_y,points)
}
choose_move_1(potential_moves)
end

# Dumber, faster
def choose_move_1(potential_moves)
potential_moves.detect { |move|
max_square_distance > move.max_square_distance
}
end

# Look for the shortest path, slower (~1.4x on average with many points)
def choose_move_2(potential_moves)
potential_moves.reject! { |move|
max_square_distance <= move.max_square_distance
}
potential_moves.min { |c1,c2|
c1.max_square_distance <=> c2.max_square_distance
}
end

# Check that the given offset is big enough to move
# the current position
# (floating point rounding errors prevent infinite precision...)
def can_be_moved_by(offset)
((x + offset) != x) || ((y + offset) != y)
end

def to_s
"(#{x}, #{y})"
end
end

def to_s
end
end

def chose_original_center_and_move_amount(points)
# Start with a good potential (around the middle of the point cloud)
max_x = points.map{|p| p.x}.max
min_x = points.map{|p| p.x}.min
x = (max_x + min_x) / 2
max_y = points.map{|p| p.y}.max
min_y = points.map{|p| p.y}.min
y = (max_y + min_y) / 2
center = PotentialCenter.new(x, y, points)
# The original displacement value is based on the fact that a truly
# random distribution gets a mean value with a precision of this order.
# The actual center is probably reachable with a few of these steps
move_amount = ((max_y - min_y) + (max_x - min_x)) / points.size
return [center, move_amount]
end

# Look for the best center by minimizing
# its maximum distance with the given points
# This is the same as minimizing its square
# (the maximum of the square distances)
# This can be seen as a dichotomy on a 2-dimensional space
def encircle(points)
center, ro = chose_original_center_and_move_amount(points)
# We'll move with polar coordinates (distance + angle : ro/theta)
# How many directions do we follow?
# each new angle creates a new (costly) max_square_distance call
# so we take the minimum needed to move around and climb the
angle_count = 3
thetas = (1..angle_count).map { |c| (2 * c * Math::PI) / angle_count}
iter = 0
# Trace our moves
center_list = []

loop do
center_list << center
iter += 1
puts "#{iter}: #{center.max_square_distance}"
# Is there a best move available ?
if new_center = center.best_move(ro, thetas)
center = new_center
else
ro /= 2
# We stop when ro becomes too small to move the center
break unless center.can_be_moved_by(ro)
end
end
puts "Circle found in #{iter} iterations"
return [ Circle.new(center, sqrt(center.max_square_distance)),
center_list ]
end

points = Point.generate_samples(POINT_COUNT)
t = Time.now
circle, center_list = encircle(points)
puts circle
puts "Found in: #{Time.now - t}s"

# Generate a GnuPlot command file to see the result
f = File.open('plotcommand.cmd', 'w')
f.print <<EOF
set parametric
set size square
set xrange [-0.2:1.2]
set yrange [-0.2:1.2]
set multiplot
plot "-"
#{points.map { |p| "#{p.x} #{p.y}" }.join("\n")}
end
plot "-" with lines
#{center_list.map { |p| "#{p.x} #{p.y}" }.join("\n")}
end
plot [0:2*pi]