The url gives me a 404 error.

Thank you very much Neil for these suggestions.

Meanwhile, I also realized that for the end-result (the conditional probability that more than 1 lineages were present at the time of the mitochondrial Eve about 200'000 years ago, knowing that at present there is only 1) I don't need to calculate 200'000/20=10'000 generations because it's converging rapidly only after 20 generations or so. In other words, calculating more than 20 generations is not necessary because the result doesn't change significantly.

What I still need is to fit the code with paleoanthropological data, especially the probability that a mother has 0 daughters. This is not the same than in modern societies. My guess is that this is only conditioned by sterility. This is very important to know because if its probability is zero, lineages do not die out at all... I'll let you know when I finish my article.

Yes, I have VisualStudio installed. And to my surprise your code compiled flawlessly... Normally, when I did this in the past, it always refused some code that I needed to render C-compatible.

The next step is now to find a probability distribution for the amount of daughters because in reality this is not a uniform probability as rendered with RandomInteger. A mother, for instance, that has 0 daughters is much less probable than if she has 3, in any case this is true for hunter-gatherers. I tested it with RandomVariate and the binomial distribution, which yields quite a different end-result than with RandomInteger.

As long as you have a C compiler installed on your system this is not difficult.

Just check if one is availible using

Needs["CCompilerDriver`"]
DefaultCCompiler[]
CCompilers[Full]

I quickly checked and compiling with CompilationTarget -> "C" of the compiled function in my previous post makes it faster.

(d3 = distributionF1C[3, 5, 3, 10000, 10]) // RepeatedTiming
(d3 = distributionF1C2[3, 5, 3, 10000, 10]) // RepeatedTiming

Out[17]= {0.12, {{1., 0.2657}, {2., 0.3424}, {3., 0.3544}, {4., 
   0.3604}, {5., 0.3613}, {6., 0.3611}, {7., 0.3586}, {8., 
   0.3576}, {9., 0.3574}}}

Out[18]= {0.083, {{1., 0.2731}, {2., 0.3305}, {3., 0.3393}, {4., 
   0.3391}, {5., 0.3397}, {6., 0.3402}, {7., 0.3395}, {8., 
   0.3399}, {9., 0.3393}}}

Select is generally slow, for small lists DeleteCases is faster for larger lists Pick will win.

time = Transpose[1000000 Table[
     b = RandomInteger[5, i];
     t1 = First@RepeatedTiming@Select[b, # > 0 &];
     t2 = First@RepeatedTiming@DeleteCases[b, 0];
     t3 = First@RepeatedTiming@Pick[b, Unitize[b], 1];
     {t1, t2, t3}
     , {i, 1, 20, 1}]];
ListLinePlot[time, PlotRange -> Full, 
 PlotLegends -> {"Select", "DeleteCases", "Pick"}, 
 AxesLabel -> {"N", 
   "time [\!\(\*SuperscriptBox[\(10\), \(-3\)]\) ms]"}]

time = Transpose[1000 Table[
     b = RandomInteger[5, Round[10^i]];
     t1 = First@RepeatedTiming@Select[b, # > 0 &];
     t2 = First@RepeatedTiming@DeleteCases[b, 0];
     t3 = First@RepeatedTiming@Pick[b, Unitize[b], 1];
     {t1, t2, t3}
     , {i, 1, 5, 1}]];
ListLinePlot[time, PlotRange -> Full, 
 PlotLegends -> {"Select", "DeleteCases", "Pick"}, 
 AxesLabel -> {"\!\(\*SuperscriptBox[\(10\), \(N\)]\)", "time [ms]"}]

Also i think the function call to r[x_] is not needed but not sure how much difference that makes. But compile still will gain you a lot.

distributionF1[daugthers_, mothers_, lin_, itt_, gen_] := 
 Block[{w, g, s, v},
  w = ConstantArray[0, gen - 1];

  Do[
   v = 1;
   g = DeleteCases[RandomInteger[daugthers, mothers], 0];
   s = Length[g];
   If[s == lin, w[[v]]++];

   While[s != 0,
    v++;
    If[v == gen, Break[]];
    g = DeleteCases[Total[Map[RandomInteger[daugthers, #] &, g], {2}],
       0];
    s = Length[g];
    If[s == lin, w[[v]]++];
    ];

   , {itt}];

  Thread[{Range[gen - 1], N[w/itt]}]
  ]

distributionF1C = 
  Compile[{{daugthers, _Integer}, {mothers, _Integer}, {lin, _Integer}, {itt, _Integer}, {gen, _Integer}},
   Block[{w, g, s, v},
    w = ConstantArray[0, gen - 1];
    Do[
     v = 1;
     g = DeleteCases[RandomInteger[daugthers, mothers], 0];
     s = Length[g];
     If[s == lin, w[[v]]++];

     While[s != 0,
      v++;
      If[v == gen, Break[]];
      g = DeleteCases[Map[Total[RandomInteger[daugthers, #]] &, g], 0];
      s = Length[g];
      If[s == lin, w[[v]]++];
      ];
     , {itt}];

    Thread[{Range[gen - 1], N[w/itt]}]
    ]
   ];

The gain is small with DeleteCases, but with compile its faster.

In[963]:= (d1 = distribution[3, 5, 3, 10000, 10]) // RepeatedTiming
(d2 = distributionF1[3, 5, 3, 10000, 10]) // RepeatedTiming
(d3 = distributionF1C[3, 5, 3, 10000, 10]) // RepeatedTiming

ListLinePlot[{d1, d2, d3}, PlotRange -> Full]

Out[963]= {1.09, {{1, 0.2633}, {2, 0.326}, {3, 0.3396}, {4, 
   0.3459}, {5, 0.3469}, {6, 0.346}, {7, 0.345}, {8, 0.3437}, {9, 
   0.3433}}}

Out[964]= {0.86, {{1, 0.2629}, {2, 0.3235}, {3, 0.3452}, {4, 
   0.3438}, {5, 0.345}, {6, 0.3447}, {7, 0.344}, {8, 0.3426}, {9, 
   0.3422}}}

Out[965]= {0.12, {{1., 0.2633}, {2., 0.3296}, {3., 0.3428}, {4., 
   0.3451}, {5., 0.3442}, {6., 0.3409}, {7., 0.3414}, {8., 
   0.342}, {9., 0.3428}}}