Distances between random points in a hypercube and statistics of exponents TL;DR: Why is $\text{avg}\left(|a-b|^k\right)=\frac{2}{(k+1)(k+2)}$? I.e. for $k=2$, as for finding average Euclidian distances, the result is $\frac{1}{6}$?
I've been reading a book about "Corobs," and having little understanding of math and statistics, I've been wondering why this particular (paraphrased) formula is true:

In a unit $n$-cube, the distance between two random points quickly approaches $\sqrt{\frac{n}{6}}$ as $n$ grows larger than 20.

In the shower one day, it dawned on me that since the distance formula adds up a bunch of positive numbers, by the law of averages, it makes sense that the expected distance is $\sqrt{\frac{n}{C}}$, where $C$ is some constant. I was trying to remember some simulations I've run. The average absolute difference between two random numbers (range=0-1) is $\frac{1}{3}$. Squaring that gives $\frac{1}{9}$. 
So why isn't the average distance $\sqrt{(\frac{1}{3})^2+...+(\frac{1}{3})^2}$? I realized it's because we're averaging the total distance, not the differences in each dimension -- which means that we need to simulate the average square of the differences. So rather than figuring $\text{avg}\left(a-b\right)^2$, figure $\text{avg}\left((a-b)^2\right)$.
I found that the average square of the difference between two random numbers between 0 and 1 is $\frac{1}{6}$, which makes sense when plugged into the distance formula: $\sqrt{\frac{1}{6}+....+\frac{1}{6}} = \sqrt{\frac{n}{6}}$. However, I still don't understand quite why it's $\frac{1}{6}$ for each dimension and not $\frac{1}{9}$.
So I did some more simulations. I quickly noticed a strange pattern for expected differences raised to some power $k$:
from numpy.random import *

def avg(numbers):
    return float(sum(numbers)/max(1, len(numbers)))

for k in range(1,13):
    print("k="+str(k)+": 1/"+str(
        round(1.0/avg([abs(random()-random())**k 
            for _ in range(1000000)]))))

Result (I added my notes in parentheses):
k=1: 1/3 (3=1*3)
k=2: 1/6 (6=2*3)
k=3: 1/10 (10=2*5)
k=4: 1/15 (15=3*5)
k=5: 1/21 (21=3*7)
k=6: 1/28 (28=4*7)
k=7: 1/36 (36=4*9)
k=8: 1/45 (45=5*9)
k=9: 1/55 (55=5*11)
k=10: 1/66 (66=6*11)
k=11: 1/78 (78=6*13)
k=12: 1/91 (91=7*13)
...

It appears that for any power $k$, the result is $\frac{2}{(k+1)(k+2)}$.
Is this intuitive to anyone? Where does this pattern come from, and what's the math behind it?
 A: The problem of measuring the average of an Euclidian distance $d_n$ between two random points in a n-dimensional hypercube is called a hyper cube line picking, and has no analytical solution for high dimensions. 
At first, I thought that there could be a recurrent solution like the one for volumes, but it turns out that it's not as easy.
The average was obtained for low dimensions, for instance:


*

*n=1, it's $E[d_1]=1/3$

*n=2, it's a nonintuitive fraction $E[d_2]=\frac{2 + \sqrt 2 + 5\mathrm{arcsinh}(1)} {15}\approx 0.52$.


Here's a plot from Wolfram web site with mean distance as a function of $n$:

For higher dimensions, the lower and upper bounds were found in Anderssen, R.S. & Brent, R.P. & Daley, D & Moran, P.A.P.. (1976). Concerning ∫ 0 1 ⋯∫ 0 1 (χ 1 2 +⋯+χ k 2 ) 1/2 dχ 1 ⋯dχ k and a Taylor series method. SIAM Journal on Applied Mathematics. 30. 10.1137/0130003. They used Taylor series expansion in this paper. Your book refers to this asymptotic upper bound: $E[d_n]\le\sqrt{n/6}$, when $n\to\infty$. This is not the mean distance itself.
Square distance
Note, that the average square distance is a very different animal. Consider this:
$$d_n^2=\sum_{i=1}^n(a_i-b_i)^2$$
If you want to take the expectation of this square distance, then it separates into a sum of simple integrals as follows:
$$E[d_n^2]=\int_0^1da_1\int_0^1db_1\int_0^1da_2\int_0^1db_2
\dots
\int_0^1da_n\int_0^1db_nf(a_1,b_1,a_2,b_2\dots,a_n,b_n)d_n^2=
\\
\sum_{i=1}^n \int_0^1da_i\int_0^1db_i f(a_i,b_i)(a_i-b_i)^2$$
This separation of one volume integral into a sum of dimension integrals for the distance itself is prevented by the square root operation which binds all the coordinate squares into one expression (nonlinearly):$\int_{V}dV\sqrt{\sum_{i=1}^n(a_i-b_i)^2}$
For n dimensions the average square distance is $n/3$
Power of absolute sequence
Your sequence for power of a difference of two random numbers has nothing to do with any of the above. Here's how you calculate it:
$$E[|a_1-b_1|^n]=2\int_0^1 (1-x) x^n dx=2\int_0^1 (x^n-x^{n+1})  dx=\\
2/(n+1)x^{n+1}-2/(n+2)x^{n+2}\mid_0^1=2(n+2-(n+1)/(n+1)/(n+2))=\\
\frac 2 {(n+1)(n+2)}$$
This is also the reciprocal of the binomial coefficient: $\frac 1 {C(n+2,n)}$, where $C(n,k)=\frac{n!}{k!(n-k)!}$
You can plug $n=12$ to see that it's $\frac 1 {78}$ as you obtained from simulation.
It's very difficult to be in this field armed only with computers. You need a little bit of mathematics if you don't want to feel lost all the time.
