C言語でガウス分布を実装したい

Question

前提・実現したいこと

C言語でガウス分布（正規分布）を出すような関数を作っているのですが，上手く作れません。

https://omitakahiro.github.io/random/random_variables_generation.html#Prepare_MT
こちらのサイトを参考にして次のコードを書きました。しかし，標準正規分布の場合にして，例えばn=10000の場合に対してガウス分布を通したものに対して合計値を出力しても，その値が0になりません。
当然その値が0にはならないと思いますが，限りなく0に近い値になるはずです。少なくとも，後述するように-135程度の比較的大きい値になるとは考えられないです。
何がおかしいのか分かりません。型はfloatで，要素数nのベクトルをガウス分布で初期化する方法を教えてください。
よろしくお願いします。

該当のソースコード

#include <stdio.h>
#include <stdlib.h>
#include <math.h>  
#include <time.h>

//[0:1]の一様乱数を生成
float uniform( void ){
    return rand()/(RAND_MAX+1.0);
}

//要素数nのoを正規分布で初期化
void Gaussian_init(int n, float *o, float mu, float sigma){
    int i;
    for(i=0;i<n;i++){
        float z =sqrt( -2.0*log(uniform()) ) * cos( 2.0*M_PI*uniform() );
        o[i] =  sigma * z + mu;
    }
}

int main() {
    srand(time(NULL));
    int n = 10000;
    float *x = malloc(sizeof(float) * n);
    float sum = 0.0;

    Gaussian_init(n, x, 0.0, 1.0);
    for (int i = 0; i < n; i++) {
        sum += x[i];
    }
    printf("sum = %lf\n", sum);

    return 0;
}

出力結果

sum = -135.292908

Answer 1

和が0に近づくというのが間違いだと思います。

互いに独立なn個の確率変数X1, X2, ..., XnがN(μ, σ^2)に従うとすると,
正規分布の再生性から、確率変数の和X1+X2+...+XnはN(nμ, nσ^2)に従います。
今回μ=0, σ=1, n=10000であるため,
N(0, 10000)となります。
なのでsumが±100(1σ区間)に収まる確率が68 %, ±200(2σ区間)に収まる確率が95 %なので,
-135となるのはおおよそ当たり前だと言えます。
個数の1/2乗に比例して標準偏差（ばらつき)が大きくなっていくので、和が限りなく0に近づくというのは間違いです。

一方、sumをnで割った値は0に近づいていきます。
なぜなら、(X1+X2+...+Xn)/nはN(μ, σ^2/n)に従うことから、標準偏差が個数の1/2乗に反比例し
個数を増やせば増やすほどばらつきが小さくなるためです。

Answer 2

float uniform()の返す値が，rand()/(RAND_MAX+1.0)であり，
分母の値が，rand()が取りえる最大値よりも1.0だけ多いことが，
結果として計算されるガウス乱数値を負の側にわずかに偏らせるのではないでしょうか．
（ uniform()の返す値の範囲が　[0:1] ではなく [0:1) になっている ）

で，それが，ざっくりとn=10000倍したスケールで最終結果に出てきているのでは．