@Markiyan Hirnyk 

When I add a new vertex to the already created route, I always check that the new vertex is not the same as the previous vertices of this route (in this fragment of the code):

if convert([seq(r[i]<>p, i=1..k-1)], `and`) then M:=[op(M),[op(r), p]]:  fi:

@Markiyan Hirnyk 

When I add a new vertex to the already created route, I always check that the new vertex is not the same as the previous vertices of this route (in this fragment of the code):

if convert([seq(r[i]<>p, i=1..k-1)], `and`) then M:=[op(M),[op(r), p]]:  fi:

@Markiyan Hirnyk

The subprocedure  Rule  is of fundamental importance in procedure  Routes. It allows you to extend the already created list of routes  K   one more step. Applying  Rule  to the list of all the vertices  L , we get a list of all routes of length 1. Applying the procedure  Rule  to the previous list, we get all routes of length 2, and so on.

@Markiyan Hirnyk

The subprocedure  Rule  is of fundamental importance in procedure  Routes. It allows you to extend the already created list of routes  K   one more step. Applying  Rule  to the list of all the vertices  L , we get a list of all routes of length 1. Applying the procedure  Rule  to the previous list, we get all routes of length 2, and so on.

nops(L1)  means the number of elements of the list  L1 . 

Fragment of the code

for i from 1 to n do

U[1, i]:=L1[i]; U[2, i]:=L2[i]; U[i, 1]:=M1[i]; U[i, 2]:=M2[i];

od;

 specifies the first two lines and columns of the matrix  U .

nops(L1)  means the number of elements of the list  L1 . 

Fragment of the code

for i from 1 to n do

U[1, i]:=L1[i]; U[2, i]:=L2[i]; U[i, 1]:=M1[i]; U[i, 2]:=M2[i];

od;

 specifies the first two lines and columns of the matrix  U .

Axel Vogt, I do not understand this statement :

"I think for that example over the _Reals_ the solution is

piecewise(x < 0, 4-4/365*(-11315*x+133225)^(1/2), 0<=x, 0) ;" 

A number or set of numbers  are  the roots of the equation, but you write a function that depends on x.

Your example in the real domain can be easily solved as above by the change:

restart;

eq1:= 2-2*(-31/365*x+1)^(1/2)-2*(2-2*(-31/365*x+1)^(1/2)-31/365*x)^(1/2)=0:

eq2:=-31/365*x+1=t^2:

simplify(algsubs(eq2, eq1)) assuming t>=0;

sol:=solve(%, useassumptions) assuming t>=0;

L:=convert(sol, list);

x:=unapply(solve(eq2, x), t);

'RealRange'(x(L[2]), x(L[1]));  # All roots of the equation eq1 in real domain

This equation can be solved by hand, if you select a perfect square under external root. If the variable  x  does not belong to the segment from 0 to 365/31, the complex numbers (in particular negative numbers) will be under the roots. If we keep in mind the principal value of the square root, it is easy to check that any number  x>=0  (in particular, x=12) is the root. Any other numbers (real or complex, including negative numbers) are not roots.

Axel Vogt, I do not understand this statement :

"I think for that example over the _Reals_ the solution is

piecewise(x < 0, 4-4/365*(-11315*x+133225)^(1/2), 0<=x, 0) ;" 

A number or set of numbers  are  the roots of the equation, but you write a function that depends on x.

Your example in the real domain can be easily solved as above by the change:

restart;

eq1:= 2-2*(-31/365*x+1)^(1/2)-2*(2-2*(-31/365*x+1)^(1/2)-31/365*x)^(1/2)=0:

eq2:=-31/365*x+1=t^2:

simplify(algsubs(eq2, eq1)) assuming t>=0;

sol:=solve(%, useassumptions) assuming t>=0;

L:=convert(sol, list);

x:=unapply(solve(eq2, x), t);

'RealRange'(x(L[2]), x(L[1]));  # All roots of the equation eq1 in real domain

This equation can be solved by hand, if you select a perfect square under external root. If the variable  x  does not belong to the segment from 0 to 365/31, the complex numbers (in particular negative numbers) will be under the roots. If we keep in mind the principal value of the square root, it is easy to check that any number  x>=0  (in particular, x=12) is the root. Any other numbers (real or complex, including negative numbers) are not roots.

VolMike ,

you are right!

VolMike ,

you are right!

@Carl Love 

Thanks for your interest and  improvements in the algorithm.

Comparison of timing of all three variants for a large number  10^11  ( Factoring1  is Joe Riel's modification).

t:=time():

Factoring(10^11):

time()-t;

     53.290

t:=time():

Factoring1(10^11):

time()-t;

      8.923

t:=time():

Factorings(10^11):

time()-t;

      3.666

@Markiyan Hirnyk

Thanks for the shown attempts! Unfortunately, the result is incorrect and does not satisfy the assigned objective . Waiting for your third and subsequent attempts.

@Markiyan Hirnyk

Thanks for the shown attempts! Unfortunately, the result is incorrect and does not satisfy the assigned objective . Waiting for your third and subsequent attempts.

@Markiyan Hirnyk

If you know of another way to solve this problem (for example, as integer programming problem), please show your code.

@Markiyan Hirnyk

If you know of another way to solve this problem (for example, as integer programming problem), please show your code.