Deviation from expected value of mean distance between 2 points 'on ' a sphere

I was trying to verify the mean distance between 2 points in various 3-D and 2-D structures by taking average of multiple random cases. Almost all the time, I was getting a pretty good accuracy except for the case of points on the surface of a sphere. My code uses Gaussian distribution inspired from this answer (see the second most up voted answer)

Here is the python code:

import math as m
from random import uniform as u
sum = 0
for i in range(10000):
    x1 = u(-1, 1)
    y1 = u(-1, 1)
    x2 = u(-1, 1)
    y2 = u(-1, 1)
    z1 = u(-1, 1)
    z2 = u(-1, 1)
    if x1 == y1 == z1 == 0:
        sum += m.sqrt((x2) ** 2 + (y2) ** 2 + (z2) ** 2)
    elif x2 == y2 == z2 == 0:
        sum += m.sqrt((x1) ** 2 + (y1) ** 2 + (z1) ** 2)
    else:
        x1 /= m.sqrt(x1 ** 2 + y1 ** 2 + z1 ** 2)
        y1 /= m.sqrt(x1 ** 2 + y1 ** 2 + z1 ** 2)
        z1 /= m.sqrt(x1 ** 2 + y1 ** 2 + z1 ** 2)
        x2 /= m.sqrt(x2 ** 2 + y2 ** 2 + z2 ** 2)
        y2 /= m.sqrt(x2 ** 2 + y2 ** 2 + z2 ** 2)
        z2 /= m.sqrt(x2 ** 2 + y2 ** 2 + z2 ** 2)
        sum += m.sqrt((x1 - x2) ** 2 + (y1 - y2) ** 2 + (z1-z2) ** 2)
print(sum/10000)

The expected value is 4/3 which is shown here

Arguably the absolute difference is not very large. But the percentage deviation from expected value is around 1% on any run. On the other hand, in all other similar programs with other shapes and same number of random cases, the % deviation is around 0.05% on average.

Also, the value that the code returns is always less than 4/3. This is my major concern.

My guess is that I have implemented the algorithm in a wrong way. Any help is appreciated.