A bug in sum  occurs for the inner sum.

x:=3:
sum(binomial(x,j)*(1/3)^j*(2/3)^(x-j), j = i+1 .. x);
       1

Even for a simpler sum:

sum(binomial(10,j), j = k .. 10);
                            1024

Using infolevel[sum]:=5;  it appears that the bug is in the summation using hypergeometric functions.

Just replace rotate with 'rotate'.

Your n-dimensional integral

J:= n -> int(add(x[i],i=1..n)^2, [seq(x[i]=0..1,i=1..n)]);

has the exact value

Jexact:= n -> n*(3*n + 1)/12;

seq('J'(n)=Jexact(n), n=1..20);
J(1) = 1/3, J(2) = 7/6, J(3) = 5/2, J(4) = 13/3, J(5) = 20/3, J(6) = 19/2, J(7) = 77/6, J(8) = 50/3, J(9) = 21, J(10) = 155/6, J(11) = 187/6, J(12) = 37, J(13) = 130/3, J(14) = 301/6, J(15) = 115/2, J(16) = 196/3, J(17) = 221/3, J(18) = 165/2, J(19) = 551/6, J(20) = 305/3

(Just in case you want to test the Monte Carlo method for a higher n.)

restart;
J:=Int((-5*ln(x)^4*Pi^4-20*ln(x)^2*Pi^4-8*Pi^4+120*MeijerG([[0, 0], [1, 1, 1]], [[0, 0, 0, 0, 0], []], x^Pi))/(120*Pi^4*(-1+x)^2), x = 4*(1/10) .. 6*(1/10)):
convert(J,Sum):
evalf[15](%);

                      -0.00555168769360751

It's a pity that fsolve fails for the equation erf(x) = 1 - 10^(-20)  when Digits = 5000.
It works when Digits = 4000.
The equation can be solved via

Digits:=5000:
x1:=RootFinding:-Analytic(erf(x) = 1 - 10^(-20), 6-I/10 .. 7+I/10);

or even directly with the good old Newton's method:

restart;
y1:=1-10^(-20):
Digits:=5100:
a:=6.;h:=1. : N:=0:
K:=evalf(sqrt(Pi)/2):
while  abs(h)>1e-5000 do
h:=(erf(a)-y1)*exp(a^2)*K;
a:=a-h:  N:=N+1;
od:  N; a;

Maple computes the limit directly:

a:=sqrt(x): b:=a:
to 3 do b:=sqrt(x+b) od;
limit(b-a, x=infinity);

         1/2

The tutor is designed for beginners. It works for simple cases. Try it e.g. for two nested square roots and it will work.

Note that I simplified your KK supposing p is real.
 

Download K.mw

for i to 3 do
  assign(f[i], x -> P(i,1,1,x))
od;

Edit. Ignore this, it's wrong!

If you use

ln(x) + ln(y) = ln(x y) 

for random complex numbers, the probability to be correct is  75%.
[supposing uniform distribution].

A possible favorable interpretation of the fact that odetest gives 0 for a _Cnn solution of the ODE and <>0 for another constant:

1. In the case of _Cnn  the 0 result says that the solution is correct, but without mentioning that the intervals of validity must be found.
2. In the case of another constant (symbolic or numeric) the result of odetest allows to effectively find the correct interval(s) where the solution is valid.

It would be interesting to know whether this interpretation is correct (i.e. the behavior is intended) or we have a simple bug here.

It's an interesting catch. Vote up!
It seems that the _Cnn  constants are treated in a special manner by odetest.
So, odetest  will return   0  for   sin(sqrt(x) + _C123) ,  but nonzero   for other constants (numeric or symbolic).

I am not able to explain this behavior, but the corect result for odetest should be the nonzero one!  Because sin(sqrt(x)+C/2)  is a solution (assuming x>=0)  only for  sqrt(x)+C/2  in [-Pi/2, Pi/2]  mod 2Pi  
which is equivalent to cos(sqrt(x)+C/2) - abs( cos(sqrt(x)+C/2)  ) = 0  given by odetest.