@sand15 

assuming positive  means that all the variables are assumed to be >0.
On my system (Maple 2016.2, 64 bit) I get FAIL  for     is( r^2 - r*s + s^2> 0) assuming r >0, s>0;

positive means real and >0. So,  is(I, positive)  is of course false (not being real).

r^2+s^2 > 0   is simple for Maple (sum of two positive numbers); s^2 - s*r + r^2  seems to be too complicated for the moment.

Edit (concerning "simple conditions"):

is( (r - s)^2 + r^2  + s^2  > 0)  assuming positive;
                              FAIL
is( (r - s/2)^2 + 3/4*s^2  > 0) assuming positive;
                              FAIL
is( (r - s/2)^2 + 3/4*s^2  > 0) assuming positive, r>s/2;
                              true

@mehdibaghaee 

My question was about the type of the problem. Where/how  those 10^6 variables are they used?
For instance, do you have a system of equations with 10^6 unknowns?

@mehdibaghaee 

It would be interesting to know the type of the problem which needs 10^6 variables and can be solved.

@Markiyan Hirnyk 

But you have spent a rather long time with this problem as I can see. So, selectively busy.

@Markiyan Hirnyk
It is obvious (as I have already said) that the procedure cannot print an infinity of points.
And I do not have enough spare time to analyse more deeply the problem (for which there are many articles in the literature).
But you are free (and invited) to write a better procedure. I will vote it up!

@Zeineb 

Using distributions (such as Dirac) is a delicate task. They are not actually functions but Maple allows algebraic operations with them. It is the user's responsibility to ckeck that these operations have a mathematical sense (for example, usually it is not possible to multiply a non-regular distribution by a function which is not C^oo).

For t=7Pi/6 AND 11Pi/6

@Markiyan Hirnyk 

@Markiyan Hirnyk 

Why don't you try to understand the problem first?
t runs in the interval [0,2Pi).
For example, the point [2,0] is double (in this approach): it is obtained for t=0 AND for t=Pi.

@Markiyan Hirnyk 

In a Lissajous curve m,n should be relatively prime i.e. igcd(m,n)=1. Otherwise all the points are double/multiple and the procedure shows only the "extreme" ones. So, use instead dblpoints(1,2,2,3).

For a non-trivial procedure see the chapter "Advanced examples" in

Char B.W. et al - First Leaves, A Tutorial Introduction to Maple V, Springer 1992

It is very fast and works for huge numbers.

@Carl Love 

Hello Carl, glad to "see" you here again.
I think that OP's confusion is due to the fact that t was omitted in 'maximize'.
So, probably it is better to write:

f:= cos(2*t/m) + cos(2*(t+5)/m):
plot('maximize'(f,t), m= 1..10);

@Christopher2222 

It is not the loop. If you recompute several type evalf(add(...)) ==> random results.
And it is not the add, the same thing happens for `+`(seq(...)).
The problem is the randomness, not necessarily the 3 digits!
This is not nice.

@FlamingFox 

In the first procedure (=J) the condition 1<n<k  is not correct in 1D math and not needed in 2D; it should be simply n<k  as in J.
In the last return you have too many ")".

@FlamingFox 

Then it would be nice to post your corrected procedure.

@Kitonum 

Of course, k>1. I have not included a check for this.
Or, one may insert
if k=1 then return n-1 fi;
but I prefer to keep it as it is.