15232 Reputation

24 Badges

12 years, 80 days

MaplePrimes Activity

These are Posts that have been published by Kitonum

Thу post is the answer to

This message contains two procedures . The first procedure called  LargestIncircle  symbolically looks for the largest circle inscribed in a convex polygon . Polygon must be specified as a list of lists of its vertices . If the polygon does not have parallel sides , the problem has a unique solution , otherwise may be several solutions . It is easy to prove that there is always a maximum circle that touches three sides of the polygon . Procedure code based on this property of the optimal solution.

The second procedure called  MinCoveringCircle  symbolically seeking smallest circle that encloses a given set of points. The set of points can be anyone not necessarily being vertices of a convex polygon. Procedure code based on the following geometric properties of the optimal solution , which can easily be proved :

1) The minimum covering circle is unique.

2) The minimum covering circle of a set Set can be determined by at most three points in  Set  which lie on the boundary of the circle. If it is determined by only two points, then the line segment joining those two points must be a diameter of the minimum circle. If it is determined by three points, then the triangle consisting of those three points is not obtuse.


Code of the 1st procedure:



local n, x, y, P, S, T, Eqs, IsInside, OC, Pt, E, R, t, Sol, P0, R1, Circle;

uses combinat, geometry;


n:=nops(L); _EnvHorizontalName := x; _EnvVerticalName := y;

P:=[seq(point(A||i, L[i]), i=1..n)];

if nops(convexhull(P))<n then error "L should be a convex polygon" fi;

S:={seq(line(l||i,[A||i, A||(i+1)]), i=1..n-1), line(l||n,[A||n, A||1])};

Eqs:=map(lhs,{seq(Equation(l||i), i=1..n)});

T:=choose({seq(l||i, i=1..n)}, 3);


IsInside:=proc(Point, OC, Eqs)

convert([seq(is(eval(Eqs[i], {x=Point[1], y=Point[2]})*eval(Eqs[i], {x=OC[1], y=OC[2]})>0), i=1..n)], `and`);

end proc;


OC:=[add(L[i,1], i=1..n)/n, add(L[i,2], i=1..n)/n]; 

R:=0; point(Pt,[x,y]);

   for t in T do

     solve([distance(Pt,t[1])=distance(Pt,t[2]),      distance(Pt,t[2])=distance(Pt,t[3])],[x,y]);

     Sol:=map(a->[rhs(a[1]), rhs(a[2])], %);

     P0:=op(select(b->IsInside(b, OC, Eqs), Sol)); if P0<>NULL then R1:=distance(point(E,P0),      t[1]);

     if convert([seq(is(R1<=distance(E, i)), i=S minus t)], `and`) and is(R1>R) then R:=R1;      Circle:=[P0, R] fi; fi;



end proc:


Code the 2nd procedure:

MinCoveringCircle:=proc(Set::{set, list}(list))

local S, T, t, d, A, B, C, P, R, R1, O1, Coord, L;

uses combinat, geometry; 

if type(Set, list) then S:=convert(Set, set) else S:=Set fi;

T:=choose(S, 2);

   for t in T do

      d:=simplify(distance(point(A, t[1]), point(B, t[2]))); midpoint(C, A, B);

      if convert([seq(is(distance(point(P, p), C)<=d/2), p in S minus {t[1], t[2]})], `and`)   then return            [coordinates(C), d/2] fi;


T:=choose(S, 3);


   for t in T do

      point(A, t[1]); point(B, t[2]); point(C, t[3]);

      if is(distance(A, B)^2<=distance(B, C)^2+distance(A, C)^2 and        distance(A,C)^2<=distance(B,        C)^2+distance(A, B)^2 and distance(B, C)^2 <= distance(A, C)^2 + distance(A, B)^2) then

      circle(c, [A,B,C],'centername'=O1); R1:=radius(c); Coord:=coordinates(O1);

      if convert([seq(is(distance(point(P, p), O1)<=R1), p in S minus {t[1], t[2], t[3]})], `and`) and is(R1<R) then    R:=R1; L:=[Coord, R] fi; fi;



end proc: 


Examples of using the first procedure:



A:=plottools[disk](op(C), color=green):

B:=plottools[polygon](L, style=line, thickness=2):

plots[display](A,B, scaling=constrained);



The procedure can be used to find the largest circle inscribed in a convex shape, bounded by the curves, if these curves are approximated by broken lines:

L:=[seq([k,k^2], k=-1..2, 0.1)]:


A:=plottools[disk](op(C), color=green):

B:=plottools[polygon](L, style=line, thickness=2):

plots[display](A,B, scaling=constrained);



Examples of using the second procedure:



A:=plottools[circle](op(C), color=red, thickness=2):

B:=plottools[polygon](L, color=cyan):




The smallest circle covering 30 random points:

> roll:=rand(1..10):

L:=[seq([roll(), roll()], i=1..30)]:


A:=plottools[circle](op(C), color=red, thickness=2):

B:=seq(plottools[disk](L[i],0.07,color=blue), i=1..nops(L)):

C1:=plottools[disk](C[1], 0.07, color=red):

T:=plots[textplot]([C[1,1]-0.2, C[1,2]-0.2, "C"], font=[TIMES,ROMAN,16]):

plots[display](A, B, C1, T);



The work consists of two independent procedures. The first procedure  IsConvex  checks the convexity of a polygon. The second procedure  IsSimple  verifies the simplicity of a polygon. Formal argument   is the list of vertices of the polygon.

Regarding the basic concepts, see



local n, Z, f, i, x, y;




   for i to n do

    if  convert([seq(is(subs(x=j[1],y=j[2],f[i])<=0), j in {op(X)} minus  {X[i],X[irem(i,n)+1]})],`or`) and      convert([seq(is(subs(x=j[1],y=j[2],f[i])>=0),

    j in {op(X)} minus {X[i],X[irem(i,n)+1]})], `or`) then break fi;


if i<=n then return false else true fi;

end proc:



local n, Z, i, j, f, T, Q, x, y;

Z:=[op(X),X[1],X[2]]; n:=nops(X);

if n>nops({op(X)}) then   return false  fi;

   for i from 2 to nops(Z)-1 do

     if is((Z[i-1,1]-Z[i,1])*(Z[i+1,1]-Z[i,1])+(Z[i-1,2]-Z[i,2])*(Z[i+1,2]-Z[i,2]) =  sqrt((Z[i-1,1] -Z[i,1])^2+(Z[i-1,2]-Z[i,2])^2)*sqrt((Z[i+1,1]-Z[i,1])^2 +(Z[i+1,2]-Z[i,2])^2)) then return false fi;



_Envsignum0:= 0: 

   for i from 1 to n do

   T[i]:=[]; Q[i]:=[];

      for j from 1 to n do

      if modp(j-i,n)<>0 and modp(j-i,n)<>1 and modp(j-i,n)<>n-1 and                  not(signum(subs(x=Z[j,1],y=Z[j,2],f[i])*subs(x=Z[j+1,1],y=Z[j+1,2],f[i]))=-1 and             signum(subs(x=Z[i,1],y=Z[i,2],f[j])*subs(x=Z[i+1,1],y=Z[i+1,2],f[j]))=-1) then

      if (subs(x=Z[j,1],y=Z[j,2],f[i])=0 implies (signum((Z[j,1]-Z[i,1])*(Z[i+1,1]-Z[j,1]))=-1 or             signum((Z[j,2]-Z[i,2])*(Z[i+1,2]-Z[j,2]))=-1)) then

      T[i]:=[op(T[i]),1]; Q[i]:=[op(Q[i]),1] else  T[i]:=[op(T[i]),1]  fi; fi;   od;


convert([seq(nops(T[i])=n-3,i=1..n), seq(nops(Q[i])=n-3,i=1..n)],`and`)  

end proc:



X:=[[0,0],[1,0],[2,1],[3,0],[4,0],[2,2]]: IsConvex(X), IsSimple(X);

X:=[[0,0],[2,0],[1,1],[1,-1]]: IsConvex(X), IsSimple(X);

X:=[[0,0],[2,0],[1,1],[1,0], [-1,-1]]: IsConvex(X), IsSimple(X);

X:=[[0,0],[1,0],[1,2],[-2,2],[-2,-2],[3,-2],[3,4],[-4,4],[-4,-4],[5,-4],[5,6],[-6,6],[-6,-6],[7,-6],[7,8],[-6,8],[-6,7],[6,7],[6,-5],[-5,-5],[-5,5],[4,5],[4,-3],[-3,-3],[-3,3],[2,3],[2,-1],[-1,-1],[-1,1],[0,1]]: IsConvex(X), IsSimple(X);

X:=[seq([cos(2*Pi*k/17), sin(2*Pi*k/17)], k=0..16)]: IsConvex(X), IsSimple(X);


Edited: The variables  x  and  y  are made local.


Well-known problem is the problem of placing eight shess queens on an 8×8 chessboard so that no two queens attack each other. In this post, we consider the same problem of placing  m  shess queens on an  n×n  chessboard. The problem has a solution if  n>3  and  m<=n .

Work consists of two procedures. The first procedure  Queens  returns the total number of solutions and saves a complete list of all solutions (global variable  S ). The second procedure  QueensPic  shows the user-defined solutions from the list  S  on the board. Formal argument  t  is the number of solutions in each row of the display. The second procedure should be used in the standard interface, rather than in the classic one, since in the latter it may not work properly.

Queens := proc (m::posint, n::posint)

local It, K, l, L, M, P;

global S, p, q;

It := proc (L)

local P, k, i, j;

M := []; k := nops(L[1]);

for i in L do

for j to n do

if convert([seq(j <> i[s, 2], s = 1 .. k)], `and`) and convert([seq(l[k+1]-i[s, 1] <> i[s, 2]-j, s = 1 .. k)], `and`) and convert([seq(l[k+1]-i[s, 1] <> j-i[s, 2], s = 1 .. k)], `and`) then M := [op(M), [op(i), [l[k+1], j]]]


od; od;


end proc;

K := combinat:-choose([`$`(1 .. n)], m);

S := [];

for l in K do P := [];

L := [seq([[l[1], i]], i = 1 .. n)];

P := [op(P), op((It@@(m-1))(L))];

S := [op(S), op(P)]


p := args[1]; q := args[2];


end proc:


QueensPic := proc (M, t::posint)

local m, n, HL, VL, T, A, N;

uses plottools, plots;

m := p; n := q; N := nops(args[1]);

HL := seq(line([.5, .5+k], [.5+n, .5+k], color = black, thickness = 2), k = 0 .. n);

VL := seq(line([.5+k, .5], [.5+k, .5+n], color = black, thickness = 2), k = 0 .. n);

T := [seq(textplot([seq([op(M[i, j]), Q], j = 1 .. m)], color = red, font = [TIMES, ROMAN, 24]), i = 1 .. N)];

if m <= n and 3 < n then

A := seq(display(HL, VL, T[k], axes = none, scaling = constrained), k = 1 .. N), seq(display(plot([[0, 0]]), axes = none, scaling = constrained), k = 1 .. t*ceil(N/t)-N);

Matrix(ceil(N/t), t, [A]);



end proc:


Examples of work:

Queens(5, 6);  

S[70], S[140], S[210];

QueensPic([%], 3); 


[[1, 5], [2, 3], [3, 6], [4, 4], [6, 1]], [[1, 3], [2, 5], [4, 1], [5, 4], [6, 2]], [[2, 1], [3, 4], [4, 2], [5, 5], [6, 3]]


Two solutions of classic problem:

Queens(8, 8); 


QueensPic(%, 2);


[[[1, 5], [2, 8], [3, 4], [4, 1], [5, 7], [6, 2], [7, 6], [8, 3]], [[1, 6], [2, 1], [3, 5], [4, 2], [5, 8], [6, 3], [7, 7], [8, 4]]]

A kryptarithm - this is an example of arithmetic, in which all or some of the digits are replaced by letters. The rule must be satisfied: the different letters represent different digits, the identical letters represent the identical digits. See the link

The following procedure, called  Ksolve , solves kryptarithms in which are used four operations ...

It's interesting that every continuous piecewise linear function can be specified by one explicit equation with absolute values​​. The procedure JoggedLine carries out such conversion.

Formal arguments of the procedure: 

A - a list of the coordinates of the vertices of the polyline or the continuous piecewise linear expression defined on the entire real axis.

B (optional) - a point on the left "tail"...

First 6 7 8 9 10 11 Page 8 of 11