@alimahmood 

n is the dimension of the square matrix Z. Maple cannot use a symbolic n, so, n must be specified (any positive integer).

Probably y(...) means  y*(...)  but  sigma[s] = ?

@AmirHosein Sadeghimanesh 

You should be aware that it is not correct to use Riemann sums for improper integrals. E.g.

Int(sin(1/t^6)/t^6, t = 0 .. 1);

Here the Riemann sums diverge but the integral converges (not absolutely).

@AmirHosein Sadeghimanesh 

I assume that you have the absolute value in the integral. Then, by Fubini-Tonelli one may swap the integrals and the inner one is +oo.
If you remove the absolute value, some computations are needed and it seems that the iterated integral still diverges (but I had not the time for complete analysis in this case).

@Rouben Rostamian  

_Z is not created at top level. It could have been local.

@Carl Love 

I agree. But it is not my solution, it is OP's! I don't know the reason why he does it.

@Carl Love 

Of course if a is not too large.
If a has say 100 digits, the last 10 digits of a^b  can be obtained, but we will be able to detect usually only the divisors a^n (n <= b). E.g. it is not possible to find the last 10 digits of the least proper divisor of a^b.

Probably the only correct approach for an "exact rank" would be to compute the singular values using range arithmetic and return FAIL if one of the singular values is in an interval containing 0. I don't think that this will be possible in the near future.

@acer 

The rank as a function of the matrix is not continuous, so anyway an optimal approach does not exist.

@mqb 

You are right the rank of a vector is of course 0 or 1. I have corrected the answer.

@ganelon 

Unfortunately no. As well known, finding a single divisor (or even its existence) for a huge number could be practically impossible.

@_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;