В упор не вижу откуда 15. Ладно, зайдем с др стороны

C++ (Qt)
#include <QtWidgets>
 
inline float RandF( void )
{
	return float(qrand()) / RAND_MAX - 0.5f;
}
 
int main(int argc, char *argv[])
{
	qsrand(0);	
	const int num = 1000 * 1024;
	double sum[5] = { 0.0 };
	for (int i = 0; i < num; ++i) {
		QVector3D vec(RandF(), RandF(), RandF());
		vec.normalize();
		double x2 = vec.x() * vec.x();
		double y2 = vec.y() * vec.y();
		double z2 = vec.z() * vec.z();
		sum[0] += x2;
		sum[1] += x2 * x2;
		sum[2] += x2 * y2;
		sum[3] += x2 * z2;
	}
	for (int i = 0; i < 4; ++i)
		qDebug() << (sum[i] / num) << "/" << (num / sum[i]); 
 
	return 0;
}
 

Печатает
Т.е. не 9 и не 15, а какое-то левое 13.07. Предполагаю потому что x, y, z зависимы (хотя и равноправны). Также среднее x2 * x2 совершенно другое, и простые соображения говорят что это верно. Тогда "один делитель на всех" неверно. Выходит методика не обобщается на 4 и более степеней? Или я где-то насвистел в тесте?

Генерация случайного вектора в сферической СК - популярные грабли, она не дает равномерного распределения

C++ (Qt)
QVector3D RandomVec( void )
{
 float theta = float(qrand()) / RAND_MAX * M_PI;   // случайный "угол возвышения"
 float phi = float(qrand()) / RAND_MAX * M_PI * 2;   // случайный "азимут"
 return QVector3D(sin(theta) * cos(phi), sin(theta) * sin(phi), cos(theta)); 
}

"На пальцах" - компонента Z выбрасывается как будто в 2D, а X и Y "придавлены" ее синусом. Поэтому

среднее(x * x) = 1/4
среднее(y * y) = 1/4
среднее(z * z) = 1/2

Впрочем часто нужно не равномерное, а именно такое распределение, такая задачка здесь была

Прошу исполнить. Спасибо

Ну с этим никто не спорит..

Однако если выбрать параметризацию, в которой случайной величиной является проекция ед. вектора на ось z и угол phi:

C++ (Qt)
static const double pi = 4.0*atan(1.0);
 
inline double sign(const double & x) { return (x >= 0.0) ? 1 : -1; }
 
 
int main()
{
    random_device rd;
    mt19937 gen(rd());
 
    uniform_real_distribution<double> phi_dist(0.0, 2.0*pi); // распределение для угла phi
    uniform_real_distribution<double> z_dist(-1.0, 1.0); // распределение для проекции вектора на ось z
 
 
    const int num = 1000 * 1024;
    double sum[5] = { 0.0 };
    for (int i = 0; i < num; ++i)
    {
        double z = z_dist(gen);
        double phi = phi_dist(gen);
 
        double x = sign(z)*sqrt(1.0 - z*z)*cos(phi);
        double y = sign(z)*sqrt(1.0 - z*z)*sin(phi);
 
 
        double x2 = x*x;
        double y2 = y*y;
        double z2 = z*z;
        sum[0] += x2;
        sum[1] += x2 * x2;
        sum[2] += x2 * y2;
        sum[3] += x2 * z2;
    }
    for (int i = 0; i < 4; ++i)
        cout << (sum[i] / num) << endl;
 
    return 0;
}
 

то всё будет корректно. Такое усреднение будет соответствовать интегрированию по телесному углу.

Выкладываю pdf с выводом формулы для среднего от произведения n_i n_j n_k n_l

Ну да, там просто const должен быть, конечно)

Надо посмотреть, что там по этому поводу стандарт говорит)

Я так понимаю, gcc все же умеет constexpr для этих древних legacy функций

А для остальных, я так понял причину - типа так как функции очень старые и раньше не были constexpr, то и не будем мы это менять, как бы чего не вышло? Но что может измениться, косинусы углов вдруг станет модно считать по другому Имхо, надуманная причина "оставить как есть"...

20) Есть 2 шарика одинакового размера и материала (напр бильярдных)

C++ (Qt)
struct Ball {
 QVector3D m_position;
 QVector3D m_speed;
 static float radius;
};

Написать ф-цию

C++ (Qt)
bool BallsCollided( Ball & ball1, Ball & ball2 )
{
// Если шарики столкнулись то вернуть true
// и рассчитать новые скорости шариков
}