The  procedure  Basis  in  the  Groebner  package  computes  the  reduced Groebner  basis,  but after  this  computation  it  takes  the  primitive  parts  of all polynomials in  the basis. (Unfortunately this is not documented).

So, for
B := Basis(J,T) ;  
the "true"  reduced Groebner basis is
B/~LeadingCoefficient(B, T);

The integral cannot be computed symbolically.

restart;
F:=t -> Int(exp(-1-1/v)*(1-exp(v^2/(2*t^2)))/v^2, v = 0 .. 1):
evalf[30](F(1/10));
              -7.16224515512186261856578554808*10^18 
plot(F, 1..5);

evalf(Int(PDF(0.3*p1+0.7*p2, z, inert), z = 0 .. 1));
                          1.000000000
(only a few seconds).

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].