I'm working on modernizing Rosetta Code's infrastructure. Starting with communications. Please accept this time-limited open invite to RC's Slack.. --Michael Mol (talk) 20:59, 30 May 2020 (UTC)

# Smallest enclosing circle problem

Smallest enclosing circle problem is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.

The smallest enclosing circle problem (aka minimum covering circle problem, bounding circle problem) is a mathematical problem of computing the smallest circle that contains all of a given set of points in the Euclidean plane.

Initially it was proposed by the English mathematician James Joseph Sylvester in 1857.

Find circle of smallest radius containing all of given points.

• Circle is defined by it's center and radius;
• Points are defined by their coordinates in n-dimensional space;
• Circle (sphere) contains point when distance between point and circle center <= circle radius.

## Go

Translation of: Wren
`package main import (    "fmt"    "log"    "math"    "math/rand"    "time") type Point struct{ x, y float64 } type Circle struct {    c Point    r float64} // string representation of a Pointfunc (p Point) String() string { return fmt.Sprintf("(%f, %f)", p.x, p.y) } // returns the square of the distance between two pointsfunc distSq(a, b Point) float64 {    return (a.x-b.x)*(a.x-b.x) + (a.y-b.y)*(a.y-b.y)} // returns the center of a circle defined by 3 pointsfunc getCircleCenter(bx, by, cx, cy float64) Point {    b := bx*bx + by*by    c := cx*cx + cy*cy    d := bx*cy - by*cx    return Point{(cy*b - by*c) / (2 * d), (bx*c - cx*b) / (2 * d)}} // returns whether a circle contains the point 'p'func (ci Circle) contains(p Point) bool {    return distSq(ci.c, p) <= ci.r*ci.r} // returns whether a circle contains a slice of pointsfunc (ci Circle) encloses(ps []Point) bool {    for _, p := range ps {        if !ci.contains(p) {            return false        }    }    return true} // string representation of a Circlefunc (ci Circle) String() string { return fmt.Sprintf("{%v, %f}", ci.c, ci.r) } // returns a unique circle that intersects 3 pointsfunc circleFrom3(a, b, c Point) Circle {    i := getCircleCenter(b.x-a.x, b.y-a.y, c.x-a.x, c.y-a.y)    i.x += a.x    i.y += a.y    return Circle{i, math.Sqrt(distSq(i, a))}} // returns smallest circle that intersects 2 pointsfunc circleFrom2(a, b Point) Circle {    c := Point{(a.x + b.x) / 2, (a.y + b.y) / 2}    return Circle{c, math.Sqrt(distSq(a, b)) / 2}} // returns smallest enclosing circle for n <= 3func secTrivial(rs []Point) Circle {    size := len(rs)    if size > 3 {        log.Fatal("There shouldn't be more than 3 points.")    }    if size == 0 {        return Circle{Point{0, 0}, 0}    }    if size == 1 {        return Circle{rs[0], 0}    }    if size == 2 {        return circleFrom2(rs[0], rs[1])    }    for i := 0; i < 2; i++ {        for j := i + 1; j < 3; j++ {            c := circleFrom2(rs[i], rs[j])            if c.encloses(rs) {                return c            }        }    }    return circleFrom3(rs[0], rs[1], rs[2])} // helper function for Welzl methodfunc welzlHelper(ps, rs []Point, n int) Circle {    rc := make([]Point, len(rs))    copy(rc, rs) // 'rs' passed by value so need to copy    if n == 0 || len(rc) == 3 {        return secTrivial(rc)    }    idx := rand.Intn(n)    p := ps[idx]    ps[idx], ps[n-1] = ps[n-1], p    d := welzlHelper(ps, rc, n-1)    if d.contains(p) {        return d    }    rc = append(rc, p)    return welzlHelper(ps, rc, n-1)} // applies the Welzl algorithm to find the SECfunc welzl(ps []Point) Circle {    var pc = make([]Point, len(ps))    copy(pc, ps)    rand.Shuffle(len(pc), func(i, j int) {        pc[i], pc[j] = pc[j], pc[i]    })    return welzlHelper(pc, []Point{}, len(pc))} func main() {    rand.Seed(time.Now().UnixNano())    tests := [][]Point{        {Point{0, 0}, Point{0, 1}, Point{1, 0}},        {Point{5, -2}, Point{-3, -2}, Point{-2, 5}, Point{1, 6}, Point{0, 2}},    }    for _, test := range tests {        fmt.Println(welzl(test))    }}`
Output:
```{(0.500000, 0.500000), 0.707107}
{(1.000000, 1.000000), 5.000000}
```

## Julia

N-dimensional solution using the Welzl algorithm. Derived from the BoundingSphere.jl module at https://github.com/JuliaFEM. See also Bernd Gärtner's paper at people.inf.ethz.ch/gaertner/subdir/texts/own_work/esa99_final.pdf.

`import Base.pop!, Base.push!, Base.length, Base.*using LinearAlgebra, Random struct ProjectorStack{P <: AbstractVector}    vs::Vector{P}endBase.push!(p::ProjectorStack, v) = (push!(p.vs, v); p)Base.pop!(p::ProjectorStack) = (pop!(p.vs); p)Base.:*(p::ProjectorStack, v) = (s = zero(v); for vi in p.vs s += vi * dot(vi, v) end; s) """    GärtnerBoundary See the passage regarding M_B in Section 4 of Gärtner's paper."""mutable struct GärtnerBoundary{P<:AbstractVector, F<:AbstractFloat}    centers::Vector{P}    square_radii::Vector{F}    # projection onto of affine space spanned by points    # shifted such that first point becomes origin    projector::ProjectorStack{P}    empty_center::P # center of nsphere spanned by empty boundaryend function GärtnerBoundary(pts)    P = eltype(pts)    F = eltype(P)    projector, centers, square_radii = ProjectorStack(P[]), P[], F[]    empty_center = F(NaN) * first(pts)    GärtnerBoundary(centers, square_radii, projector, empty_center)end function push_if_stable!(b::GärtnerBoundary, pt)    if isempty(b)        push!(b.square_radii, zero(eltype(pt)))        push!(b.centers, pt)        dim = length(pt)        return true    end    q0, center = first(b.centers), b.centers[end]    C, r2  = center - q0, b.square_radii[end]    Qm, M = pt - q0, b.projector    Qm_bar = M*Qm    residue, e = Qm - Qm_bar, sqdist(Qm, C) - r2    z, tol = 2*sqnorm(residue), eps(eltype(pt)) * max(r2, one(r2))    isstable = abs(z) > tol    if isstable        center_new  = center + (e/z) * residue        r2new = r2 + (e^2)/(2z)        push!(b.projector, residue / norm(residue))        push!(b.centers, center_new)        push!(b.square_radii, r2new)    end    return isstableend function Base.pop!(b::GärtnerBoundary)    n = length(b)    pop!(b.centers)    pop!(b.square_radii)    if n >= 2        pop!(b.projector)    end    return bend struct NSphere{P,F}    center::P    sqradius::Fend function isinside(pt, nsphere::NSphere; atol=0, rtol=0)    r2, R2 = sqdist(pt, center(nsphere)), sqradius(nsphere)    return r2 <= R2 || isapprox(r2, R2;atol=atol^2,rtol=rtol^2)endallinside(pts, nsphere; kw...) = all(pt -> isinside(pt, nsphere; kw...), pts) function move_to_front!(pts, i)    pt = pts[i]    for j in eachindex(pts)        pts[j], pt = pt, pts[j]        j == i && break    end    ptsend prefix(pts, i) = view(pts, 1:i)Base.length(b::GärtnerBoundary) = length(b.centers)Base.isempty(b::GärtnerBoundary) = length(b) == 0center(b::NSphere) = b.centerradius(b::NSphere) = sqrt(b.sqradius)sqradius(b::NSphere) = b.sqradiusdist(p1,p2) = norm(p1-p2)sqdist(p1::AbstractVector, p2::AbstractVector) = sqnorm(p1-p2)sqnorm(p) = sum(abs2,p)ismaxlength(b::GärtnerBoundary) = length(b) == length(b.empty_center) + 1 function NSphere(b::GärtnerBoundary)    return isempty(b) ? NSphere(b.empty_center, 0.0) :        NSphere(b.centers[end], b.square_radii[end])end function welzl!(pts, bdry)    bdry_len, support_count, nsphere = length(bdry), 0, NSphere(bdry)    ismaxlength(bdry) && return nsphere, 0    for i in eachindex(pts)        pt = pts[i]        if !isinside(pt, nsphere)            pts_i = prefix(pts, i-1)            isstable = push_if_stable!(bdry, pt)            if isstable                nsphere, s = welzl!(pts_i, bdry)                pop!(bdry)                move_to_front!(pts, i)                support_count = s + 1            end        end    end    return nsphere, support_countend function find_max_excess(nsphere, pts, k1)    T = eltype(first(pts))    e_max, k_max = T(-Inf), k1 -1    for k in k1:length(pts)        pt = pts[k]        e = sqdist(pt, center(nsphere)) - sqradius(nsphere)        if  e > e_max            e_max = e            k_max = k        end    end    return e_max, k_maxend function welzl(points, maxiterations=2000)    pts = deepcopy(points)    bdry, t = GärtnerBoundary(pts), 1    nsphere, s = welzl!(prefix(pts, t), bdry)    for i in 1:maxiterations        e, k = find_max_excess(nsphere, pts, t + 1)        P = eltype(pts)        F = eltype(P)        e <= eps(F) && break        pt = pts[k]        push_if_stable!(bdry, pt)        nsphere_new, s_new = welzl!(prefix(pts, s), bdry)        pop!(bdry)        move_to_front!(pts, k)        nsphere = nsphere_new        t, s = s + 1, s_new + 1    end    return nsphereend function testwelzl()    testdata =[        [[0.0, 0.0], [0.0, 1.0], [1.0, 0.0]],        [[5.0, -2.0], [-3.0, -2.0], [-2.0, 5.0], [1.0, 6.0], [0.0, 2.0]],        [[2.0, 4.0, -1.0], [1.0, 5.0, -3.0], [8.0, -4.0, 1.0], [3.0, 9.0, -5.0]],        [randn(5) for _ in 1:8]    ]    for test in testdata        nsphere = welzl(test)        println("For points: ", test)        println("    Center is at: ", nsphere.center)        println("    Radius is: ", radius(nsphere), "\n")    endend testwelzl() `
Output:
```For points: [[0.0, 0.0], [0.0, 1.0], [1.0, 0.0]]
Center is at: [0.5, 0.5]

For points: [[5.0, -2.0], [-3.0, -2.0], [-2.0, 5.0], [1.0, 6.0], [0.0, 2.0]]
Center is at: [1.0, 1.0]

For points: [[2.0, 4.0, -1.0], [1.0, 5.0, -3.0], [8.0, -4.0, 1.0], [3.0, 9.0, -5.0]]
Center is at: [5.5, 2.5, -2.0]

For points: [[-0.6400900432782123, 2.643703255134232, 0.4016122094706093, 1.8959601399652273, -1.1624046608380516], [0.5632393652621324, -0.015981105190064373, -2.193725940351997, -0.9094586577358262, 0.7165036664470906], [-1.7704367632976061, 0.2518283698686299, -1.3810444289625348, -0.597516704360172, 1.089645656962647], [1.3448578652803103, -0.18140877132223002, -0.4288734015080915, 1.53271973321691, 0.41896461833399573], [0.2132336397213029, 0.07659442168765788, 0.148636431531099, 2.3386893481333795, -2.3000459709300927], [0.6023153188617328, 0.3051735340025314, 1.0732600437151525, 1.1148388039984898, 0.047605838564167786], [1.3655523661720959, 0.5461416420929995, -0.09321951900362889, -1.004771137760985, 1.6532914656050117], [-0.14974382165751837, -0.5375672589202939, -0.15845596754003466, -0.2750720340454088, -0.441247015836271]]
Center is at: [0.0405866077655439, 0.5556683897981481, -0.2313678300507679, 0.6074066023194586, -0.2003463948612026]
```

## Phix

Based on the same code as Wren, and likewise limited to 2D circles - I believe (but cannot prove) the main barrier to coping with more than two dimensions lies wholly within the circle_from3() routine.

`type point(sequence) return true end typeenum CENTRE, RADIUStype circle(sequence) return true end type function distance(point a, b)    return sqrt(sum(sq_power(sq_sub(a,b),2)))end function function in_circle(point p, circle c)    return distance(p,c[CENTRE]) <= c[RADIUS]end function function circle_from2(point a, b)    -- return the smallest circle that intersects 2 points:    point midpoint = sq_div(sq_add(a,b),2)    atom halfdiameter = distance(a,b)/2    circle res = { midpoint, halfdiameter }    return resend function function circle_from3(point a, b, c)    -- return a unique circle that intersects three points     point bma = sq_sub(b,a),          cma = sq_sub(c,a)    atom {{aX,aY},{bX,bY},{cX,cY}} = {a,bma,cma},         B = sum(sq_power(bma,2)),         C = sum(sq_power(cma,2)),         D = (bX*cY - bY*cX)*2     point centre = {(cY*B - bY*C)/D + aX,                     (bX*C - cX*B)/D + aY }    atom radius = distance(centre,a) -- (=== b,c)    circle res = { centre, radius }    return resend function function trivial(sequence r)    integer l = length(r)    switch l do        case 0: return {{0,0},0}        case 1: return {r[1],0}        case 2: return circle_from2(r[1],r[2])        case 3: return circle_from3(r[1],r[2],r[3])    end switchend function function welzlr(sequence p, r)    if p={} or length(r)=3 then return trivial(r) end if    point p1 = p[1]    p = p[2..\$]    circle d = welzlr(p, r)    if in_circle(p1,d) then return d end if    return welzlr(p, append(r,p1))end function procedure welzl(sequence p)    circle c = welzlr(shuffle(p),{})    string s = sprintf("centre %v radius %.14g",c)    printf(1,"Points %v ==> %s\n",{p,s})end procedure constant tests = {{{0, 0},{ 0, 1},{ 1,0}},                  {{5,-2},{-3,-2},{-2,5},{1,6},{0,2}}}papply(tests,welzl)`
Output:
```Points {{0,0},{0,1},{1,0}} ==> centre {0.5,0.5} radius 0.70710678118655
Points {{5,-2},{-3,-2},{-2,5},{1,6},{0,2}} ==> centre {1,1} radius 5
```

### n-dimensional

Translation of: Python

Uses the last test from Julia, however, since that's given more accurately.

`constant inf = 1e300*1e300,         nan = -(inf/inf),         eps = 1e-8 function sqnorm(sequence p)    return sum(sq_power(p,2))end function function sqdist(sequence p1, p2)    return sqnorm(sq_sub(p1,p2))end function class ProjectorStack    --    -- Stack of points that are shifted / projected to put first one at origin.    --    sequence vs = {}     procedure push(sequence v)        vs = append(vs, v)    end procedure     function pop()        sequence v = vs[\$]        vs = vs[1..\$-1]        return v    end function     function mult(sequence v)        sequence s = repeat(0,length(v))        for i=1 to length(vs) do            sequence vi = vs[i]            s = sq_add(s,sq_mul(vi,sum(sq_mul(vi, v))))        end for        return s    end function end class class GaertnerBoundary  public:    sequence empty_center,             centers = {},             square_radii = {}    ProjectorStack projector = new()end class function push_if_stable(GaertnerBoundary bound, sequence pt)    if length(bound.centers) == 0 then        bound.square_radii = append(bound.square_radii, 0)        bound.centers = {pt}        return true    end if    ProjectorStack M = bound.projector    sequence q0 = bound.centers[1],             center = bound.centers[\$],             C = sq_sub(center,q0),             Qm = sq_sub(pt,q0),             Qm_bar = M.mult(Qm),             residue = sq_sub(Qm,Qm_bar)    atom r2 = bound.square_radii[\$],         e = sqdist(Qm, C) - r2,         z = 2 * sqnorm(residue),         tol = eps * max(r2, 1),         isstable = abs(z) > tol    if isstable then        sequence center_new = sq_add(center,sq_mul(e/z,residue))        atom r2new = r2 + (e * e) / (2 * z)        ProjectorStack p = bound.projector        p.push(sq_div(residue,sqrt(sqnorm(residue))))        bound.centers = append(bound.centers, center_new)        bound.square_radii = append(bound.square_radii, r2new)    end if    return isstableend function function pop(GaertnerBoundary bound)    integer n = length(bound.centers)    bound.centers = bound.centers[1..\$-1]    bound.square_radii = bound.square_radii[1..\$-1]    if n >= 2 then        ProjectorStack p = bound.projector        {} = p.pop()    end if    return bound end function class NSphere  public:    sequence center    atom sqradiusend class function isinside(sequence pt, NSphere nsphere)    atom r2 = sqdist(pt, nsphere.center),         R2 = nsphere.sqradius    if r2=nan then return false end if    return r2<=R2 or abs(r2-R2)<epsend function function move_to_front(sequence pts, integer i)    sequence pt = pts[i]    for j=1 to i do        {pts[j], pt} = {pt, pts[j]}    end for    return ptsend function function ismaxlength(GaertnerBoundary bound)    return length(bound.centers) == length(bound.empty_center) + 1end function function makeNSphere(GaertnerBoundary bound)    NSphere res    if length(bound.centers) == 0 then        res = new({bound.empty_center, 0.0})    else        res = new({bound.centers[\$], bound.square_radii[\$]})    end if    return resend function function _welzl(sequence pts, integer pos, GaertnerBoundary bdry)    integer support_count = 0    NSphere nsphere = makeNSphere(bdry)    if ismaxlength(bdry) then        return {nsphere, 0}    end if    for i=1 to pos do        if not isinside(pts[i], nsphere) then            bool isstable = push_if_stable(bdry, pts[i])            if isstable then                {nsphere, integer s} = _welzl(pts, i, bdry)                bdry = pop(bdry)                pts = move_to_front(pts, i)                support_count = s + 1            end if        end if    end for    return {nsphere, support_count}end function function find_max_excess(NSphere nsphere, sequence pts, integer k1)    atom err_max  = -inf    integer k_max = k1 - 1    for k=k_max to length(pts) do        sequence pt = pts[k]        atom err = sqdist(pt, nsphere.center) - nsphere.sqradius        if err > err_max then            err_max = err            k_max = k + 1        end if    end for    return {err_max, k_max - 1}end function procedure welzl(sequence points, integer maxiterations=2000)    sequence pts = points    GaertnerBoundary bdry = new({repeat(nan,length(pts[1]))})    integer t = 1    {NSphere nsphere, integer s} = _welzl(pts, t, bdry)    for i=1 to maxiterations do        atom {e, k} = find_max_excess(nsphere, pts, t + 1)        if e <= eps then exit end if        sequence pt = pts[k]        {} = push_if_stable(bdry, pt)        {NSphere nsphere_new, integer s_new} = _welzl(pts, s, bdry)        bdry = pop(bdry)        pts = move_to_front(pts, k)        nsphere = nsphere_new        t = s + 1        s = s_new + 1    end for    printf(1,"For points: ") pp(points,{pp_Indent,12,pp_FltFmt,"%11.8f"})    printf(1,"    Center is at: %v\n", {nsphere.center})    printf(1,"    Radius is: %.16g\n", {sqrt(nsphere.sqradius)})end procedure constant tests = {{{0, 0}, { 0, 1}, { 1, 0}},                  {{5,-2}, {-3,-2}, {-2, 5}, {1, 6}, {0, 2}},                  {{2, 4, -1}, {1, 5, -3}, {8, -4, 1}, {3, 9, -5}},                  {{-0.6400900432782123, 2.643703255134232,    0.4016122094706093, 1.8959601399652273,-1.1624046608380516},                   { 0.5632393652621324,-0.015981105190064373,-2.193725940351997, -0.9094586577358262, 0.7165036664470906},                   {-1.7704367632976061, 0.2518283698686299,  -1.3810444289625348,-0.597516704360172,  1.089645656962647},                   { 1.3448578652803103,-0.18140877132223002, -0.4288734015080915, 1.53271973321691,   0.41896461833399573},                   { 0.2132336397213029, 0.07659442168765788,  0.1486364315310990, 2.3386893481333795,-2.3000459709300927},                   { 0.6023153188617328, 0.3051735340025314,   1.0732600437151525, 1.1148388039984898, 0.047605838564167786},                   { 1.3655523661720959, 0.5461416420929995, -0.09321951900362889,-1.004771137760985,  1.6532914656050117},                   {-0.14974382165751837,-0.5375672589202939,-0.15845596754003466,-0.2750720340454088,-0.441247015836271}}}papply(tests,welzl)`
Output:
```For points: {{0,0}, {0,1}, {1,0}}
Center is at: {0.5,0.5}
For points: {{5,-2}, {-3,-2}, {-2,5}, {1,6}, {0,2}}
Center is at: {1,1}
For points: {{2,4,-1}, {1,5,-3}, {8,-4,1}, {3,9,-5}}
Center is at: {5.5,2.5,-2}
For points: {{-0.64009004, 2.64370326, 0.40161221, 1.89596014,-1.16240466},
{ 0.56323937,-0.01598111,-2.19372594,-0.90945866, 0.71650367},
{-1.77043676, 0.25182837,-1.38104443,-0.59751670, 1.08964566},
{ 1.34485787,-0.18140877,-0.42887340, 1.53271973, 0.41896462},
{ 0.21323364, 0.07659442, 0.14863643, 2.33868935,-2.30004597},
{ 0.60231532, 0.30517353, 1.07326004, 1.11483880, 0.04760584},
{ 1.36555237, 0.54614164,-0.09321952,-1.00477114, 1.65329147},
{-0.14974382,-0.53756726,-0.15845597,-0.27507203,-0.44124702}}
Center is at: {0.0405866078,0.5556683898,-0.2313678301,0.6074066023,-0.2003463949}
```

## Python

Translation of: Julia
`import numpy as np class ProjectorStack:    """    Stack of points that are shifted / projected to put first one at origin.    """    def __init__(self, vec):        self.vs = np.array(vec)     def push(self, v):        if len(self.vs) == 0:            self.vs = np.array([v])        else:            self.vs = np.append(self.vs, [v], axis=0)        return self     def pop(self):        if len(self.vs) > 0:            ret, self.vs = self.vs[-1], self.vs[:-1]            return ret     def __mul__(self, v):        s = np.zeros(len(v))        for vi in self.vs:            s = s + vi * np.dot(vi, v)        return s class GaertnerBoundary:    """        GärtnerBoundary     See the passage regarding M_B in Section 4 of Gärtner's paper.    """    def __init__(self, pts):        self.projector = ProjectorStack([])        self.centers, self.square_radii = np.array([]), np.array([])        self.empty_center = np.array([np.NaN for _ in pts[0]])  def push_if_stable(bound, pt):    if len(bound.centers) == 0:        bound.square_radii = np.append(bound.square_radii, 0.0)        bound.centers = np.array([pt])        return True    q0, center = bound.centers[0], bound.centers[-1]    C, r2  = center - q0, bound.square_radii[-1]    Qm, M = pt - q0, bound.projector    Qm_bar = M * Qm    residue, e = Qm - Qm_bar, sqdist(Qm, C) - r2    z, tol = 2 * sqnorm(residue), np.finfo(float).eps * max(r2, 1.0)    isstable = np.abs(z) > tol    if isstable:        center_new  = center + (e / z) * residue        r2new = r2 + (e * e) / (2 * z)        bound.projector.push(residue / np.linalg.norm(residue))        bound.centers = np.append(bound.centers, np.array([center_new]), axis=0)        bound.square_radii = np.append(bound.square_radii, r2new)    return isstable def pop(bound):    n = len(bound.centers)    bound.centers = bound.centers[:-1]    bound.square_radii = bound.square_radii[:-1]    if n >= 2:        bound.projector.pop()    return bound  class NSphere:    def __init__(self, c, sqr):        self.center = np.array(c)        self.sqradius = sqr def isinside(pt, nsphere, atol=1e-6, rtol=0.0):    r2, R2 = sqdist(pt, nsphere.center), nsphere.sqradius    return r2 <= R2 or np.isclose(r2, R2, atol=atol**2,rtol=rtol**2) def allinside(pts, nsphere, atol=1e-6, rtol=0.0):    for p in pts:        if not isinside(p, nsphere, atol, rtol):            return False    return True def move_to_front(pts, i):    pt = pts[i]    for j in range(len(pts)):        pts[j], pt = pt, np.array(pts[j])        if j == i:            break    return pts def dist(p1, p2):    return np.linalg.norm(p1 - p2) def sqdist(p1, p2):    return sqnorm(p1 - p2) def sqnorm(p):    return np.sum(np.array([x * x for x in p])) def ismaxlength(bound):    len(bound.centers) == len(bound.empty_center) + 1 def makeNSphere(bound):    if len(bound.centers) == 0:         return NSphere(bound.empty_center, 0.0)    return NSphere(bound.centers[-1], bound.square_radii[-1]) def _welzl(pts, pos, bdry):    support_count, nsphere = 0, makeNSphere(bdry)    if ismaxlength(bdry):        return nsphere, 0    for i in range(pos):        if not isinside(pts[i], nsphere):            isstable = push_if_stable(bdry, pts[i])            if isstable:                nsphere, s = _welzl(pts, i, bdry)                pop(bdry)                move_to_front(pts, i)                support_count = s + 1    return nsphere, support_count def find_max_excess(nsphere, pts, k1):    err_max, k_max = -np.Inf, k1 - 1    for (k, pt) in enumerate(pts[k_max:]):        err = sqdist(pt, nsphere.center) - nsphere.sqradius        if  err > err_max:            err_max, k_max = err, k + k1    return err_max, k_max - 1 def welzl(points, maxiterations=2000):    pts, eps = np.array(points, copy=True), np.finfo(float).eps    bdry, t = GaertnerBoundary(pts), 1    nsphere, s = _welzl(pts, t, bdry)    for i in range(maxiterations):        e, k = find_max_excess(nsphere, pts, t + 1)        if e <= eps:            break        pt = pts[k]        push_if_stable(bdry, pt)        nsphere_new, s_new = _welzl(pts, s, bdry)        pop(bdry)        move_to_front(pts, k)        nsphere = nsphere_new        t, s = s + 1, s_new + 1    return nsphere if __name__ == '__main__':    TESTDATA =[        np.array([[0.0, 0.0], [0.0, 1.0], [1.0, 0.0]]),        np.array([[5.0, -2.0], [-3.0, -2.0], [-2.0, 5.0], [1.0, 6.0], [0.0, 2.0]]),        np.array([[2.0, 4.0, -1.0], [1.0, 5.0, -3.0], [8.0, -4.0, 1.0], [3.0, 9.0, -5.0]]),        np.random.normal(size=(8, 5))    ]    for test in TESTDATA:        nsphere = welzl(test)        print("For points: ", test)        print("    Center is at: ", nsphere.center)        print("    Radius is: ", np.sqrt(nsphere.sqradius), "\n") `
Output:
```For points:  [[0. 0.]
[0. 1.]
[1. 0.]]
Center is at:  [0.5 0.5]

For points:  [[ 5. -2.]
[-3. -2.]
[-2.  5.]
[ 1.  6.]
[ 0.  2.]]
Center is at:  [1. 1.]

For points:  [[ 2.  4. -1.]
[ 1.  5. -3.]
[ 8. -4.  1.]
[ 3.  9. -5.]]
Center is at:  [ 5.5  2.5 -2. ]

For points:  [[-0.28550471 -0.12787614 -0.65666165 -0.57319826 -0.09253894]
[ 0.08947276  0.64778564 -0.54953015  0.05332253 -0.93055218]
[ 1.26721546 -1.26410984 -0.79937865 -0.45096792 -1.23946668]
[-0.62653779  0.56561466  0.62403237  0.23226903 -0.95820264]
[ 0.90973949 -0.9821474   0.41400032 -1.11268937  0.19568717]
[-0.4931165   0.94778404 -0.10534124  0.96431358  0.36504087]
[ 1.43895269 -0.69957774 -0.31486014 -1.14980913  0.42550193]
[ 1.4722404   1.19711551 -0.18735466  1.69208268  1.04594225]]
Center is at:  [ 1.00038333  0.08534583 -0.22414667  0.50927606 -0.15936026]
```

## Raku

Translation of: Go
`# 20201208 Raku programming solution class P { has (\$.x, \$.y) is rw;                # Point   method gist { "({\$.x~", "~\$.y})" }} class C { has (P \$.c, Numeric \$.r);            # Circle    method gist { "Center = " ~\$.c.gist ~ " and Radius = \$.r\n" }    # returns whether a circle contains the point 'p'   method contains(P \p --> Bool) { return distSq(\$.c, p) ≤ \$.r² }    method encloses(@ps --> Bool) {      @ps.map: { return False unless \$.contains(\$_) }      return True   } # returns whether a circle contains a slice of point} # returns the square of the distance between two pointssub distSq (P \a, P \b) { return (a.x - b.x)² + (a.y - b.y)² } sub getCenter (\bx, \by, \cx, \cy --> P) {   my (\b,\c,\d) =  bx²+by² , cx²+cy², bx*cy - by*cx;   P.new: x => (cy*b - by*c) / (2 * d), y => (bx*c - cx*b) / (2 * d)} # returns the center of a circle defined by 3 points sub circleFrom3 (P \a, P \b, P \c --> C)  {   my \k = \$ = getCenter(b.x - a.x, b.y - a.y, c.x - a.x, c.y - a.y);   k.x += a.x;   k.y += a.y;   return C.new: c => k, r => distSq(k, a).sqrt} # returns a unique circle that intersects 3 points sub circleFrom2 (P \a, P \b --> C ) {   my \center = P.new: x => ((a.x + b.x) / 2), y => ((a.y + b.y) / 2) ;   return C.new: c => center, r => (distSq(a, b).sqrt / 2)} # returns smallest circle that intersects 2 points sub secTrivial( @rs --> C ) {   given +@rs {      when *  > 3 { die "There shouldn't be more than 3 points." }      when * == 0 { return C.new: c => (P.new: x => 0, y => 0) , r => 0 }      when * == 1 { return C.new: c => @rs[0], r => 0 }      when * == 2 { return circleFrom2 @rs[0], @rs[1] }   }   for 0, 1 X 1, 2 -> ( \i, \j ) {      given circleFrom2(@rs[i], @rs[j]) { return \$_ if .encloses(@rs) }   }   return circleFrom3 @rs[0], @rs[1], @rs[2]} # returns smallest enclosing circle for n ≤ 3 sub welzlHelper( @ps is copy, @rs is copy , \n --> C ) {   return secTrivial(@rs) if ( n == 0 || +@rs == 3 );   my \p = @ps.shift;   given  welzlHelper(@ps, @rs, n-1) { return \$_ if .contains(p) }   return welzlHelper(@ps, @rs.append(p), n-1)} # helper function for Welzl method # applies the Welzl algorithm to find the SECsub welzl(@ps --> C) { return welzlHelper(@ps.pick(*), [], +@ps) } my @tests = (   [ (0,0), (0,1), (1,0) ],   [ (5,-2), (-3,-2), (-2,5), (1,6), (0,2) ]).map: {   @_.map: { P.new: x => \$_[0], y => \$_[1] }} say "Solution for smallest circle enclosing {\$_.gist} :\n", welzl \$_ for @tests;`
Output:
```Solution for smallest circle enclosing ((0, 0) (0, 1) (1, 0)) :
Center = (0.5, 0.5) and Radius = 0.7071067811865476

Solution for smallest circle enclosing ((5, -2) (-3, -2) (-2, 5) (1, 6) (0, 2)) :
Center = (1, 1) and Radius = 5```

## Wren

Well a circle is a two dimensional figure and so, despite any contradictory indications in the task description, that's what this solution provides.

It is based on Welzl's algorithm and follows closely the C++ code here.

`import "random" for Random var Rand = Random.new() class Point {    // returns the square of the distance between two points    static distSq(a, b) { (a.x - b.x)*(a.x - b.x) + (a.y - b.y)*(a.y - b.y) }     // returns the center of a circle defined by 3 points    static getCircleCenter(bx, by, cx, cy) {        var b = bx*bx + by*by        var c = cx*cx + cy*cy        var d = bx*cy - by*cx        return Point.new((cy*b - by*c) / (2 * d), (bx*c - cx*b) / (2 * d))    }     // constructs a new Point object    construct new(x, y) {        _x = x        _y = y    }     // basic properties    x { _x }    x=(o) { _x = o }    y { _y }    y=(o) { _y = o }     // returns a copy of the current instance    copy() { Point.new(_x, _y) }     // string representation    toString { "(%(_x), %(_y))" }} class Circle {    // returns a unique circle that intersects 3 points    static from(a, b, c) {        var i = Point.getCircleCenter(b.x - a.x, b.y - a.y, c.x - a.x, c.y - a.y)        i.x = i.x + a.x        i.y = i.y + a.y        return Circle.new(i, Point.distSq(i, a).sqrt)    }     // returns smallest circle that intersects 2 points    static from(a, b) {        var c = Point.new((a.x + b.x) / 2, (a.y + b.y) / 2)        return Circle.new(c, Point.distSq(a, b).sqrt/2)    }     // returns smallest enclosing circle for n <= 3    static secTrivial(rs) {        var size = rs.count        if (size > 3) Fiber.abort("There shouldn't be more than 3 points.")        if (size == 0) return Circle.new(Point.new(0, 0), 0)        if (size == 1) return Circle.new(rs[0], 0)        if (size == 2) return Circle.from(rs[0], rs[1])        for (i in 0..1) {            for (j in i+1..2) {                var c = Circle.from(rs[i], rs[j])                if (c.encloses(rs)) return c            }        }        return Circle.from(rs[0], rs[1], rs[2])    }     // private helper method for Welzl method    static welzl_(ps, rs, n) {        rs = rs.toList // passed by value so need to copy        if (n == 0 || rs.count == 3) return secTrivial(rs)        var idx = Rand.int(n)        var p = ps[idx]        ps[idx] = ps[n-1]        ps[n-1] = p        var d = welzl_(ps, rs, n-1)        if (d.contains(p)) return d        rs.add(p)        return welzl_(ps, rs, n-1)    }     // applies the Welzl algorithm to find the SEC    static welzl(ps) {        var pc = List.filled(ps.count, null)        for (i in 0...pc.count) pc[i] = ps[i].copy()        Rand.shuffle(pc)        return welzl_(pc, [], pc.count)    }     // constructs a new Circle object from its center point and its radius    construct new(c, r) {        _c = c.copy()        _r = r    }     // basic properties    c { _c }    r { _r }     // returns whether the current instance contains the point 'p'    contains(p) { Point.distSq(_c, p) <= _r * _r }     // returns whether the current instance contains a list of points    encloses(ps) {        for (p in ps) if (!contains(p)) return false        return true    }     // string representation    toString { "Center %(_c), Radius %(_r)" }} var tests = [    [Point.new(0, 0), Point.new(0, 1), Point.new(1, 0)],    [Point.new(5, -2), Point.new(-3, -2), Point.new(-2, 5), Point.new(1, 6), Point.new(0, 2)]] for (test in tests) System.print(Circle.welzl(test))`
Output:
```Center (0.5, 0.5), Radius 0.70710678118655