I suspect that knowbody knows the exact definition of a singularity used in singular.
Seems to be a combination of "maths", "syntax" and "heuristics".
The singularity at infinity is even more vague.

@tomleslie  It's obviously a bug in solve. The CleanInfty procedure was designed mainly for SolveTools:-SemiAlgebraic where the simple presence of infinity produces an error. Because in SemiAlgebraic the expressions must be polynomials and the variables are real, the presence of infinity is simply not needed (e.g. x < infinity is superfluous).
In solve, a variable is considered real if it appears in an inequality; so here x<infinity could make sense if x is not a subject to other inequality; but here too one may avoid infinity using e.g. x^2 > -1.
In solve the presence of infinity does not always produce an error because the commands from SolveTools are called after some manipulations of the input.

@nm 0<0 was introduced to replace conditions like x < -infinity  or  x=infinity. Actually CleanInfty should not be used; the idea is to avoid infinity in solve because the results seem to be unpredictible.

Note. In your example the second set of r was ignored because CleanInfty was called with 3 actual parameters [it was not defined using $].

@digerdiga 

The things are a bit different, but essentially f(u) = (1-(-1)^floor(u))/(2*u^2)  was considered 0 for u>20;
(Actually  the integral was transformed with a change of variables to [0,1]  and evaluated as 0 in [0,1/20].)
You may see the exact steps executing

infolevel[`evalf/int`]:=5;

@kuwait1 

It seems that you don't read or you are not interested in the provided answers/comments.

Seems to be OK:

restart;
Digits := 15: Z := 4: r := 4:
INT := (n, i) ->
Int(-(2*r*p)^(2*i)*Z*hypergeom([3/2-n-2*i, 1/2], [3/2], Z^2/(Z^2+p^2+r^2))/((Z^2+p^2+r^2)^(1/2)*(p^2+r^2)^(n+2*i-1/2)),
     p = -2 .. 2, method = _d01ajc):
evalf(INT(6,30));
evalf(INT(6,100));
evalf(INT(6,1000));
                     -1.30904529806182 e-15  
                     -5.98504320471557 e-30   
                     -3.61141364348684 e-206   
QQ := proc(n,imax)
local i;
evalf(-2*Pi*r*(add(pochhammer((1/2)*n, i)*pochhammer((1/2)*n+1/2, i)*INT(n, i)/factorial(i)^2, i = 0 .. imax)))
end:
QQ(6,70);
QQ(6,100);
QQ(6,200);
                      0.000211915930003326
                      0.000211915930014491
                      0.000211915930014491

Actually it would be nice to be able to increase Digits, but this is complicated because non-NAG methods are too slow and need fine-tuning for this problem.

A good apprimation would be nice too but it must be done with great care.
 

@digerdiga 

If I were a computer, I'd do the same  :-)

Edit.

= 155685007/232792560 = 0.668771403...      [so, the numeric Int stopped at u=20].

@kuwait1 

You cannot obtain solid numerical results like this. You must find qmax properly (using estimates).
It is not clear how the second QQ was derived. BTW it contains some stange expressions such as
sin(abs(arctan((-1)^s*k/(2*r))))^2   which is constant and recomputed repeatedly.

If you really want to compute the first QQ using your own method then you must learn and use the numerical analysis principles!

@kuwait1 

In this file you only have a finite (triple) sum  containing special functions, so no convergence involved.

If one or two of the inner sums are intended to approximate some hypergeometric functions then you must choose carefully the upper limit of the sums as discussed in my answer (some estimates will be needed for this).
Actually you should say clearly what are you trying to obtain.

@tomleslie 

The content of the files can be viewed:

`#mover(mi("Ξ"),mo("&uminus0;"))`;   # Array ( named Xi_hat - atomic )

TKtm;                                   # expression

@Markiyan Hirnyk 

What about the OP's: " there is the obvious solution X=[e_3,e_1,e_2] where (e_i) is the canonical basis of the space".
Corollary ..., emptyness ...
 

@Markiyan Hirnyk 

The existence was known. It's about the method; it does not work.
If you like your solution (including warnings and/or errors) , it's ok.
 

@Markiyan Hirnyk 

Solid people should check their own "product". Anyway, for your convenience:

restart; with(LinearAlgebra):
X := Matrix(3, 3, symbol = x):
Y := Transpose(X) . X-IdentityMatrix(3):
F := NULL: for i to 3 do for j from i to 3 do F := F, Y[i, j] = 0 end do end do:
sol := [solve([F], explicit)]:
nops(sol);
                               64
solve({x[1, 1] < 1/2, 3/10 < x[2, 3], sol[62][]});
Warning, solutions may have been lost
 ...

@Markiyan Hirnyk 
It produces several "Warning, solutions may have been lost".  And lost they are.

@petit loup 

In my opinion solve or  SolveTools:-SemiAlgebraic  (which is called by solve) are not good tools for your problem. 
Just for fun I have reduced tour toy example to a system

[0 <= 3*a^2+3*b^2-c^2-1, 0 <= 3*a^2-20*a*b+3*b^2+3*c^2-20*c+3]

and  both commands failed for it (I had to interrupt them).

If your specific problem is to find an orthogonal matrix in a neighborhood of a given matrix in GL(n,R)  then I'd suggest to use the Cayley transform (see wiki). This way, the problem reduces to find an antisymmetric matrix in the nbd of a given matrix, which is a simple task (for any n).

@Markiyan Hirnyk 

What explanation are you expecting? We have a good approximate solution and you have checked it. Of course ANY approximate solution could be invalid (in very rare situations).
BTW, you probably know that solve itself uses approximations to isolate roots etc.