@Kitonum  

My first program was unsuccessful. I forgot that break keyword only interrupts the inner loop, but the program continues to run. Therefore returned not the first but the last solution.
One solution (the first one) can be found immediately when the program as a procedure to issue and uses  return keyword:

restart;
P := proc ()
local a, b, c, d, e, f, g, h;
for c from -10 to 10 do
for d from -10 to 10 do
for e from -10 to 10 do
for f from -10 to 10 do
for g from -10 to 10 do
for h from -10 to 10 do
a := c+2*d+4*e+6*f+10*g+22*h-16; b := -2*c-3*d-5*e-7*f-11*g-23*h+20;
if -10 <= a and a <= 10 and -10 <= b and b <= 10 then return [a, b, c, d, e, f, g, h]
end if end do end do end do end do end do end do
end proc:

P();

                                               [10, 10, -10, -10, -10, 8, 7, -1]

@Mikhail Drugov  Yes,  I think there is no automatic way.

@brian bovril  Of course I can. This is really easy.

Let the member of the sequence, standing on the  n-th position is  m. Then it is obvious that

1+2+...+(m-1)< n <=1+2+...+m   or   (m-1)*m/2< n <=m*(m+1)/2

Then, solving this double inequality for m> = 1, we obtain

-(1/2) + 1/2*sqrt(1 + 8*n) <= m < 1/2 + 1/2*sqrt(1 + 8*n)

Because left and right parts of the solution differ by 1, we immediately obtain the formula.

The question arises - what of it relates to Maple?

your question more specifically?

@Carl Love   

I think that the use of functions  ln  and  evalf  is impractical,  because it slows down the calculation. As discrepancy for a given partition I take the simple sum of absolute differences.

Compare two approaches for the small set {$ 1 .. 15}:

restart;

ts := time():

S := {$ 1 .. 15}:
AllP := [seq(P, P = Iterator:-SetPartitions(S, [[5, 3]], compile = false))]:
lnp := evalf(ln((`*`(S[]))^(1/3))):
Var := proc (P::(list(set)))
local r; r := evalf(`+`(map(proc (b) options operator, arrow; abs(ln(`*`(b[]))-lnp) end proc, P)[])) end proc: 
sort(AllP, (x, y) ->Var(x) < Var(y) )[1];

time()-ts;

                       [{1, 8, 9, 10, 15}, {2, 5, 7, 12, 13}, {3, 4, 6, 11, 14}]
                                                             322.328

restart;

ts := time():

S := {$ 1 .. 15}: Var := infinity:

T := combinat[choose](S, 5):

for t in T do

S1 := `minus`(S, t):

U := combinat[choose](S1, 5):

for u in U do

S2 := `minus`(S1, u):

Var1 := abs(convert(t, `*`)-convert(u, `*`))+abs(convert(u, `*`)-convert(S2, `*`))+abs(convert(t, `*`)-convert(S2, `*`)):

if Var1 < Var then Var := Var1:

L := {Var, {S2, t, u}} end if

end do end do:

L;

time()-ts;

                                {576, {{1, 4, 13, 14, 15}, {2, 5, 9, 10, 12}, {3, 6, 7, 8, 11}}}
                                                                                 14.516

The results are different, but the discrepancy is actually one and the same, that just shows the non-uniqueness of solution:

[{1, 8, 9, 10, 15}, {2, 5, 7, 12, 13}, {3, 4, 6, 11, 14}];

abs(convert(%[1], `*`)-convert(%[2], `*`))+abs(convert(%[2], `*`)-convert(%[3], `*`))+abs(convert(%[3], `*`)-convert(%[1], `*`));

                                                                576

@Markiyan Hirnyk  In my example, do the rotation by the angle  arcsin (a/2)  about the axis Oy and the homothety with coefficient  sqrt (3) and the center at the origin. Then you will get the same coordinates as in Wiki.

@Markiyan Hirnyk 

1. Both procedures are very simple conceptually and based on brute force technique of two or three vertices and then checking (using geometric properties listed above). So I do not think there is a need for separate comments.

2. To be more precise, that used O (n ^ 3) algorithm, because  binomial(n, 3)=n*(n-1)*(n-2)/6 . Of course, in terms of time, these procedures are not optimal and appropriate for relatively small polygons (n<100). But their advantages in simplicity and exact solutions. Therefore they can be used in the proofs in geometry problems.

3. I do not agree with this statement "Exact solutions and close-form-expressions are from somewhat outdated math" . The advantage of an exact solution at least in the fact that from it is easy to get  an approximate solution with any precision, but not vice versa.

@toandhsp  I was looking for it of right triangles with legs   x+3  and  y+3  and the hypotenuse  x+y . Radius of the inscribed circle is  3  for any such triangle.

isolve((x+3)^2+(y+3)^2=(x+y)^2):
op(select(s->rhs(s[1])>0 and rhs(s[2])>0,[%]));

{x = 4, y = 21}, {x = 5, y = 12}, {x = 6, y = 9}, {x = 9, y = 6}, {x = 12, y = 5}, {x = 21, y = 4}

@toandhsp   [[4,-2], [4,7], [-8,-2]]

see  http://www.mapleprimes.com/posts/200122-The-Largest-Incircle-And-The-Smallest-Covering-Circle-

@Markiyan Hirnyk  Nice and simple solution. Now I am preparing two procedures in which  both problems will be solved exactly (symbolically). For example, the exact answer to this example is

[[57/22, 43/22], 1/242*1625^(1/2)*242^(1/2)];
evalf[15](%);

@sainimu78  Syntax error! Under the symbol  ∑  the command  add  should not be used.

@isbah It is impossible to obtain solution explicitly in the form of formulas. The reason is that the solutions are the roots of the polynomial of high degree  RootOf(...) . It is well known that if the degree of a polynomial more than 4, then the roots are not expressed as a function of its coefficients. But the procedure

f:=a->fsolve({0.5704+0.5211*exp(-0.6569*x)-0.03849*y=a, 0.561+0.5284*exp(-0.6479*x)-0.03929*y=a}):

can be used also successfully as  explicit formulas, for examples, to find the values ​​at certain points, build plots, etc.

@Carl Love  In initial task, we don't need use these commands, because they slow down the program and we can do the work without them.

Example:

A:=[1,4,8]: B:=[8,9,12]:
type(sqrt(A[1]^2+A[2]^2+A[3]^2), posint);
type(sqrt(B[1]^2+B[2]^2+B[3]^2), posint);
type(sqrt((A[1]-B[1])^2+(A[2]-B[2])^2+(A[3]-B[3])^2), posint);

                                               true

                                               true

                                               false

The main disadvantages of the source code - it is slow and duplicates occurs.