Rosetta Code

Smallest enclosing circle problem: Difference between revisions

Page Discussion

Smallest enclosing circle problem (view source)

Revision as of 20:32, 7 April 2022

55,421 bytes added ,  2 years ago
m
→‎{{header|Phix}}: syntax coloured, made p2js compatible
m (→‎{{header|Phix}}: syntax coloured, made p2js compatible)
Line 554:
=={{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.
<!--<lang Phix>type point(sequencephixonline) 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 = (bX*cY - bY*cX)*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 = {(cY*B - bY*C)/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>
(bX*C - cX*B)/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 ==&gt; %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>
<!--</lang>-->
{{out}}
<pre>
Line 626 ⟶ 629:
{{trans|Python}}
Uses the last test from Julia, however, since that's given more accurately.
<!--<lang Phix>constant inf = 1e300*1e300,(phixonline)-->
<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
<!--</lang>-->
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>
Petelomax
7,794

edits

Retrieved from "https://rosettacode.org/wiki/Special:MobileDiff/303456"