Smallest enclosing circle problem: Difference between revisions

← Older edit

Smallest enclosing circle problem (view source)

Revision as of 17:02, 22 February 2024

62,208 bytes added , 2 months ago

Realize in F#

Nigel Galloway

2,171

edits

Revision as of 18:11, 12 September 2021 (view source) Shuisman (talk \| contribs) (→‎{{header\|Nim}}) ← Older edit		Latest revision as of 17:02, 22 February 2024 (view source) Nigel Galloway (talk \| contribs) (Realize in F#)
(7 intermediate revisions by 5 users not shown)
Line 12: =={{header\|11l}}== {{trans\|Nim}} {{trans\|Go}} <syntaxhighlight lang="11l">T Point = (Float, Float) T Circle Point c Float r F (c, r) .c = c .r = r F String() R ‘Center ’(.c)‘, Radius ’(.r) F getCircleCenter(Float bx, by, cx, cy) V b = bx * bx + by * by V c = cx * cx + cy * cy V d = bx * cy - by * cx R Point((cy * b - by * c) / (2 * d), (bx * c - cx * b) / (2 * d)) F contains(Circle ci, Point p) R sqlen(ci.c - p) <= ci.r ^ 2 F encloses(Circle ci, [Point] ps) L(p) ps I !contains(ci, p) R 0B R 1B F circleFrom3(Point a, b, c) V 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 R Circle(i, distance(i, a)) F circleFrom2(Point a, b) V c = Point((a.x + b.x) * 0.5, (a.y + b.y) * 0.5) R Circle(c, distance(a, b) / 2) F secTrivial([Point] rs) I rs.empty R Circle(Point(0, 0), 0) I rs.len == 1 R Circle(rs[0], 0) I rs.len == 2 R circleFrom2(rs[0], rs[1]) assert(rs.len == 3, ‘There shouldn't be more than three points.’) L(i) 2 L(j) i + 1 .< 3 V c = circleFrom2(rs[i], rs[j]) I encloses(c, rs) R c R circleFrom3(rs[0], rs[1], rs[2]) F welzlHelper([Point] &ps, [Point] rs, Int n) V rc = copy(rs) I n == 0 \| rc.len == 3 R secTrivial(rc) V idx = random:(n) V p = ps[idx] swap(&ps[idx], &ps[n - 1]) V d = welzlHelper(&ps, rc, n - 1) I contains(d, p) R d rc.append(p) R welzlHelper(&ps, rc, n - 1) F welzl([Point] ps) V pc = copy(ps) random:shuffle(&pc) R welzlHelper(&pc, [Point](), pc.len) V Tests = [[Point(0.0, 0.0), Point(0.0, 1.0), Point(1.0, 0.0)], [Point(5.0, -2.0), Point(-3.0, -2.0), Point(-2.0, 5.0), Point(1.0, 6.0), Point(0.0, 2.0)]] L(test) Tests print(welzl(test))</syntaxhighlight> {{out}} <pre> Center (0.5, 0.5), Radius 0.707106781 Center (1, 1), Radius 5 </pre> =={{header\|F_Sharp\|F#}}== This task uses [https://rosettacode.org/wiki/Centre_and_radius_of_a_circle_passing_through_3_points_in_a_plane#F# Centre and radius of a circle passing through 3 points in a plane (F#)] <syntaxhighlight lang="fsharp"> // Smallest enclosing circle problem. Nigel Galloway: February 22nd., 2024 let mcc points= let fN g=points\|>List.map(fun n->(n,d g n))\|>List.maxBy snd let P1,P2=Seq.allPairs points points\|>Seq.maxBy(fun(n,g)->d n g) let c1,r1=let ((nx,ny) as n),((gx,gy) as g)=P1,P2 in (((nx+gx)/2.0,(ny+gy)/2.0),(d n g)/2.0) match fN c1 with (_,d) when d=r1->(c1,r1) \|(P3,_)->let c2,r2=circ P1 P2 P3 match fN c2 with (_,d) when d=r2->(c2,r2) \|(P4,_)->[circ P1 P3 P4;circ P2 P3 P4]\|>List.minBy snd let testMCC n=let centre,radius=mcc n in printfn "Minimum Covering Circle is centred at %A with a radius of %f" centre radius testMCC [(0.0, 0.0);(0.0, 1.0);(1.0, 0.0)] testMCC [(5.0, -2.0);(-3.0, -2.0);(-2.0, 5.0);(1.0, 6.0);(0.0, 2.0)] testMCC [(15.85,6.05);(29.82,22.11);(4.84,20.32);(25.47,29.46);(33.65,17.31);(7.30,10.43);(28.15,6.67);(20.85,11.47);(22.83,2.07);(26.28,23.15);(14.39,30.24);(30.24,25.23)] testMCC [(15.85,6.05);(29.82,22.11);(4.84,20.32);(25.47,29.46);(33.65,17.31);(7.30,10.43);(28.15,6.67);(20.85,11.47);(22.83,2.08);(26.28,23.15);(14.39,30.24);(30.24,25.23)] </syntaxhighlight> {{out}} <pre> Minimum Covering Circle is centred at (0.5, 0.5) with a radius of 0.707107 Minimum Covering Circle is centred at (1.0, 1.0) with a radius of 5.000000 Minimum Covering Circle is centred at (18.97851566, 16.2654108) with a radius of 14.708624 Minimum Covering Circle is centred at (18.97952365, 16.27401209) with a radius of 14.707010 </pre> =={{header\|Go}}== {{trans\|Wren}} <~~lang~~syntaxhighlight lang="go">package main import ( Line 142 ⟶ 256: fmt.Println(welzl(test)) } }</~~lang~~syntaxhighlight> {{out}} Line 148 ⟶ 262: {(0.500000, 0.500000), 0.707107} {(1.000000, 1.000000), 5.000000} </pre> =={{header\|jq}}== {{works with\|jq}} '''Works with gojq, the Go implementation of jq''' '''Adapted from [[#Wren\|Wren]]''' In the following, a Point (x,y) is represented by a JSON array [x,y], and a Circle by a JSON object of the form {center, radius}. <syntaxhighlight lang="jq"> # Vector sum of the input array of Points def sum: transpose \| map(add); # The square of the distance between two points def distSq: def sq: . * .; . as [$a, $b] \| (($a[0] - $b[0]) \| sq) +(($a[1] - $b[1]) \| sq) ; # The distance between two points def distance: distSq \| sqrt; # Does the input Circle contain the point $p? def contains($p): ([.center, $p] \| distSq) <= .radius * .radius; # Does the input Circle contain the given array of points? def encloses($ps): . as $c \| all($ps[]; . as $p \| $c \| contains($p)); # The center of a circle defined by 3 points def getCircleCenter(bx; by; cx; cy): (bxbx + byby) as $b \| (cxcx + cycy) as $c \| (bxcy - bycx) as $d \| [(cy$b - by$c) / (2 * $d), (bx$c - cx$b) / (2 * $d)]; # The (smallest) circle that intersects 1, 2 or 3 distinct points def Circle: if length == 1 then {center: .[0], radius: 0} elif length == 2 then {center: (sum \| map(./2)), radius: (distance/2) } elif length == 3 then . as [$a, $b, $c] \| ([getCircleCenter($b[0] - $a[0]; $b[1] - $a[1]; $c[0] - $a[0]; $c[1] - $a[1]), $a] \| sum) as $i \| {center: $i, radius: ([$i, $a]\|distance)} else "Circle/0 expects an input array of from 1 to 3 Points" \| error end; # The smallest enclosing circle for n <= 3 def secTrivial: . as $rs \| length as $length \| if $length > 3 then "Internal error: secTrivial given more than 3 points." \| error elif $length == 0 then [[0, 0]] \| Circle elif $length <= 2 then Circle else first( range(0; 2) as $i \| range($i+1; 3) as $j \| ([$rs[$i], $rs[$j]] \| Circle) \| select(encloses($rs)) ) // Circle end; # Use the Welzl algorithm to find the SEC of the incoming array of points def welzl: # Helper function: # $rs is an array of points on the boundary; # $n keeps track of how many points remain: def welzl_($rs; $n): if $n == 0 or ($rs\|length) == 3 then $rs \| secTrivial else 0 as $idx # or Rand($n) \| .[$idx] as $p \| .[$idx] = .[$n-1] \| .[$n-1] = $p \| welzl_($rs; $n-1) as $d \| if ($d \| contains($p)) then $d else welzl_($rs + [$p]; $n-1) end end ; welzl_([]; length); def tests: [ [0, 0], [0, 1], [1, 0] ], [ [5, -2], [-3, -2], [-2, 5], [1, 6], [0, 2] ] ; tests \| welzl \| "Circle(\(.center), \(.radius))" </syntaxhighlight> {{output}} <pre> Circle([0.5,0.5], 0.7071067811865476) Circle([1,1], 5) </pre> =={{header\|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. <~~lang~~syntaxhighlight lang="julia">import Base.pop!, Base.push!, Base.length, Base.* using LinearAlgebra, Random Line 323 ⟶ 542: testwelzl() </~~lang~~syntaxhighlight>{{out}} <pre> For points: [[0.0, 0.0], [0.0, 1.0], [1.0, 0.0]] Line 343 ⟶ 562: =={{header\|Mathematica}}/{{header\|Wolfram Language}}== <~~lang~~syntaxhighlight ~~Mathematica~~lang="mathematica">f=BoundingRegion[#,"MinBall"]&; f[{{0,0},{0,1},{1,0}}] f[{{5,-2},{-3,-2},{-2,5},{1,6},{0,2}}] f[{{2,4,-1},{1,5,-3},{8,-4,1},{3,9,-5}}] f[{{-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`}}]</~~lang~~syntaxhighlight> {{out}} <pre>Ball[{0.5,0.5},0.707107] Line 357 ⟶ 576: {{trans\|Go}} {{trans\|Wren}} <~~lang~~syntaxhighlight ~~Nim~~lang="nim">import math, random, strutils type Line 458 ⟶ 677: for test in Tests: echo welzl(test)</~~lang~~syntaxhighlight> {{out}} Line 466 ⟶ 685: =={{header\|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. <!--<syntaxhighlight lang="phix">(phixonline)--> ~~<lang Phix>type point(sequence) return true end type~~ <span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span> ~~enum CENTRE, RADIUS~~ <span style="color: #008080;">type</span> <span style="color: #000000;">point</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">s</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">return</span> <span style="color: #004600;">true</span> <span style="color: #008080;">end</span> <span style="color: #008080;">type</span> ~~type circle(sequence) return true end type~~ <span style="color: #008080;">enum</span> <span style="color: #000000;">CENTRE</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">RADIUS</span> <span style="color: #008080;">type</span> <span style="color: #000000;">circle</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">s</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">return</span> <span style="color: #004600;">true</span> <span style="color: #008080;">end</span> <span style="color: #008080;">type</span> ~~function distance(point a, b)~~ ~~return sqrt(sum(sq_power(sq_sub(a,b),2)))~~ <span style="color: #008080;">function</span> <span style="color: #000000;">distance</span><span style="color: #0000FF;">(</span><span style="color: #000000;">point</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">b</span><span style="color: #0000FF;">)</span> ~~end function~~ <span style="color: #008080;">return</span> <span style="color: #7060A8;">sqrt</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">sum</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">sq_power</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">sq_sub</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">,</span><span style="color: #000000;">b</span><span style="color: #0000FF;">),</span><span style="color: #000000;">2</span><span style="color: #0000FF;">)))</span> <span style="color: #008080;">end</span> <span style="color: #008080;">function</span> ~~function in_circle(point p, circle c)~~ ~~return distance(p,c[CENTRE]) <= c[RADIUS]~~ <span style="color: #008080;">function</span> <span style="color: #000000;">in_circle</span><span style="color: #0000FF;">(</span><span style="color: #000000;">point</span> <span style="color: #000000;">p</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">circle</span> <span style="color: #000000;">c</span><span style="color: #0000FF;">)</span> ~~end function~~ <span style="color: #008080;">return</span> <span style="color: #000000;">distance</span><span style="color: #0000FF;">(</span><span style="color: #000000;">p</span><span style="color: #0000FF;">,</span><span style="color: #000000;">c</span><span style="color: #0000FF;">[</span><span style="color: #000000;">CENTRE</span><span style="color: #0000FF;">])</span> <span style="color: #0000FF;"><=</span> <span style="color: #000000;">c</span><span style="color: #0000FF;">[</span><span style="color: #000000;">RADIUS</span><span style="color: #0000FF;">]</span> <span style="color: #008080;">end</span> <span style="color: #008080;">function</span> ~~function circle_from2(point a, b)~~ ~~-- return the smallest circle that intersects 2 points:~~ <span style="color: #008080;">function</span> <span style="color: #000000;">circle_from2</span><span style="color: #0000FF;">(</span><span style="color: #000000;">point</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">b</span><span style="color: #0000FF;">)</span> ~~point midpoint = sq_div(sq_add(a,b),2)~~ <span style="color: #000080;font-style:italic;">-- return the smallest circle that intersects 2 points:</span> ~~atom halfdiameter = distance(a,b)/2~~ <span style="color: #000000;">point</span> <span style="color: #000000;">midpoint</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sq_div</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">sq_add</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">,</span><span style="color: #000000;">b</span><span style="color: #0000FF;">),</span><span style="color: #000000;">2</span><span style="color: #0000FF;">)</span> ~~circle res = { midpoint, halfdiameter }~~ <span style="color: #004080;">atom</span> <span style="color: #000000;">halfdiameter</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">distance</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">,</span><span style="color: #000000;">b</span><span style="color: #0000FF;">)/</span><span style="color: #000000;">2</span> ~~return res~~ <span style="color: #000000;">circle</span> <span style="color: #000000;">res</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{</span> <span style="color: #000000;">midpoint</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">halfdiameter</span> <span style="color: #0000FF;">}</span> ~~end function~~ <span style="color: #008080;">return</span> <span style="color: #000000;">res</span> <span style="color: #008080;">end</span> <span style="color: #008080;">function</span> ~~function circle_from3(point a, b, c)~~ ~~-- return a unique circle that intersects three points~~ <span style="color: #008080;">function</span> <span style="color: #000000;">circle_from3</span><span style="color: #0000FF;">(</span><span style="color: #000000;">point</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">b</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">c</span><span style="color: #0000FF;">)</span> ~~point bma = sq_sub(b,a),~~ <span style="color: #000080;font-style:italic;">-- return a unique circle that intersects three points </span> ~~cma = sq_sub(c,a)~~ <span style="color: #000000;">point</span> <span style="color: #000000;">bma</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sq_sub</span><span style="color: #0000FF;">(</span><span style="color: #000000;">b</span><span style="color: #0000FF;">,</span><span style="color: #000000;">a</span><span style="color: #0000FF;">),</span> ~~atom {{aX,aY},{bX,bY},{cX,cY}} = {a,bma,cma},~~ <span style="color: #000000;">cma</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sq_sub</span><span style="color: #0000FF;">(</span><span style="color: #000000;">c</span><span style="color: #0000FF;">,</span><span style="color: #000000;">a</span><span style="color: #0000FF;">)</span> ~~B = sum(sq_power(bma,2)),~~ <span style="color: #004080;">atom</span> <span style="color: #0000FF;">{{</span><span style="color: #000000;">aX</span><span style="color: #0000FF;">,</span><span style="color: #000000;">aY</span><span style="color: #0000FF;">},{</span><span style="color: #000000;">bX</span><span style="color: #0000FF;">,</span><span style="color: #000000;">bY</span><span style="color: #0000FF;">},{</span><span style="color: #000000;">cX</span><span style="color: #0000FF;">,</span><span style="color: #000000;">cY</span><span style="color: #0000FF;">}}</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">a</span><span style="color: #0000FF;">,</span><span style="color: #000000;">bma</span><span style="color: #0000FF;">,</span><span style="color: #000000;">cma</span><span style="color: #0000FF;">},</span> ~~C = sum(sq_power(cma,2)),~~ <span style="color: #000000;">B</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sum</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">sq_power</span><span style="color: #0000FF;">(</span><span style="color: #000000;">bma</span><span style="color: #0000FF;">,</span><span style="color: #000000;">2</span><span style="color: #0000FF;">)),</span> ~~D = (bXcY - bYcX)2~~ <span style="color: #000000;">C</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sum</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">sq_power</span><span style="color: #0000FF;">(</span><span style="color: #000000;">cma</span><span style="color: #0000FF;">,</span><span style="color: #000000;">2</span><span style="color: #0000FF;">)),</span> ~~point centre = {(cYB - bYC)/D + aX,~~ <span style="color: #000000;">D</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">bX</span><span style="color: #0000FF;"></span><span style="color: #000000;">cY</span> <span style="color: #0000FF;">-</span> <span style="color: #000000;">bY</span><span style="color: #0000FF;"></span><span style="color: #000000;">cX</span><span style="color: #0000FF;">)</span><span style="color: #000000;">2</span> ~~(bXC - cXB)/D + aY }~~ <span style="color: #000000;">point</span> <span style="color: #000000;">centre</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{(</span><span style="color: #000000;">cY</span><span style="color: #0000FF;"></span><span style="color: #000000;">B</span> <span style="color: #0000FF;">-</span> <span style="color: #000000;">bY</span><span style="color: #0000FF;"></span><span style="color: #000000;">C</span><span style="color: #0000FF;">)/</span><span style="color: #000000;">D</span> <span style="color: #0000FF;">+</span> <span style="color: #000000;">aX</span><span style="color: #0000FF;">,</span> ~~atom radius = distance(centre,a) -- (=== b,c)~~ <span style="color: #0000FF;">(</span><span style="color: #000000;">bX</span><span style="color: #0000FF;"></span><span style="color: #000000;">C</span> <span style="color: #0000FF;">-</span> <span style="color: #000000;">cX</span><span style="color: #0000FF;"></span><span style="color: #000000;">B</span><span style="color: #0000FF;">)/</span><span style="color: #000000;">D</span> <span style="color: #0000FF;">+</span> <span style="color: #000000;">aY</span> <span style="color: #0000FF;">}</span> ~~circle res = { centre, radius }~~ <span style="color: #004080;">atom</span> <span style="color: #000000;">radius</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">distance</span><span style="color: #0000FF;">(</span><span style="color: #000000;">centre</span><span style="color: #0000FF;">,</span><span style="color: #000000;">a</span><span style="color: #0000FF;">)</span> <span style="color: #000080;font-style:italic;">-- (=== b,c)</span> ~~return res~~ <span style="color: #000000;">circle</span> <span style="color: #000000;">res</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{</span> <span style="color: #000000;">centre</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">radius</span> <span style="color: #0000FF;">}</span> ~~end function~~ <span style="color: #008080;">return</span> <span style="color: #000000;">res</span> <span style="color: #008080;">end</span> <span style="color: #008080;">function</span> ~~function trivial(sequence r)~~ ~~integer l = length(r)~~ <span style="color: #008080;">function</span> <span style="color: #000000;">trivial</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">r</span><span style="color: #0000FF;">)</span> ~~switch l do~~ <span style="color: #004080;">integer</span> <span style="color: #000000;">l</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">r</span><span style="color: #0000FF;">)</span> ~~case 0: return {{0,0},0}~~ <span style="color: #008080;">switch</span> <span style="color: #000000;">l</span> <span style="color: #008080;">do</span> ~~case 1: return {r[1],0}~~ <span style="color: #008080;">case</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">:</span> <span style="color: #008080;">return</span> <span style="color: #0000FF;">{{</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span><span style="color: #000000;">0</span><span style="color: #0000FF;">},</span><span style="color: #000000;">0</span><span style="color: #0000FF;">}</span> ~~case 2: return circle_from2(r[1],r[2])~~ <span style="color: #008080;">case</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">:</span> <span style="color: #008080;">return</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">r</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">],</span><span style="color: #000000;">0</span><span style="color: #0000FF;">}</span> ~~case 3: return circle_from3(r[1],r[2],r[3])~~ <span style="color: #008080;">case</span> <span style="color: #000000;">2</span><span style="color: #0000FF;">:</span> <span style="color: #008080;">return</span> <span style="color: #000000;">circle_from2</span><span style="color: #0000FF;">(</span><span style="color: #000000;">r</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">],</span><span style="color: #000000;">r</span><span style="color: #0000FF;">[</span><span style="color: #000000;">2</span><span style="color: #0000FF;">])</span> ~~end switch~~ <span style="color: #008080;">case</span> <span style="color: #000000;">3</span><span style="color: #0000FF;">:</span> <span style="color: #008080;">return</span> <span style="color: #000000;">circle_from3</span><span style="color: #0000FF;">(</span><span style="color: #000000;">r</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">],</span><span style="color: #000000;">r</span><span style="color: #0000FF;">[</span><span style="color: #000000;">2</span><span style="color: #0000FF;">],</span><span style="color: #000000;">r</span><span style="color: #0000FF;">[</span><span style="color: #000000;">3</span><span style="color: #0000FF;">])</span> ~~end function~~ <span style="color: #008080;">end</span> <span style="color: #008080;">switch</span> <span style="color: #008080;">end</span> <span style="color: #008080;">function</span> ~~function welzlr(sequence p, r)~~ ~~if p={} or length(r)=3 then return trivial(r) end if~~ <span style="color: #008080;">function</span> <span style="color: #000000;">welzlr</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">p</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">r</span><span style="color: #0000FF;">)</span> ~~point p1 = p[1]~~ <span style="color: #008080;">if</span> <span style="color: #000000;">p</span><span style="color: #0000FF;">={}</span> <span style="color: #008080;">or</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">r</span><span style="color: #0000FF;">)=</span><span style="color: #000000;">3</span> <span style="color: #008080;">then</span> <span style="color: #008080;">return</span> <span style="color: #000000;">trivial</span><span style="color: #0000FF;">(</span><span style="color: #000000;">r</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span> ~~p = p[2..$]~~ <span style="color: #000000;">point</span> <span style="color: #000000;">p1</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">p</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">]</span> ~~circle d = welzlr(p, r)~~ <span style="color: #000000;">p</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">p</span><span style="color: #0000FF;">[</span><span style="color: #000000;">2</span><span style="color: #0000FF;">..$]</span> ~~if in_circle(p1,d) then return d end if~~ <span style="color: #000000;">circle</span> <span style="color: #000000;">d</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">welzlr</span><span style="color: #0000FF;">(</span><span style="color: #000000;">p</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">r</span><span style="color: #0000FF;">)</span> ~~return welzlr(p, append(r,p1))~~ <span style="color: #008080;">if</span> <span style="color: #000000;">in_circle</span><span style="color: #0000FF;">(</span><span style="color: #000000;">p1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">d</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">then</span> <span style="color: #008080;">return</span> <span style="color: #000000;">d</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span> ~~end function~~ <span style="color: #008080;">return</span> <span style="color: #000000;">welzlr</span><span style="color: #0000FF;">(</span><span style="color: #000000;">p</span><span style="color: #0000FF;">,</span> <span style="color: #7060A8;">append</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">deep_copy</span><span style="color: #0000FF;">(</span><span style="color: #000000;">r</span><span style="color: #0000FF;">),</span><span style="color: #000000;">p1</span><span style="color: #0000FF;">))</span> <span style="color: #008080;">end</span> <span style="color: #008080;">function</span> ~~procedure welzl(sequence p)~~ ~~circle c = welzlr(shuffle(p),{})~~ <span style="color: #008080;">procedure</span> <span style="color: #000000;">welzl</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">p</span><span style="color: #0000FF;">)</span> ~~string s = sprintf("centre %v radius %.14g",c)~~ <span style="color: #000000;">circle</span> <span style="color: #000000;">c</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">welzlr</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">shuffle</span><span style="color: #0000FF;">(</span><span style="color: #000000;">p</span><span style="color: #0000FF;">),{})</span> ~~printf(1,"Points %v ==> %s\n",{p,s})~~ <span style="color: #004080;">string</span> <span style="color: #000000;">s</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sprintf</span><span style="color: #0000FF;">(</span><span style="color: #008000;">"centre %v radius %.14g"</span><span style="color: #0000FF;">,</span><span style="color: #000000;">c</span><span style="color: #0000FF;">)</span> ~~end procedure~~ <span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"Points %v ==> %s\n"</span><span style="color: #0000FF;">,{</span><span style="color: #000000;">p</span><span style="color: #0000FF;">,</span><span style="color: #000000;">s</span><span style="color: #0000FF;">})</span> <span style="color: #008080;">end</span> <span style="color: #008080;">procedure</span> ~~constant tests = {{{0, 0},{ 0, 1},{ 1,0}},~~ ~~{{5,-2},{-3,-2},{-2,5},{1,6},{0,2}}}~~ <span style="color: #008080;">constant</span> <span style="color: #000000;">tests</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{{{</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">},{</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">},{</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">0</span><span style="color: #0000FF;">}},</span> ~~papply(tests,welzl)</lang>~~ <span style="color: #0000FF;">{{</span><span style="color: #000000;">5</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">2</span><span style="color: #0000FF;">},{-</span><span style="color: #000000;">3</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">2</span><span style="color: #0000FF;">},{-</span><span style="color: #000000;">2</span><span style="color: #0000FF;">,</span><span style="color: #000000;">5</span><span style="color: #0000FF;">},{</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">6</span><span style="color: #0000FF;">},{</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span><span style="color: #000000;">2</span><span style="color: #0000FF;">}}}</span> <span style="color: #7060A8;">papply</span><span style="color: #0000FF;">(</span><span style="color: #000000;">tests</span><span style="color: #0000FF;">,</span><span style="color: #000000;">welzl</span><span style="color: #0000FF;">)</span> <!--</syntaxhighlight>--> {{out}} <pre> Line 538 ⟶ 760: {{trans\|Python}} Uses the last test from Julia, however, since that's given more accurately. <!--<syntaxhighlight lang="phix">(phixonline)--> ~~<lang Phix>constant inf = 1e3001e300,~~ <span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span> ~~nan = -(inf/inf),~~ <span style="color: #008080;">constant</span> <span style="color: #000000;">inf</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">1e300</span><span style="color: #0000FF;"></span><span style="color: #000000;">1e300</span><span style="color: #0000FF;">,</span> ~~eps = 1e-8~~ <span style="color: #000000;">nan</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">-(</span><span style="color: #000000;">inf</span><span style="color: #0000FF;">/</span><span style="color: #000000;">inf</span><span style="color: #0000FF;">),</span> <span style="color: #000000;">eps</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">1e-8</span> ~~function sqnorm(sequence p)~~ ~~return sum(sq_power(p,2))~~ <span style="color: #008080;">function</span> <span style="color: #000000;">sqnorm</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">p</span><span style="color: #0000FF;">)</span> ~~end function~~ <span style="color: #008080;">return</span> <span style="color: #7060A8;">sum</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">sq_power</span><span style="color: #0000FF;">(</span><span style="color: #000000;">p</span><span style="color: #0000FF;">,</span><span style="color: #000000;">2</span><span style="color: #0000FF;">))</span> <span style="color: #008080;">end</span> <span style="color: #008080;">function</span> ~~function sqdist(sequence p1, p2)~~ ~~return sqnorm(sq_sub(p1,p2))~~ <span style="color: #008080;">function</span> <span style="color: #000000;">sqdist</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">p1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">p2</span><span style="color: #0000FF;">)</span> ~~end function~~ <span style="color: #008080;">return</span> <span style="color: #000000;">sqnorm</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">sq_sub</span><span style="color: #0000FF;">(</span><span style="color: #000000;">p1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">p2</span><span style="color: #0000FF;">))</span> <span style="color: #008080;">end</span> <span style="color: #008080;">function</span> ~~class ProjectorStack~~ -- <span style="color: #008080;">function</span> <span style="color: #000000;">projector_mult</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">vs</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">v</span><span style="color: #0000FF;">)</span> ~~-- Stack of points that are shifted / projected to put first one at origin.~~ <span style="color: #004080;">sequence</span> <span style="color: #000000;">s</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span><span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">v</span><span style="color: #0000FF;">))</span> -- <span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">vs</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span> ~~sequence vs = {}~~ <span style="color: #004080;">sequence</span> <span style="color: #000000;">vi</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">vs</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]</span> <span style="color: #000000;">s</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sq_add</span><span style="color: #0000FF;">(</span><span style="color: #000000;">s</span><span style="color: #0000FF;">,</span><span style="color: #7060A8;">sq_mul</span><span style="color: #0000FF;">(</span><span style="color: #000000;">vi</span><span style="color: #0000FF;">,</span><span style="color: #7060A8;">sum</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">sq_mul</span><span style="color: #0000FF;">(</span><span style="color: #000000;">vi</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">v</span><span style="color: #0000FF;">))))</span> <span style="color: #008080;">end</span> <span style="color: #008080;">for</span> <span style="color: #008080;">return</span> <span style="color: #000000;">s</span> <span style="color: #008080;">end</span> <span style="color: #008080;">function</span> <span style="color: #000080;font-style:italic;">-- Gaertner Boundary:</span> <span style="color: #004080;">sequence</span> <span style="color: #000000;">empty_center</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">centers</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">square_radii</span><span style="color: #0000FF;">,</span> <span style="color: #000080;font-style:italic;">-- Stack of points that are shifted / projected to put first one at origin:</span> <span style="color: #000000;">projector</span> <span style="color: #008080;">function</span> <span style="color: #000000;">push_if_stable</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">pt</span><span style="color: #0000FF;">)</span> ~~procedure push(sequence v)~~ <span style="color: #008080;">if</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">centers</span><span style="color: #0000FF;">)</span> <span style="color: #0000FF;">==</span> <span style="color: #000000;">0</span> <span style="color: #008080;">then</span> ~~vs = append(vs, v)~~ <span style="color: #000000;">square_radii</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">append</span><span style="color: #0000FF;">(</span><span style="color: #000000;">square_radii</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">)</span> ~~end procedure~~ <span style="color: #000000;">centers</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">pt</span><span style="color: #0000FF;">}</span> <span style="color: #008080;">return</span> <span style="color: #004600;">true</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span> <span style="color: #004080;">sequence</span> <span style="color: #000000;">q0</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">centers</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">],</span> <span style="color: #000000;">center</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">centers</span><span style="color: #0000FF;">[$],</span> <span style="color: #000000;">C</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sq_sub</span><span style="color: #0000FF;">(</span><span style="color: #000000;">center</span><span style="color: #0000FF;">,</span><span style="color: #000000;">q0</span><span style="color: #0000FF;">),</span> <span style="color: #000000;">Qm</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sq_sub</span><span style="color: #0000FF;">(</span><span style="color: #000000;">pt</span><span style="color: #0000FF;">,</span><span style="color: #000000;">q0</span><span style="color: #0000FF;">),</span> <span style="color: #000000;">Qm_bar</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">projector_mult</span><span style="color: #0000FF;">(</span><span style="color: #000000;">projector</span><span style="color: #0000FF;">,</span><span style="color: #000000;">Qm</span><span style="color: #0000FF;">),</span> <span style="color: #000000;">residue</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sq_sub</span><span style="color: #0000FF;">(</span><span style="color: #000000;">Qm</span><span style="color: #0000FF;">,</span><span style="color: #000000;">Qm_bar</span><span style="color: #0000FF;">)</span> <span style="color: #004080;">atom</span> <span style="color: #000000;">r2</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">square_radii</span><span style="color: #0000FF;">[$],</span> <span style="color: #000000;">e</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">sqdist</span><span style="color: #0000FF;">(</span><span style="color: #000000;">Qm</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">C</span><span style="color: #0000FF;">)</span> <span style="color: #0000FF;">-</span> <span style="color: #000000;">r2</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">z</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">2</span> <span style="color: #0000FF;"></span> <span style="color: #000000;">sqnorm</span><span style="color: #0000FF;">(</span><span style="color: #000000;">residue</span><span style="color: #0000FF;">),</span> <span style="color: #000000;">tol</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">eps</span> <span style="color: #0000FF;"></span> <span style="color: #7060A8;">max</span><span style="color: #0000FF;">(</span><span style="color: #000000;">r2</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">),</span> <span style="color: #000000;">isstable</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">abs</span><span style="color: #0000FF;">(</span><span style="color: #000000;">z</span><span style="color: #0000FF;">)</span> <span style="color: #0000FF;">></span> <span style="color: #000000;">tol</span> <span style="color: #008080;">if</span> <span style="color: #000000;">isstable</span> <span style="color: #008080;">then</span> <span style="color: #004080;">sequence</span> <span style="color: #000000;">center_new</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sq_add</span><span style="color: #0000FF;">(</span><span style="color: #000000;">center</span><span style="color: #0000FF;">,</span><span style="color: #7060A8;">sq_mul</span><span style="color: #0000FF;">(</span><span style="color: #000000;">e</span><span style="color: #0000FF;">/</span><span style="color: #000000;">z</span><span style="color: #0000FF;">,</span><span style="color: #000000;">residue</span><span style="color: #0000FF;">))</span> <span style="color: #004080;">atom</span> <span style="color: #000000;">r2new</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">r2</span> <span style="color: #0000FF;">+</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">e</span> <span style="color: #0000FF;"></span> <span style="color: #000000;">e</span><span style="color: #0000FF;">)</span> <span style="color: #0000FF;">/</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">2</span> <span style="color: #0000FF;"></span> <span style="color: #000000;">z</span><span style="color: #0000FF;">)</span> <span style="color: #000000;">projector</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">append</span><span style="color: #0000FF;">(</span><span style="color: #000000;">projector</span><span style="color: #0000FF;">,</span><span style="color: #7060A8;">sq_div</span><span style="color: #0000FF;">(</span><span style="color: #000000;">residue</span><span style="color: #0000FF;">,</span><span style="color: #7060A8;">sqrt</span><span style="color: #0000FF;">(</span><span style="color: #000000;">sqnorm</span><span style="color: #0000FF;">(</span><span style="color: #000000;">residue</span><span style="color: #0000FF;">))))</span> <span style="color: #000000;">centers</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">append</span><span style="color: #0000FF;">(</span><span style="color: #000000;">centers</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">center_new</span><span style="color: #0000FF;">)</span> <span style="color: #000000;">square_radii</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">append</span><span style="color: #0000FF;">(</span><span style="color: #000000;">square_radii</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">r2new</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span> <span style="color: #008080;">return</span> <span style="color: #000000;">isstable</span> <span style="color: #008080;">end</span> <span style="color: #008080;">function</span> <span style="color: #008080;">procedure</span> <span style="color: #000000;">pop_gb</span><span style="color: #0000FF;">()</span> <span style="color: #004080;">integer</span> <span style="color: #000000;">n</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">centers</span><span style="color: #0000FF;">)</span> <span style="color: #000000;">centers</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">centers</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">..$-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">]</span> <span style="color: #000000;">square_radii</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">square_radii</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">..$-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">]</span> <span style="color: #008080;">if</span> <span style="color: #000000;">n</span> <span style="color: #0000FF;">>=</span> <span style="color: #000000;">2</span> <span style="color: #008080;">then</span> <span style="color: #000000;">projector</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">projector</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">..$-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">]</span> <span style="color: #000080;font-style:italic;">-- (pop/discard)</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span> <span style="color: #008080;">end</span> <span style="color: #008080;">procedure</span> <span style="color: #008080;">function</span> <span style="color: #000000;">isinside</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">pt</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">center</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">atom</span> <span style="color: #000000;">sqradius</span><span style="color: #0000FF;">)</span> <span style="color: #004080;">atom</span> <span style="color: #000000;">r2</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">sqdist</span><span style="color: #0000FF;">(</span><span style="color: #000000;">pt</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">center</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">return</span> <span style="color: #000000;">r2</span><span style="color: #0000FF;">!=</span><span style="color: #000000;">nan</span> <span style="color: #008080;">and</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">r2</span><span style="color: #0000FF;"><=</span><span style="color: #000000;">sqradius</span> <span style="color: #008080;">or</span> <span style="color: #7060A8;">abs</span><span style="color: #0000FF;">(</span><span style="color: #000000;">r2</span><span style="color: #0000FF;">-</span><span style="color: #000000;">sqradius</span><span style="color: #0000FF;">)<</span><span style="color: #000000;">eps</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">end</span> <span style="color: #008080;">function</span> <span style="color: #008080;">function</span> <span style="color: #000000;">move_to_front</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">pts</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">integer</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">)</span> <span style="color: #004080;">sequence</span> <span style="color: #000000;">pt</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">pts</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]</span> <span style="color: #008080;">for</span> <span style="color: #000000;">j</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">i</span> <span style="color: #008080;">do</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">pts</span><span style="color: #0000FF;">[</span><span style="color: #000000;">j</span><span style="color: #0000FF;">],</span> <span style="color: #000000;">pt</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">pt</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">pts</span><span style="color: #0000FF;">[</span><span style="color: #000000;">j</span><span style="color: #0000FF;">]}</span> <span style="color: #008080;">end</span> <span style="color: #008080;">for</span> <span style="color: #008080;">return</span> <span style="color: #000000;">pts</span> <span style="color: #008080;">end</span> <span style="color: #008080;">function</span> <span style="color: #008080;">function</span> <span style="color: #000000;">ismaxlength</span><span style="color: #0000FF;">()</span> <span style="color: #008080;">return</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">centers</span><span style="color: #0000FF;">)</span> <span style="color: #0000FF;">==</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">empty_center</span><span style="color: #0000FF;">)</span> <span style="color: #0000FF;">+</span> <span style="color: #000000;">1</span> <span style="color: #008080;">end</span> <span style="color: #008080;">function</span> <span style="color: #008080;">function</span> <span style="color: #000000;">_welzl</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">pts</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">integer</span> <span style="color: #000000;">pos</span><span style="color: #0000FF;">)</span> <span style="color: #004080;">integer</span> <span style="color: #000000;">support_count</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">s</span> <span style="color: #004080;">sequence</span> <span style="color: #000000;">center</span> <span style="color: #0000FF;">=</span> <span style="color: #008080;">iff</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">centers</span><span style="color: #0000FF;">)?</span><span style="color: #000000;">centers</span><span style="color: #0000FF;">[$]:</span><span style="color: #000000;">empty_center</span><span style="color: #0000FF;">)</span> <span style="color: #004080;">atom</span> <span style="color: #000000;">sqradius</span> <span style="color: #0000FF;">=</span> <span style="color: #008080;">iff</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">centers</span><span style="color: #0000FF;">)?</span><span style="color: #000000;">square_radii</span><span style="color: #0000FF;">[$]:</span><span style="color: #000000;">0</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">if</span> <span style="color: #000000;">ismaxlength</span><span style="color: #0000FF;">()</span> <span style="color: #008080;">then</span> <span style="color: #008080;">return</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">center</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">sqradius</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">}</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span> <span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">pos</span> <span style="color: #008080;">do</span> <span style="color: #008080;">if</span> <span style="color: #008080;">not</span> <span style="color: #000000;">isinside</span><span style="color: #0000FF;">(</span><span style="color: #000000;">pts</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">],</span> <span style="color: #000000;">center</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">sqradius</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">then</span> <span style="color: #004080;">bool</span> <span style="color: #000000;">isstable</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">push_if_stable</span><span style="color: #0000FF;">(</span><span style="color: #000000;">pts</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">])</span> <span style="color: #008080;">if</span> <span style="color: #000000;">isstable</span> <span style="color: #008080;">then</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">center</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">sqradius</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">s</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">_welzl</span><span style="color: #0000FF;">(</span><span style="color: #000000;">pts</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">)</span> <span style="color: #000000;">pop_gb</span><span style="color: #0000FF;">()</span> <span style="color: #000000;">pts</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">move_to_front</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">deep_copy</span><span style="color: #0000FF;">(</span><span style="color: #000000;">pts</span><span style="color: #0000FF;">),</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">)</span> <span style="color: #000000;">support_count</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">s</span> <span style="color: #0000FF;">+</span> <span style="color: #000000;">1</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span> <span style="color: #008080;">end</span> <span style="color: #008080;">for</span> <span style="color: #008080;">return</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">center</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">sqradius</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">support_count</span><span style="color: #0000FF;">}</span> <span style="color: #008080;">end</span> <span style="color: #008080;">function</span> <span style="color: #008080;">function</span> <span style="color: #000000;">find_max_excess</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">center</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">pts</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">atom</span> <span style="color: #000000;">sqradius</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">integer</span> <span style="color: #000000;">k1</span><span style="color: #0000FF;">)</span> <span style="color: #004080;">atom</span> <span style="color: #000000;">err_max</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">inf</span> <span style="color: #004080;">integer</span> <span style="color: #000000;">k_max</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">k1</span> <span style="color: #0000FF;">-</span> <span style="color: #000000;">1</span> <span style="color: #008080;">for</span> <span style="color: #000000;">k</span><span style="color: #0000FF;">=</span><span style="color: #000000;">k_max</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">pts</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span> <span style="color: #004080;">sequence</span> <span style="color: #000000;">pt</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">pts</span><span style="color: #0000FF;">[</span><span style="color: #000000;">k</span><span style="color: #0000FF;">]</span> <span style="color: #004080;">atom</span> <span style="color: #000000;">err</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">sqdist</span><span style="color: #0000FF;">(</span><span style="color: #000000;">pt</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">center</span><span style="color: #0000FF;">)</span> <span style="color: #0000FF;">-</span> <span style="color: #000000;">sqradius</span> <span style="color: #008080;">if</span> <span style="color: #000000;">err</span> <span style="color: #0000FF;">></span> <span style="color: #000000;">err_max</span> <span style="color: #008080;">then</span> <span style="color: #000000;">err_max</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">err</span> <span style="color: #000000;">k_max</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">k</span> <span style="color: #0000FF;">+</span> <span style="color: #000000;">1</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span> <span style="color: #008080;">end</span> <span style="color: #008080;">for</span> <span style="color: #008080;">return</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">err_max</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">k_max</span> <span style="color: #0000FF;">-</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">}</span> <span style="color: #008080;">end</span> <span style="color: #008080;">function</span> <span style="color: #008080;">procedure</span> <span style="color: #000000;">welzl</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">points</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">integer</span> <span style="color: #000000;">maxiterations</span><span style="color: #0000FF;">=</span><span style="color: #000000;">2000</span><span style="color: #0000FF;">)</span> <span style="color: #004080;">sequence</span> <span style="color: #000000;">pts</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">points</span> <span style="color: #000000;">empty_center</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #000000;">nan</span><span style="color: #0000FF;">,</span><span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">pts</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">]))</span> <span style="color: #000000;">centers</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{}</span> <span style="color: #000000;">square_radii</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{}</span> <span style="color: #000000;">projector</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{}</span> <span style="color: #004080;">integer</span> <span style="color: #000000;">t</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">s_new</span> ~~function pop()~~ <span style="color: #0000FF;">{</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">center</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">atom</span> <span style="color: #000000;">sqradius</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">integer</span> <span style="color: #000000;">s</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">_welzl</span><span style="color: #0000FF;">(</span><span style="color: #000000;">pts</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">t</span><span style="color: #0000FF;">)</span> ~~sequence v = vs[$]~~ <span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">maxiterations</span> <span style="color: #008080;">do</span> ~~vs = vs[1..$-1]~~ <span style="color: #004080;">atom</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">e</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">k</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">find_max_excess</span><span style="color: #0000FF;">(</span><span style="color: #000000;">center</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">pts</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">sqradius</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">t</span> <span style="color: #0000FF;">+</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">)</span> ~~return v~~ <span style="color: #008080;">if</span> <span style="color: #000000;">e</span> <span style="color: #0000FF;"><=</span> <span style="color: #000000;">eps</span> <span style="color: #008080;">then</span> <span style="color: #008080;">exit</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span> ~~end function~~ <span style="color: #004080;">sequence</span> <span style="color: #000000;">pt</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">pts</span><span style="color: #0000FF;">[</span><span style="color: #000000;">k</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">{}</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">push_if_stable</span><span style="color: #0000FF;">(</span><span style="color: #000000;">pt</span><span style="color: #0000FF;">)</span> ~~function mult(sequence v)~~ <span style="color: #0000FF;">{</span><span style="color: #000000;">center</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">sqradius</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">s_new</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">_welzl</span><span style="color: #0000FF;">(</span><span style="color: #000000;">pts</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">s</span><span style="color: #0000FF;">)</span> ~~sequence s = repeat(0,length(v))~~ <span style="color: #000000;">pop_gb</span><span style="color: #0000FF;">()</span> ~~for i=1 to length(vs) do~~ <span style="color: #000000;">pts</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">move_to_front</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">deep_copy</span><span style="color: #0000FF;">(</span><span style="color: #000000;">pts</span><span style="color: #0000FF;">),</span> <span style="color: #000000;">k</span><span style="color: #0000FF;">)</span> ~~sequence vi = vs[i]~~ <span style="color: #000000;">t</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">s</span> <span style="color: #0000FF;">+</span> <span style="color: #000000;">1</span> ~~s = sq_add(s,sq_mul(vi,sum(sq_mul(vi, v))))~~ <span style="color: #000000;">s</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">s_new</span> <span style="color: #0000FF;">+</span> <span style="color: #000000;">1</span> ~~end for~~ <span style="color: #008080;">end</span> <span style="color: #008080;">for</span> ~~return s~~ <span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"For points: "</span><span style="color: #0000FF;">)</span> <span style="color: #7060A8;">pp</span><span style="color: #0000FF;">(</span><span style="color: #000000;">points</span><span style="color: #0000FF;">,{</span><span style="color: #004600;">pp_Indent</span><span style="color: #0000FF;">,</span><span style="color: #000000;">12</span><span style="color: #0000FF;">,</span><span style="color: #004600;">pp_FltFmt</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"%11.8f"</span><span style="color: #0000FF;">})</span> ~~end function~~ <span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">" Center is at: %v\n"</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">center</span><span style="color: #0000FF;">})</span> <span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">" Radius is: %.16g\n"</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">{</span><span style="color: #7060A8;">sqrt</span><span style="color: #0000FF;">(</span><span style="color: #000000;">sqradius</span><span style="color: #0000FF;">)})</span> ~~end class~~ <span style="color: #008080;">end</span> <span style="color: #008080;">procedure</span> ~~class GaertnerBoundary~~ <span style="color: #008080;">constant</span> <span style="color: #000000;">tests</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{{{</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">},</span> <span style="color: #0000FF;">{</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">},</span> <span style="color: #0000FF;">{</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">}},</span> ~~public:~~ <span style="color: #0000FF;">{{</span><span style="color: #000000;">5</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">2</span><span style="color: #0000FF;">},</span> <span style="color: #0000FF;">{-</span><span style="color: #000000;">3</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">2</span><span style="color: #0000FF;">},</span> <span style="color: #0000FF;">{-</span><span style="color: #000000;">2</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">5</span><span style="color: #0000FF;">},</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">6</span><span style="color: #0000FF;">},</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">2</span><span style="color: #0000FF;">}},</span> ~~sequence empty_center,~~ <span style="color: #0000FF;">{{</span><span style="color: #000000;">2</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">4</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">},</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">5</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">3</span><span style="color: #0000FF;">},</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">8</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">4</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">},</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">3</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">9</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">5</span><span style="color: #0000FF;">}},</span> ~~centers = {},~~ <span style="color: #0000FF;">{{-</span><span style="color: #000000;">0.6400900432782123</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">2.643703255134232</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0.4016122094706093</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.8959601399652273</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">1.1624046608380516</span><span style="color: #0000FF;">},</span> ~~square_radii = {}~~ <span style="color: #0000FF;">{</span> <span style="color: #000000;">0.5632393652621324</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">0.015981105190064373</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">2.193725940351997</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">0.9094586577358262</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0.7165036664470906</span><span style="color: #0000FF;">},</span> ~~ProjectorStack projector = new()~~ <span style="color: #0000FF;">{-</span><span style="color: #000000;">1.7704367632976061</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0.2518283698686299</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">1.3810444289625348</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">0.597516704360172</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.089645656962647</span><span style="color: #0000FF;">},</span> ~~end class~~ <span style="color: #0000FF;">{</span> <span style="color: #000000;">1.3448578652803103</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">0.18140877132223002</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">0.4288734015080915</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.53271973321691</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0.41896461833399573</span><span style="color: #0000FF;">},</span> <span style="color: #0000FF;">{</span> <span style="color: #000000;">0.2132336397213029</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0.07659442168765788</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0.1486364315310990</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">2.3386893481333795</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">2.3000459709300927</span><span style="color: #0000FF;">},</span> ~~function push_if_stable(GaertnerBoundary bound, sequence pt)~~ <span style="color: #0000FF;">{</span> <span style="color: #000000;">0.6023153188617328</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0.3051735340025314</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.0732600437151525</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.1148388039984898</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0.047605838564167786</span><span style="color: #0000FF;">},</span> ~~if length(bound.centers) == 0 then~~ <span style="color: #0000FF;">{</span> <span style="color: #000000;">1.3655523661720959</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0.5461416420929995</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">0.09321951900362889</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">1.004771137760985</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1.6532914656050117</span><span style="color: #0000FF;">},</span> ~~bound.square_radii = append(bound.square_radii, 0)~~ <span style="color: #0000FF;">{-</span><span style="color: #000000;">0.14974382165751837</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">0.5375672589202939</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">0.15845596754003466</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">0.2750720340454088</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">0.441247015836271</span><span style="color: #0000FF;">}}}</span> ~~bound.centers = {pt}~~ <span style="color: #7060A8;">papply</span><span style="color: #0000FF;">(</span><span style="color: #000000;">tests</span><span style="color: #0000FF;">,</span><span style="color: #000000;">welzl</span><span style="color: #0000FF;">)</span> ~~return true~~ <!--</syntaxhighlight>--> ~~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 isstable~~ ~~end 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 sqradius~~ ~~end 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)<eps~~ ~~end 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 pts~~ ~~end function~~ ~~function ismaxlength(GaertnerBoundary bound)~~ ~~return length(bound.centers) == length(bound.empty_center) + 1~~ ~~end 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 res~~ ~~end 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)</lang>~~ {{out}} <pre> Line 753 ⟶ 941: =={{header\|Python}}== {{trans\|Julia}} <~~lang~~syntaxhighlight lang="python">import numpy as np class ProjectorStack: Line 912 ⟶ 1,100: print(" Center is at: ", nsphere.center) print(" Radius is: ", np.sqrt(nsphere.sqradius), "\n") </~~lang~~syntaxhighlight>{{out}} <pre> For points: [[0. 0.] Line 949 ⟶ 1,137: =={{header\|Raku}}== {{trans\|Go}} <syntaxhighlight lang="raku" ~~perl6~~line>class P { has ($.x, $.y) is rw; # Point method gist { "({$.x~", "~$.y})" } } Line 1,015 ⟶ 1,203: } say "Solution for smallest circle enclosing {$_.gist} :\n", welzl $_ for @tests;</~~lang~~syntaxhighlight> {{out}} <pre>Solution for smallest circle enclosing ((0, 0) (0, 1) (1, 0)) : Line 1,027 ⟶ 1,215: It is based on Welzl's algorithm and follows closely the C++ code [https://www.geeksforgeeks.org/minimum-enclosing-circle-set-2-welzls-algorithm/?ref=rp here]. <~~lang~~syntaxhighlight ~~ecmascript~~lang="wren">import "random" for Random var Rand = Random.new() Line 1,143 ⟶ 1,331: ] for (test in tests) System.print(Circle.welzl(test))</~~lang~~syntaxhighlight> {{out}}