@_Maxim_ 

This is simply a bug. In `simplify/commonpow`  the case a=0  in  a^b  was forgotten such that

simplify(0^(k-1)) assuming k>1;

produces an error.

@kainmuth 

It's not correct.

int(w(x-y), y=0..2*Pi) assuming 0<x, x<2*Pi;
     infinity

(obvious without Maple).
By Fubini, the double integral is +oo.

( It's like trying to use Newton-Leibniz for 1/|x|,  x in [-1,1]. )

@Adam Ledger 

Try a simple experiment.
n:=3703703951851853;

a) The direct computation is out of the question.
b) Now write yourself the simplest procedure for prime decomposition using only irem
You will obtain easily n=p*q  ==> phi(n) = ...

@mnovaes 

for x in Elements(S3) do 
  print('x'=x, 
        orbits=Orbits( PermutationGroup(x, degree=3 )),
        numorbits=numelems(Orbits(PermutationGroup(x, degree=3)))  ) 
end do;

Download yourU.mw

@vv 

OK, I was looking only for symmetric solutions and I see now that your U is not symmetric.
Changing

U:=< u11,u12,0,  0;     # look for a block-diagonal U, such as M and C
     u21,u22,0,  0;
     0,  0,  u33,u34;
     0,  0,  u43,u44>;

==> 64 solutions, some of them depending on a parameter.
Most probably your U is one of them. Just check.

@Adam Ledger 

It is faster, and maybe e.g. the user only wants to have the duplicates near each other.

@Kitonum 

A not so obvious fact is that the sequence X(n) is dense in the interval [1,2], i.e. for each 1 <= t <= 2  there is a subsequence converging to t.

@das1404 

The lines can be removed. E,g.
display(polygon(A), axes=none, color=gold, style=patchnogrid);

@acer 

The real option in solve does not seem to be reliable.

restart;
solve((x^2-1)*(exp(x)+1),{x});
                 {x = I Pi}, {x = 1}, {x = -1}
solve((x^2-1)*(exp(x)+1),{x},real);
                 {x = I Pi}, {x = 1}, {x = -1}
solve((x^2-1)*(exp(x)+1),{x},real,allsolutions);
                 {x = I Pi}, {x = 1}, {x = -1}
solve((x^2-1)*(exp(x)+1),{x},allsolutions);
           {x = I Pi + 2 I Pi _Z1}, {x = 1}, {x = -1}

identify is far for being reliable. It fails here for some polynomials.

@Mariusz Iwaniuk 

@nyarko 

After all, sqrt is the simplest function after polynomials.
Just think: can the sqrt be removed in  2 + sqrt(3) ?

The method was posted just for fun and to show a (very inefficient!) possibility. The polygonal approach is much more efficient but designing nice letter shapes is time consuming. The ideal solution would be to have an angle parameter in textplot.

Or, including the asymptotics:

with(IntegrationTools):
A:=Int(exp(-x*t)/sqrt(t^2+t),t=1..infinity):
to 6 do
  A:=Parts(A,GetIntegrand(A)/exp(-x*t)) assuming x>0;
od:
eval(A,Int=0);

@rlopez 

@Math Pi Euler 

Probably you have included the prompt.
Here is the worksheet: asy.mw