Dear Mr. Hirnik! Thank you very much for the detailed comment and the time spent! Unfortunately, your code for the solution of this example as an integer linear programming problem has an error. Look at the example and all will be clear:

But if the error would be corrected , the package still does not solve this problem! Look what I got:

The problem was solved after i fundamentally rewrote my previous procedure P1. The basic idea: now I do not form immediately the whole list of possible arrangements of {-1, 1}. For your list of 25 of the eleven-digit numbers, this list will have 2^25  elements, and it's more 32 million. So my previous procedures P and P1 did not work, because in Maple, apparently, there is a limit on the length of the generated lists. Now at first I build two much shorter lists, and then concatenate them in a loop, and  by bubble's method i find the minimum, and the appropriate placement of signs.

The code of the new procedure:

P2:=proc(L)

local It,n1, n2,M0, M1, M2, Min, i, j, S,T ;

It:=proc(K)

[seq([-1,op(K[i])],i=1..nops(K)),seq([1,op(K[i])],i=1..nops(K))];

end proc;

n1:=floor(nops(L)/2); n2:=ceil(nops(L)/2)-1;

M1:=(It@@n1)([[ ]]); M0:=(It@@n2)([[ ]]); M2:=[seq([1,op(M0[i])],i=1..nops(M0))]; 

Min:=infinity;

for i from 1 to nops(M1) do

for j from 1 to nops(M2) do

S:=abs(add(L[k]*M1[i,k],k=1..n1)+add(L[n1+k]*M2[j,k],k=1..n2+1));

if S<Min then Min:=S; T:=[seq(L[k]*M1[i,k],k=1..n1),seq(L[n1+k]*M2[j,k],k=1..n2+1)];

fi;

od; od;

print(Minimum=Min);

T;

end proc:

Solving of the problem:

st:= time():  P2([42924197338, 42184753307, 4832623176, 84884659951, 6402031177, 88749509929,  26310525700, 48147481357, 77529384202, 60286056128, 36115225338,

 49211688701, 45043765556, 98444407624, 37828057713, 42961778036,  97541730413, 32830562321, 23433482548, 72043425495, 6759224825, 89664990345,  37830825668, 77385749629, 25000320804]); 

cat(convert(time() - st,symbol),`  seconds`)=cat(convert(evalf[4]((time() - st)/60),symbol),`  minutes`);

You have two solutions. I would have written this line without "op" by way:

for i from 1 to 2 do `M`[i]:=subs(sol[i],[0,a,0]) end do;

Dear Anton!  There is a standard way of saving the results obtained at each step of the loop of your procedure. Before the loop you create an empty list, and then at each step of the loop adding to this list the results according the scheme:  

1) Before the loop  L:=[ ]:

2) In the loop  L:= [ [х(t[i]), y(t[i]), ... ], op(L) ] , where  х(t[i]), y(t[i]), ...  are the results obtained on i_th step of the work of your loop.

If you speak Russian, you can ask your questions on the forum exponenta.ru

combine(sqrt(a+sqrt(b))*sqrt(a-sqrt(b))-sqrt((a+sqrt(b))*(a-sqrt(b)))) assuming a>0,b>0,a>sqrt(b);

                       0

This is called a parametric form of the equation of the line. Read  ?plot/details 

plot([1, 2-x, [0.466667, t, t=-2..1.5]], x=0..3, -2..2);

A:=[1,2]:

B:= [3,4]:

C:= [op(A), op(B)];

That's right (just missing a colon in a row)! As it should be found two solutions.

Normal vector  (a, b, c)  to your plane can be obtained as the solution of the system

solve({a^2+b^2+c^2=1, a*(0-2)+b*(1-(-1))+c*(-2-0)=0, abs(a+2*b+c)/sqrt(1+2^2+1)=1/2});

Normal vector  (a, b, c)  to your plane can be obtained as the solution of the system

solve({a^2+b^2+c^2=1, 5*a-2*b+5*c=0, abs(a-4*b-8*c)/sqrt(1+4^2+8^2)=sqrt(2)/2});

Your statement  2)  " Vector(MH) perpendicular to direction vector a = (2, 1, -1) of (d)"  is false!

Your problem  is no different from the previous plane (3). The procedure P solves the problem for any points A, B, M and numeric d (points should not be collinear).

Procedure code:

restart:

P:=proc(A,B,M,d)

local P,Sol,L,f,L1;

uses RealDomain, LinearAlgebra, ListTools;

if Equal(simplify(CrossProduct(convert(B,Vector)-convert(A,Vector),convert(M,Vector)-

convert(A,Vector))),<0,0,0>) then error `Points A, B, M should not be collinear`; fi;

P:=a*(x-A[1])+b*(y-A[2])+c*(z-A[3]);

Sol:=[solve({subs(x=B[1],y=B[2],z=B[3],P)=0,a^2+b^2+c^2=1,

abs(subs(x=M[1],y=M[2],z=M[3],P))=d},{a,b,c})];

L:=[seq(<rhs(Sol[i,1]),rhs(Sol[i,2]),rhs(Sol[i,3])>,i=1..nops(Sol))];

f:=(x,y)->Equal(simplify(CrossProduct(x,y)),<0,0,0>);

L1:=[Categorize(f, L)]; 

if nops(L1)=0 then print(`The problem has no solutions`); fi;

if nops(L1)=1 then print(`The problem has 1 solution`); print(collect(subs(a=L1[1][1][1],b=L1[1][1][2],c=L1[1][1][3],P),[x,y,z])=0); fi;

if nops(L1)=2 then print(`The problem has 2 solutions`); print(collect(subs(a=L1[1][1][1],b=L1[1][1][2],c=L1[1][1][3],P),[x,y,z])=0); print(collect(subs(a=L1[2][1][1],b=L1[2][1][2],c=L1[2][1][3],P),[x,y,z])=0); fi; 

end proc;

As an exapmle see solution the problem  plane (3)

P([-1,3,-6],[2,2,-10],[1,-1,7],3);

Similarly, you can find equation of a plane (4)

That's right! Probably it would look better if you write

 [seq(M[i]=(A + u*t)[i],i=1..3)];

or

{seq(M[i]=(A + u*t)[i],i=1..3)};

or

op([seq(M[i]=(A + u*t)[i],i=1..3)]);

Once posted the previous message, noticed in it one essential disadvantage. The code will not work if the points A, B, N1 are collinear. Therefore, rewrote the code using your idea. The code is simpler, and the result , of course, is the same:

with(LinearAlgebra):

N1:=convert(N1,Vector):

n:=o-N1:  e:=n/simplify(Norm(n,2)):

T:=<x,y,z>-N1:

collect(expand(DotProduct(e,T,conjugate=false)),[x,y,z])=0; #  Equation of the plane P

You have 2 planes: ABN1 and ABN2. Find the equation ABN1. The simplest way to use the determinant. The code (the continuation of the previous code):

N1:=convert(N1,Vector):

T:=<x,y,z> - N1: N1-A: N1-B:

P1:=collect(LinearAlgebra[Determinant](< T | N1-A | N1-B >),[x,y,z]):

k:=simplify(sqrt(coeff(P1,x)^2+coeff(P1,y)^2+coeff(P1,z)^2)):

P2:=simplify(collect(P1/k,[x,y,z])):

P:=collect(P2,[x,y,z])=0;  # Equation of the plane P

simplify(coeff(lhs(P),x)^2+coeff(lhs(P),y)^2+coeff(lhs(P),z)^2); # Checking

Similarly, we find a second plane.