Distribution of 0 digits in factorial series

Distribution of 0 digits in factorial series
You are encouraged to solve this task according to the task description, using any language you may know.

Large Factorials and the Distribution of '0' in base 10 digits.

About the task

We can see that some features of factorial numbers (the series of numbers 1!, 2!, 3!, ...) come about because such numbers are the product of a series of counting numbers, and so those products have predictable factors. For example, all factorials above 1! are even numbers, since they have 2 as a factor. Similarly, all factorials from 5! up end in a 0, because they have 5 and 2 as factors, and thus have 10 as a factor. In fact, the factorial integers add another 0 at the end of the factorial for every step of 5 upward: 5! = 120, 10! = 3628800, 15! = 1307674368000, 16! = 20922789888000 and so on.

Because factorial numbers, which quickly become quite large, continue to have another terminal 0 on the right hand side of the number for every factor of 5 added to the factorial product, one might think that the proportion of zeros in a base 10 factorial number might be close to 1/5. However, though the factorial products add another terminating 0 every factor of 5 multiplied into the product, as the numbers become quite large, the number of digits in the factorial product expands exponentially, and so the number above the terminating zeros tends toward 10% of each digit from 0 to 1 as the factorial becomes larger. Thus, as the factorials become larger, the proportion of 0 digits in the factorial products shifts slowly from around 1/5 toward 1/10, since the number of terminating zeros in n! increases only in proportion to n, whereas the number of digits of n! in base 10 increases exponentially.

The task

Create a function to calculate the mean of the proportions of 0 digits out of the total digits found in each factorial product from 1! to N!. This proportion of 0 digits in base 10 should be calculated using the number as printed as a base 10 integer.

Example: for 1 to 6 we have 1!, 2!, 3!, 4!, 5!, 6!, or (1, 2, 6, 24, 120, 720), so we need the mean of (0/1, 0/1, 0/1, 0/2, 1/3, 1/3) = (2/3) (totals of each proportion) / 6 (= N), or 0.1111111...

Example: for 1 to 25 the mean of the proportions of 0 digits in the factorial products series of N! with N from 1 to 25 is 0.26787.

Do this task for 1 to N where N is in (100, 1000, and 10000), so, compute the mean of the proportion of 0 digits for each product in the series of each of the factorials from 1 to 100, 1 to 1000, and 1 to 10000.

Stretch task

Find the N in 10000 < N < 50000 where the mean of the proportions of 0 digits in the factorial products from 1 to N permanently falls below 0.16. This task took many hours in the Python example, though I wonder if there is a faster algorithm out there.

11l

Translation of: Python

F facpropzeros(n, verbose = 1B)
   V proportions = [0.0] * n
   V (fac, psum) = (BigInt(1), 0.0)
   L(i) 0 .< n
      fac *= i + 1
      V d = String(fac)
      psum += sum(d.map(x -> Int(x == ‘0’))) / Float(d.len)
      proportions[i] = psum / (i + 1)

   I verbose
      print(‘The mean proportion of 0 in factorials from 1 to #. is #..’.format(n, psum / n))

   R proportions

L(n) [100, 1000, 10000]
   facpropzeros(n)

Output:

The mean proportion of 0 in factorials from 1 to 100 is 0.246753186.
The mean proportion of 0 in factorials from 1 to 1000 is 0.203544551.
The mean proportion of 0 in factorials from 1 to 10000 is 0.173003848.

Base 1000 version

F zinit()
   V zc = [0] * 999
   L(x) 1..9
      zc[x - 1] = 2
      zc[10 * x - 1] = 2
      zc[100 * x - 1] = 2
      L(y) (10.<100).step(10)
         zc[y + x - 1] = 1
         zc[10 * y + x - 1] = 1
         zc[10 * (y + x) - 1] = 1

   R zc

F meanfactorialdigits()
   V zc = zinit()
   V rfs = [1]
   V (total, trail, first) = (0.0, 1, 0)
   L(f) 2 .< 50000
      V (carry, d999, zeroes) = (0, 0, (trail - 1) * 3)
      V (j, l) = (trail, rfs.len)
      L j <= l | carry != 0
         I j <= l
            carry = rfs[j - 1] * f + carry

         d999 = carry % 1000
         I j <= l
            rfs[j - 1] = d999
         E
            rfs.append(d999)

         zeroes += I d999 == 0 {3} E zc[d999 - 1]
         carry I/= 1000
         j++

      L rfs[trail - 1] == 0
         trail++

      d999 = rfs.last
      d999 = I d999 >= 100 {0} E I d999 < 10 {2} E 1

      zeroes -= d999
      V digits = rfs.len * 3 - d999
      total += Float(zeroes) / digits
      V ratio = total / f
      I f C [100, 1000, 10000]
         print(‘The mean proportion of zero digits in factorials to #. is #.’.format(f, ratio))

      I ratio >= 0.16
         first = 0
      E I first == 0
         first = f

   print(‘The mean proportion dips permanently below 0.16 at ’first‘.’)

meanfactorialdigits()

Output:

The mean proportion of zero digits in factorials to 100 is 0.246753186
The mean proportion of zero digits in factorials to 1000 is 0.203544551
The mean proportion of zero digits in factorials to 10000 is 0.173003848
The mean proportion dips permanently below 0.16 at 47332.

Arturo

su: 0.0
f: 1
lim: 100
loop 1..10000 'n [
    'f * n
    str: to :string f
    'su + (enumerate str 'c -> c = `0`) // size str
    if n = lim [
        print [n ":" su // n]
        'lim * 10
    ]
]

Output:

100 : 0.2467531861674322 
1000 : 0.2035445511031646 
10000 : 0.1730038482418671

C#

Translation of: Java

using System;
using System.Collections.Generic;
using System.Numerics;

public class DistributionInFactorials
{
    public static void Main(string[] args)
    {
        List<int> limits = new List<int> { 100, 1_000, 10_000 };
        foreach (int limit in limits)
        {
            MeanFactorialDigits(limit);
        }
    }

    private static void MeanFactorialDigits(int limit)
    {
        BigInteger factorial = BigInteger.One;
        double proportionSum = 0.0;
        double proportionMean = 0.0;

        for (int n = 1; n <= limit; n++)
        {
            factorial = factorial * n;
            string factorialString = factorial.ToString();
            int digitCount = factorialString.Length;
            long zeroCount = factorialString.Split('0').Length - 1;
            proportionSum += (double)zeroCount / digitCount;
            proportionMean = proportionSum / n;
        }

        string result = string.Format("{0:F8}", proportionMean);
        Console.WriteLine("Mean proportion of zero digits in factorials from 1 to " + limit + " is " + result);
    }
}

Output:

Mean proportion of zero digits in factorials from 1 to 100 is 0.24675319
Mean proportion of zero digits in factorials from 1 to 1000 is 0.20354455
Mean proportion of zero digits in factorials from 1 to 10000 is 0.17300385

C++

Translation of: Phix

#include <array>
#include <chrono>
#include <iomanip>
#include <iostream>
#include <vector>

auto init_zc() {
    std::array<int, 1000> zc;
    zc.fill(0);
    zc[0] = 3;
    for (int x = 1; x <= 9; ++x) {
        zc[x] = 2;
        zc[10 * x] = 2;
        zc[100 * x] = 2;
        for (int y = 10; y <= 90; y += 10) {
            zc[y + x] = 1;
            zc[10 * y + x] = 1;
            zc[10 * (y + x)] = 1;
        }
    }
    return zc;
}

template <typename clock_type>
auto elapsed(const std::chrono::time_point<clock_type>& t0) {
    auto t1 = clock_type::now();
    auto duration =
        std::chrono::duration_cast<std::chrono::milliseconds>(t1 - t0);
    return duration.count();
}

int main() {
    auto zc = init_zc();
    auto t0 = std::chrono::high_resolution_clock::now();
    int trail = 1, first = 0;
    double total = 0;
    std::vector<int> rfs{1};
    std::cout << std::fixed << std::setprecision(10);
    for (int f = 2; f <= 50000; ++f) {
        int carry = 0, d999, zeroes = (trail - 1) * 3, len = rfs.size();
        for (int j = trail - 1; j < len || carry != 0; ++j) {
            if (j < len)
                carry += rfs[j] * f;
            d999 = carry % 1000;
            if (j < len)
                rfs[j] = d999;
            else
                rfs.push_back(d999);
            zeroes += zc[d999];
            carry /= 1000;
        }
        while (rfs[trail - 1] == 0)
            ++trail;
        d999 = rfs.back();
        d999 = d999 < 100 ? (d999 < 10 ? 2 : 1) : 0;
        zeroes -= d999;
        int digits = rfs.size() * 3 - d999;
        total += double(zeroes) / digits;
        double ratio = total / f;
        if (ratio >= 0.16)
            first = 0;
        else if (first == 0)
            first = f;
        if (f == 100 || f == 1000 || f == 10000) {
            std::cout << "Mean proportion of zero digits in factorials to " << f
                      << " is " << ratio << ". (" << elapsed(t0) << "ms)\n";
        }
    }
    std::cout << "The mean proportion dips permanently below 0.16 at " << first
              << ". (" << elapsed(t0) << "ms)\n";
}

Output:

Mean proportion of zero digits in factorials to 100 is 0.2467531862. (0ms)
Mean proportion of zero digits in factorials to 1000 is 0.2035445511. (1ms)
Mean proportion of zero digits in factorials to 10000 is 0.1730038482. (152ms)
The mean proportion dips permanently below 0.16 at 47332. (4598ms)

Go

Brute force

Library: GMP(Go wrapper)

Library: Go-rcu

Timings here are 2.8 seconds for the basic task and 182.5 seconds for the stretch goal.

package main

import (
    "fmt"
    big "github.com/ncw/gmp"
    "rcu"
)

func main() {
    fact  := big.NewInt(1)
    sum   := 0.0
    first := int64(0)
    firstRatio := 0.0    
    fmt.Println("The mean proportion of zero digits in factorials up to the following are:")
    for n := int64(1); n <= 50000; n++  {
        fact.Mul(fact, big.NewInt(n))
        bytes  := []byte(fact.String())
        digits := len(bytes)
        zeros  := 0
        for _, b := range bytes {
            if b == '0' {
                zeros++
            }
        }
        sum += float64(zeros)/float64(digits)
        ratio := sum / float64(n)
        if n == 100 || n == 1000 || n == 10000 {
            fmt.Printf("%6s = %12.10f\n", rcu.Commatize(int(n)), ratio)
        } 
        if first > 0 && ratio >= 0.16 {
            first = 0
            firstRatio = 0.0
        } else if first == 0 && ratio < 0.16 {
            first = n
            firstRatio = ratio           
        }
    }
    fmt.Printf("%6s = %12.10f", rcu.Commatize(int(first)), firstRatio)
    fmt.Println(" (stays below 0.16 after this)")
    fmt.Printf("%6s = %12.10f\n", "50,000", sum / 50000)
}

Output:

The mean proportion of zero digits in factorials up to the following are:
   100 = 0.2467531862
 1,000 = 0.2035445511
10,000 = 0.1730038482
47,332 = 0.1599999958 (stays below 0.16 after this)
50,000 = 0.1596200546

'String math' and base 1000

Translation of: Phix

Much quicker than before with 10,000 now being reached in 0.35 seconds and the stretch goal in about 5.5 seconds.

package main

import (
    "fmt"
    "rcu"
)

var rfs = []int{1} // reverse factorial(1) in base 1000
var zc = make([]int, 999)

func init() {
    for x := 1; x <= 9; x++ {
        zc[x-1] = 2     // 00x
        zc[10*x-1] = 2  // 0x0
        zc[100*x-1] = 2 // x00
        var y = 10
        for y <= 90 {
            zc[y+x-1] = 1      // 0yx
            zc[10*y+x-1] = 1   // y0x
            zc[10*(y+x)-1] = 1 // yx0
            y += 10
        }
    }
}

func main() {
    total := 0.0
    trail := 1
    first := 0
    firstRatio := 0.0
    fmt.Println("The mean proportion of zero digits in factorials up to the following are:")
    for f := 2; f <= 10000; f++ {
        carry := 0
        d999 := 0
        zeros := (trail - 1) * 3
        j := trail
        l := len(rfs)
        for j <= l || carry != 0 {
            if j <= l {
                carry = rfs[j-1]*f + carry
            }
            d999 = carry % 1000
            if j <= l {
                rfs[j-1] = d999
            } else {
                rfs = append(rfs, d999)
            }
            if d999 == 0 {
                zeros += 3
            } else {
                zeros += zc[d999-1]
            }
            carry /= 1000
            j++
        }
        for rfs[trail-1] == 0 {
            trail++
        }
        // d999 = quick correction for length and zeros
        d999 = rfs[len(rfs)-1]
        if d999 < 100 {
            if d999 < 10 {
                d999 = 2
            } else {
                d999 = 1
            }
        } else {
            d999 = 0
        }
        zeros -= d999
        digits := len(rfs)*3 - d999
        total += float64(zeros) / float64(digits)
        ratio := total / float64(f)
        if ratio >= 0.16 {
            first = 0
            firstRatio = 0.0
        } else if first == 0 {
            first = f
            firstRatio = ratio
        }
        if f == 100 || f == 1000 || f == 10000 {
            fmt.Printf("%6s = %12.10f\n", rcu.Commatize(f), ratio)
        }
    }
    fmt.Printf("%6s = %12.10f", rcu.Commatize(first), firstRatio)
    fmt.Println(" (stays below 0.16 after this)")
    fmt.Printf("%6s = %12.10f\n", "50,000", total/50000)
}

Output:

Same as 'brute force' version.

Java

import java.math.BigInteger;
import java.util.List;

public final class DistributionInFactorials {

	public static void main(String[] aArgs) {
		List<Integer> limits = List.of( 100, 1_000, 10_000 );
		for ( Integer limit : limits ) {
			meanFactorialDigits(limit);
		}
	}
	
	private static void meanFactorialDigits(Integer aLimit) {
		BigInteger factorial = BigInteger.ONE;
		double proportionSum = 0.0;
		double proportionMean = 0.0;
		
		for ( int n = 1; n <= aLimit; n++ ) {
	        factorial = factorial.multiply(BigInteger.valueOf(n));
	        String factorialString = factorial.toString();
	        int digitCount = factorialString.length();	        
	        long zeroCount = factorialString.chars().filter( ch -> ch == '0' ).count();      
	        proportionSum += (double) zeroCount / digitCount;
	        proportionMean = proportionSum / n;	        
		}
		
		String result = String.format("%.8f", proportionMean);
		System.out.println("Mean proportion of zero digits in factorials from 1 to " + aLimit + " is " + result);
	}	

}

Output:

Mean proportion of zero digits in factorials from 1 to 100 is 0.24675319
Mean proportion of zero digits in factorials from 1 to 1000 is 0.20354455
Mean proportion of zero digits in factorials from 1 to 10000 is 0.17300385

jq

Works with jq

The precision of jq's integer arithmetic is not up to this task, so in the following we borrow from the "BigInt" library and use a string representation of integers.

Unfortunately, although gojq (the Go implementation of jq) does support unbounded-precision integer arithmetic, it is unsuited for the task because of memory management issues.

From BigInt.jq

# multiply two decimal strings, which may be signed (+ or -)
def long_multiply(num1; num2):

  def stripsign:
    .[0:1] as $a
    | if $a == "-" then [ -1, .[1:]] 
    elif $a == "+" then [  1, .[1:]] 
    else [1, .]
    end;

  def adjustsign(sign):
     if sign == 1 then . else "-" + . end;

  # mult/2 assumes neither argument has a sign
  def mult(num1;num2):
      (num1 | explode | map(.-48) | reverse) as $a1
    | (num2 | explode | map(.-48) | reverse) as $a2
    | reduce range(0; num1|length) as $i1
        ([];  # result
         reduce range(0; num2|length) as $i2
           (.;
            ($i1 + $i2) as $ix
            | ( $a1[$i1] * $a2[$i2] + (if $ix >= length then 0 else .[$ix] end) ) as $r
            | if $r > 9 # carrying
              then
                .[$ix + 1] = ($r / 10 | floor) +  (if $ix + 1 >= length then 0 else .[$ix + 1] end )
                | .[$ix] = $r - ( $r / 10 | floor ) * 10
              else
                .[$ix] = $r
              end
         )
        ) 
    | reverse | map(.+48) | implode;

  (num1|stripsign) as $a1
  | (num2|stripsign) as $a2
  | if $a1[1] == "0" or  $a2[1] == "0" then "0"
    elif $a1[1] == "1" then $a2[1]|adjustsign( $a1[0] * $a2[0] )
    elif $a2[1] == "1" then $a1[1]|adjustsign( $a1[0] * $a2[0] )
    else mult($a1[1]; $a2[1]) | adjustsign( $a1[0] * $a2[0] )
    end;

The task

def count(s): reduce s as $x (0; .+1);

def meanfactorialdigits:
   def digits: tostring | explode;
   def nzeros: count( .[] | select(. == 48) ); # "0" is 48
   
   . as $N
   | 0.16 as $goal
   | label $out
   | reduce range( 1; 1+$N ) as $i ( {factorial: "1", proportionsum: 0.0, first: null };
        .factorial = long_multiply(.factorial; $i|tostring)
        | (.factorial|digits) as $d
        | .proportionsum += ($d | (nzeros / length)) 
        | (.proportionsum / $i) as $propmean
	| if .first
	  then if $propmean > $goal then .first = null else . end
	  elif $propmean <= $goal then .first = $i
	  else .
	  end)
    | "Mean proportion of zero digits in factorials to \($N) is \(.proportionsum/$N);" +
       (if .first then " mean <= \($goal) from N=\(.first) on." else " goal (\($goal)) unmet." end);

# The task:
100, 1000, 10000 | meanfactorialdigits

Output:

Mean proportion of zero digits in factorials to 100 is 0.24675318616743216; goal (0.16) unmet.
Mean proportion of zero digits in factorials to 1000 is 0.20354455110316458; goal (0.16) unmet.
Mean proportion of zero digits in factorials to 10000 is 0.17300384824186707; goal (0.16) unmet.

Julia

function meanfactorialdigits(N, goal = 0.0)
    factoril, proportionsum = big"1", 0.0
    for i in 1:N
        factoril *= i
        d = digits(factoril)
        zero_proportion_in_fac = count(x -> x == 0, d) / length(d)
        proportionsum += zero_proportion_in_fac
        propmean = proportionsum / i
        if i > 15 && propmean <= goal
            println("The mean proportion dips permanently below $goal at $i.")
            break
        end
        if i == N
            println("Mean proportion of zero digits in factorials to $N is ", propmean)
        end
    end
end

@time foreach(meanfactorialdigits, [100, 1000, 10000])

@time meanfactorialdigits(50000, 0.16)

Output:

Mean proportion of zero digits in factorials to 100 is 0.24675318616743216
Mean proportion of zero digits in factorials to 1000 is 0.20354455110316458
Mean proportion of zero digits in factorials to 10000 is 0.17300384824186707
  3.030182 seconds (297.84 k allocations: 1.669 GiB, 0.83% gc time, 0.28% compilation time)
The mean proportion dips permanently below 0.16 at 47332.
179.157788 seconds (3.65 M allocations: 59.696 GiB, 1.11% gc time)

Base 1000 version

Translation of: Pascal, Phix

function init_zc()
    zc = zeros(Int, 999)
    for x in 1:9
        zc[x] = 2       # 00x
        zc[10*x] = 2    # 0x0
        zc[100*x] = 2   # x00
        for y in 10:10:90
            zc[y+x] = 1         # 0yx
            zc[10*y+x] = 1      # y0x
            zc[10*(y+x)] = 1    # yx0
        end
    end
    return zc
end

function meanfactorialzeros(N = 50000, verbose = true)
    zc = init_zc()
    rfs = [1]

    total, trail, first, firstratio = 0.0, 1, 0, 0.0

    for f in 2:N
        carry, d999, zeroes = 0, 0, (trail - 1) * 3
        j, l = trail, length(rfs)
        while j <= l || carry != 0
            if j <= l
                carry = (rfs[j]) * f + carry
            end
            d999 = carry % 1000
            if j <= l
                rfs[j] = d999
            else
                push!(rfs, d999)
            end
            zeroes += (d999 == 0) ? 3 : zc[d999]
            carry ÷= 1000
            j += 1
        end
        while rfs[trail] == 0
            trail += 1
        end
        # d999 = quick correction for length and zeroes:
        d999 = rfs[end]
        d999 = d999 < 100 ? d999 < 10 ? 2 : 1 : 0
        zeroes -= d999
        digits = length(rfs) * 3 - d999
        total += zeroes / digits
        ratio = total / f
        if ratio >= 0.16
           first = 0
           firstratio = 0.0
        elseif first == 0
            first = f
            firstratio = ratio
        end
        if f in [100, 1000, 10000]
            verbose && println("Mean proportion of zero digits in factorials to $f is $ratio")
        end
    end
    verbose && println("The mean proportion dips permanently below 0.16 at $first.")
end

meanfactorialzeros(100, false)
@time meanfactorialzeros()

Output:

Mean proportion of zero digits in factorials to 100 is 0.24675318616743216
Mean proportion of zero digits in factorials to 1000 is 0.20354455110316458
Mean proportion of zero digits in factorials to 10000 is 0.17300384824186707
The mean proportion dips permanently below 0.16 at 47332.
  4.638323 seconds (50.08 k allocations: 7.352 MiB)

Mathematica/Wolfram Language

ClearAll[ZeroDigitsFractionFactorial]
ZeroDigitsFractionFactorial[n_Integer] := Module[{m},
  m = IntegerDigits[n!];
  Count[m, 0]/Length[m]
  ]
ZeroDigitsFractionFactorial /@ Range[6] // Mean // N
ZeroDigitsFractionFactorial /@ Range[25] // Mean // N
ZeroDigitsFractionFactorial /@ Range[100] // Mean // N
ZeroDigitsFractionFactorial /@ Range[1000] // Mean // N
ZeroDigitsFractionFactorial /@ Range[10000] // Mean // N

fracs = ParallelMap[ZeroDigitsFractionFactorial, Range[50000], Method -> ("ItemsPerEvaluation" -> 100)];
means = Accumulate[N@fracs]/Range[Length[fracs]];
len = LengthWhile[Reverse@means, # < 0.16 &];
50000 - len + 1

Output:

0.111111
0.267873
0.246753
0.203545
0.173004
47332

Nim

Task

Library: bignum

import strutils, std/monotimes
import bignum

let t0 = getMonoTime()
var sum = 0.0
var f = newInt(1)
var lim = 100
for n in 1..10_000:
  f *= n
  let str = $f
  sum += str.count('0') / str.len
  if n == lim:
    echo n, ":\t", sum / float(n)
    lim *= 10
echo()
echo getMonoTime() - t0

Output:

100:    0.2467531861674322
1000:   0.2035445511031646
10000:  0.1730038482418671

(seconds: 2, nanosecond: 857794404)

Stretch task

Library: bignum

At each step, we eliminate the trailing zeroes to reduce the length of the number and save some time. But this is not much, about 8%.

import strutils, std/monotimes
import bignum

let t0 = getMonoTime()
var sum = 0.0
var first = 0
var f = newInt(1)
var count0 = 0
for n in 1..<50_000:
  f *= n
  while f mod 10 == 0:    # Reduce the length of "f".
    f = f div 10
    inc count0
  let str = $f
  sum += (str.count('0') + count0) / (str.len + count0)
  if sum / float(n) < 0.16:
    if first == 0: first = n
  else:
    first = 0

echo "Permanently below 0.16 at n = ", first
echo "Execution time: ", getMonoTime() - t0

Output:

Permanently below 0.16 at n = 47332
Execution time: (seconds: 190, nanosecond: 215845101)

Pascal

Doing the calculation in Base 1,000,000,000 like in Primorial_numbers#alternative.
The most time consuming is converting to string and search for zeros.
Therefor I do not convert to string.I divide the base in sections of 3 digits with counting zeros in a lookup table.

program Factorial;
{$IFDEF FPC} {$MODE DELPHI} {$Optimization ON,ALL} {$ENDIF}
uses
  sysutils;
type
  tMul = array of LongWord;
  tpMul = pLongWord;
const
  LongWordDec = 1000*1000*1000;
  LIMIT = 50000;
var
  CountOfZero : array[0..999] of byte;
  SumOfRatio :array[0..LIMIT] of extended;


procedure OutMul(pMul:tpMul;Lmt :NativeInt);
// for testing
Begin
  write(pMul[lmt]);
  For lmt := lmt-1  downto 0 do
    write(Format('%.9d',[pMul[lmt]]));
  writeln;
end;

procedure InitCoZ;
//Init Lookup table for 3 digits
var
  x,y : integer;
begin
  fillchar(CountOfZero,SizeOf(CountOfZero),#0);
  CountOfZero[0] := 3; //000
  For x := 1 to 9 do
  Begin
    CountOfZero[x] := 2;     //00x
    CountOfZero[10*x] := 2;  //0x0
    CountOfZero[100*x] := 2; //x00
    y := 10;
    repeat
      CountOfZero[y+x] := 1;      //0yx
      CountOfZero[10*y+x] := 1;   //y0x
      CountOfZero[10*(y+x)] := 1; //yx0
      inc(y,10)
    until y > 100;
  end;
end;

function getFactorialDecDigits(n:NativeInt):NativeInt;
var
  res: extended;
Begin
  result := -1;
  IF (n > 0) AND (n <= 1000*1000) then
  Begin
    res := 0;
    repeat res := res+ln(n); dec(n); until n < 2;
    result := trunc(res/ln(10))+1;
  end;
end;

function CntZero(pMul:tpMul;Lmt :NativeInt):NativeUint;
//count zeros in Base 1,000,000,000 number
var
  q,r : LongWord;
  i : NativeInt;
begin
  result := 0;
  For i := Lmt-1 downto 0 do
  Begin
    q := pMul[i];
    r := q DIV 1000;
    result +=CountOfZero[q-1000*r];//q-1000*r == q mod 1000
    q := r;
    r := q DIV 1000;
    result +=CountOfZero[q-1000*r];
    q := r;
    r := q DIV 1000;
    result +=CountOfZero[q-1000*r];
  end;
//special case first digits no leading '0'
  q := pMul[lmt];
  while q >= 1000 do
  begin
    r := q DIV 1000;
    result +=CountOfZero[q-1000*r];
    q := r;
  end;
  while q > 0 do
  begin
    r := q DIV 10;
    result += Ord( q-10*r= 0);
    q := r;
  end;
end;

function GetCoD(pMul:tpMul;Lmt :NativeInt):NativeUint;
//count of decimal digits
var
  i : longWord;
begin
  result := 9*Lmt;
  i := pMul[Lmt];
  while i > 1000 do
  begin
    i := i DIV 1000;
    inc(result,3);
  end;
  while i > 0 do
  begin
    i := i DIV 10;
    inc(result);
  end;
end;

procedure DoChecks(pMul:tpMul;Lmt,i :NativeInt);
//(extended(1.0)* makes TIO.RUN faster // only using FPU?
Begin
  SumOfRatio[i] := SumOfRatio[i-1] + (extended(1.0)*CntZero(pMul,Lmt))/GetCoD(pMul,Lmt);
end;

function MulByI(pMul:tpMul;UL,i :NativeInt):NativeInt;
var
  prod  : Uint64;
  j     : nativeInt;
  carry : LongWord;
begin
  result := UL;
  carry := 0;
  For j := 0 to result do
  Begin
    prod  := i*pMul[0]+Carry;
    Carry := prod Div LongWordDec;
    pMul[0] := Prod - LongWordDec*Carry;
    inc(pMul);
  end;

  IF Carry <> 0 then
  Begin
    inc(result);
    pMul[0]:= Carry;
  End;
end;

procedure getFactorialExact(n:NativeInt);
var
  MulArr : tMul;
  pMul : tpMul;
  i,ul : NativeInt;
begin
  i := getFactorialDecDigits(n) DIV 9 +10;
  Setlength(MulArr,i);
  pMul := @MulArr[0];
  Ul := 0;
  pMul[Ul]:= 1;
  i := 1;
  repeat
    UL := MulByI(pMul,UL,i);
    //Now do what you like to do with i!
    DoChecks(pMul,UL,i);
    inc(i);
  until i> n;
end;

procedure Out_(i: integer);
begin
  if i > LIMIT then
    EXIT;
  writeln(i:8,SumOfRatio[i]/i:18:15);
end;

var
  i : integer;
Begin
  InitCoZ;
  SumOfRatio[0]:= 0;
  getFactorialExact(LIMIT);
  Out_(100);
  Out_(1000);
  Out_(10000);
  Out_(50000);
  i := limit;
  while i >0 do
  Begin
    if SumOfRatio[i]/i >0.16 then
      break;
    dec(i);
  end;
  inc(i);
  writeln('First ratio < 0.16 ', i:8,SumOfRatio[i]/i:20:17);
end.

Output:

     100 0.246753186167432
    1000 0.203544551103165
   10000 0.173003848241866
   50000 0.159620054602269
First ratio < 0.16    47332 0.15999999579985665 
Real time: 4.898 s  CPU share: 99.55 % // 2.67s on 2200G freepascal 3.2.2

Perl

Library: ntheory

use strict;
use warnings;
use ntheory qw/factorial/;

for my $n (100, 1000, 10000) {
    my($sum,$f) = 0;
    $f = factorial $_ and $sum += ($f =~ tr/0//) / length $f for 1..$n;
    printf "%5d: %.5f\n", $n, $sum/$n;
}

Output:

  100: 0.24675
 1000: 0.20354
10000: 0.17300

Phix

Using "string math" to create reversed factorials, for slightly easier skipping of "trailing" zeroes, but converted to base 1000 and with the zero counting idea from Pascal, which sped it up threefold.

with javascript_semantics
sequence rfs = {1}  -- reverse factorial(1) in base 1000
         
function init_zc()
    sequence zc = repeat(0,999)
    for x=1 to 9 do
        zc[x] = 2       -- 00x
        zc[10*x] = 2    -- 0x0
        zc[100*x] = 2   -- x00
        for y=10 to 90 by 10 do
            zc[y+x] = 1         -- 0yx
            zc[10*y+x] = 1      -- y0x
            zc[10*(y+x)] = 1    -- yx0
        end for
    end for
    return zc
end function
constant zc = init_zc()

atom t0 = time(),
     total = 0
integer trail = 1,
        first = 0
for f=2 to iff(platform()=JS?10000:50000) do
    integer carry = 0, d999, 
            zeroes = (trail-1)*3, 
            j = trail, l = length(rfs)
    while j<=l or carry do
        if j<=l then
            carry = (rfs[j])*f+carry
        end if
        d999 = remainder(carry,1000)
        if j<=l then
            rfs[j] = d999
        else
            rfs &= d999
        end if
        zeroes += iff(d999=0?3:zc[d999])
        carry = floor(carry/1000)
        j += 1
    end while
    while rfs[trail]=0 do trail += 1 end while
    -- d999 := quick correction for length and zeroes:
    d999 = rfs[$]
    d999 = iff(d999<100?iff(d999<10?2:1):0)
    zeroes -= d999
    integer digits = length(rfs)*3-d999

    total += zeroes/digits
    atom ratio = total/f
    if ratio>=0.16 then
        first = 0
    elsif first=0 then
        first = f
    end if
    if find(f,{100,1000,10000}) then
        string e = elapsed(time()-t0)
        printf(1,"Mean proportion of zero digits in factorials to %d is %.10f (%s)\n",{f,ratio,e})
    end if
end for
if platform()!=JS then
    string e = elapsed(time()-t0)
    printf(1,"The mean proportion dips permanently below 0.16 at %d. (%s)\n",{first,e})
end if

Output:

Mean proportion of zero digits in factorials to 100 is 0.2467531862 (0s)
Mean proportion of zero digits in factorials to 1000 is 0.2035445511 (0.2s)
Mean proportion of zero digits in factorials to 10000 is 0.1730038482 (2.3s)
The mean proportion dips permanently below 0.16 at 47332. (1 minute and 2s)

(stretch goal removed under pwa/p2js since otherwise you'd get a blank screen for 2 or 3 minutes)

trailing zeroes only

Should you only be interested in the ratio of trailing zeroes, you can do that much faster:

with javascript_semantics
atom t0 = time(),
     f10 = log10(1),
     total = 0
integer first = 0
for f=2 to 50000 do
    f10 += log10(f)
    integer digits = ceil(f10),
            zeroes = 0,
            v = 5
    while v<=f do
        zeroes += floor(f/v)
        v *= 5
    end while
    total += zeroes/digits
    atom ratio = total/f
    if ratio>=0.07 then
        first = 0
    elsif first=0 then
        first = f
    end if
    if find(f,{100,1000,10000}) then
        printf(1,"Mean proportion of trailing zeroes in factorials to %d is %f\n",{f,ratio})
    end if
end for
string e = elapsed(time()-t0)
printf(1,"The mean proportion dips permanently below 0.07 at %d. (%s)\n",{first,e})

Output:

Mean proportion of trailing zeroes in factorials to 100 is 0.170338
Mean proportion of trailing zeroes in factorials to 1000 is 0.116334
Mean proportion of trailing zeroes in factorials to 10000 is 0.081267
The mean proportion dips permanently below 0.07 at 31549. (0.1s)

Python

def facpropzeros(N, verbose = True):
    proportions = [0.0] * N
    fac, psum = 1, 0.0
    for i in range(N):
        fac *= i + 1
        d = list(str(fac))
        psum += sum(map(lambda x: x == '0', d)) / len(d)
        proportions[i] = psum / (i + 1)

    if verbose:
        print("The mean proportion of 0 in factorials from 1 to {} is {}.".format(N, psum / N))

    return proportions


for n in [100, 1000, 10000]:
    facpropzeros(n)

props = facpropzeros(47500, False)
n = (next(i for i in reversed(range(len(props))) if props[i] > 0.16))

print("The mean proportion dips permanently below 0.16 at {}.".format(n + 2))

Output:

The mean proportion of 0 in factorials from 1 to 100 is 0.24675318616743216.
The mean proportion of 0 in factorials from 1 to 1000 is 0.20354455110316458.
The mean proportion of 0 in factorials from 1 to 10000 is 0.17300384824186707.
The mean proportion dips permanently below 0.16 at 47332.

The means can be plotted, showing a jump from 0 to over 0.25, followed by a slowly dropping curve:

import matplotlib.pyplot as plt
plt.plot([i+1 for i in range(len(props))], props)

Base 1000 version

Translation of: Go via Phix via Pascal

def zinit():
    zc = [0] * 999
    for x in range(1, 10):
        zc[x - 1] = 2        # 00x
        zc[10 * x - 1] = 2   # 0x0
        zc[100 * x - 1] = 2  # x00
        for y in range(10, 100, 10):
            zc[y + x - 1] = 1           # 0yx
            zc[10 * y + x - 1] = 1      # y0x
            zc[10 * (y + x) - 1] = 1    # yx0

    return zc

def meanfactorialdigits():
    zc = zinit()
    rfs = [1]
    total, trail, first = 0.0, 1, 0
    for f in range(2, 50000):
        carry, d999, zeroes = 0, 0, (trail - 1) * 3
        j, l = trail, len(rfs)
        while j <= l or carry != 0:
            if j <= l:
                carry = rfs[j-1] * f + carry

            d999 = carry % 1000
            if j <= l:
                rfs[j-1] = d999
            else:
                rfs.append(d999)

            zeroes += 3 if d999 == 0 else zc[d999-1]
            carry //= 1000
            j += 1

        while rfs[trail-1] == 0:
            trail += 1

        # d999 is a quick correction for length and zeros
        d999 = rfs[-1]
        d999 = 0 if d999 >= 100 else 2 if d999 < 10 else 1

        zeroes -= d999
        digits = len(rfs) * 3 - d999
        total += zeroes / digits
        ratio = total / f
        if f in [100, 1000, 10000]:
            print("The mean proportion of zero digits in factorials to {} is {}".format(f, ratio))
            
        if ratio >= 0.16:
            first = 0
        elif first == 0:
            first = f

    print("The mean proportion dips permanently below 0.16 at {}.".format(first))



import time
TIME0 = time.perf_counter()
meanfactorialdigits()
print("\nTotal time:", time.perf_counter() - TIME0, "seconds.")

Output:

The mean proportion of zero digits in factorials to 100 is 0.24675318616743216
The mean proportion of zero digits in factorials to 1000 is 0.20354455110316458
The mean proportion of zero digits in factorials to 10000 is 0.17300384824186707
The mean proportion dips permanently below 0.16 at 47332.

Total time: 648.3583232999999 seconds.

Raku

Works, but depressingly slow for 10000.

sub postfix:<!> (Int $n) { ( constant factorial = 1, 1, |[\*] 2..* )[$n] }
sink 10000!; # prime the iterator to allow multithreading

sub zs ($n) { ( constant zero-share = (^Inf).race(:32batch).map: { (.!.comb.Bag){'0'} / .!.chars } )[$n+1] }

.say for (
     100
    ,1000
    ,10000
).map:  -> \n { "{n}: {([+] (^n).map: *.&zs) / n}" }

Output:

100: 0.24675318616743216
1000: 0.20354455110316458
10000: 0.17300384824186605

REXX

/*REXX program computes the mean of the proportion of "0" digits a series of factorials.*/
parse arg $                                      /*obtain optional arguments from the CL*/
if $='' | $=","  then $= 100 1000 10000          /*not specified?  Then use the default.*/
#= words($)                                      /*the number of ranges to be used here.*/
numeric digits 100                               /*increase dec. digs, but only to 100. */
big= word($, #);  != 1                           /*obtain the largest number in ranges. */
                                do i=1  for big  /*calculate biggest  !  using 100 digs.*/
                                != ! * i         /*calculate the factorial of  BIG.     */
                                end   /*i*/
if pos('E', !)>0  then do                        /*In exponential format?  Then get EXP.*/
                       parse var !  'E'  x       /*parse the exponent from the number.  */
                       numeric digits    x+1     /*set the decimal digits to  X  plus 1.*/
                       end                       /* [↑]  the  +1  is for the dec. point.*/

title= ' mean proportion of zeros in the (decimal) factorial products for  N'
say '     N     │'center(title, 80)              /*display the title for the output.    */
say '───────────┼'center(""   , 80, '─')         /*   "     a   sep   "   "     "       */

  do j=1  for #;  n= word($, j)                  /*calculate some factorial ranges.     */
  say center( commas(n), 11)'│' left(0dist(n), 75)...    /*show results for above range.*/
  end   /*j*/

say '───────────┴'center(""   , 80, '─')         /*display a foot sep for the output.   */
exit 0                                           /*stick a fork in it,  we're all done. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
commas: parse arg ?;  do jc=length(?)-3  to 1  by -3; ?=insert(',', ?, jc); end;  return ?
/*──────────────────────────────────────────────────────────────────────────────────────*/
0dist:  procedure; parse arg z;        != 1;         y= 0
                     do k=1  for z;    != ! * k;     y= y   +   countstr(0, !) / length(!)
                     end   /*k*/
        return y/z

output   when using the default inputs:

     N     │       mean proportion of zeros in the (decimal) factorial products for  N
───────────┼────────────────────────────────────────────────────────────────────────────────
    100    │ 0.2467531861674322177784158871973526991129407033266153063813195937196095976...
   1,000   │ 0.2035445511031646356400438031711455302985741167890402203486699704599684047...
  10,000   │ 0.1730038482418660531800366428930706156810278809057883361518852958446868172...
───────────┴────────────────────────────────────────────────────────────────────────────────

Rust

Translation of: Phix

fn init_zc() -> Vec<usize> {
    let mut zc = vec![0; 1000];
    zc[0] = 3;
    for x in 1..=9 {
        zc[x] = 2;
        zc[10 * x] = 2;
        zc[100 * x] = 2;
        let mut y = 10;
        while y <= 90 {
            zc[y + x] = 1;
            zc[10 * y + x] = 1;
            zc[10 * (y + x)] = 1;
            y += 10;
        }
    }
    zc
}

fn main() {
    use std::time::Instant;
    let zc = init_zc();
    let t0 = Instant::now();
    let mut trail = 1;
    let mut first = 0;
    let mut total: f64 = 0.0;
    let mut rfs = vec![1];

    for f in 2..=50000 {
        let mut carry = 0;
        let mut d999: usize;
        let mut zeroes = (trail - 1) * 3;
        let len = rfs.len();
        let mut j = trail - 1;
        while j < len || carry != 0 {
            if j < len {
                carry += rfs[j] * f;
            }
            d999 = carry % 1000;
            if j < len {
                rfs[j] = d999;
            } else {
                rfs.push(d999);
            }
            zeroes += zc[d999];
            carry /= 1000;
            j += 1;
        }
        while rfs[trail - 1] == 0 {
            trail += 1;
        }
        d999 = rfs[rfs.len() - 1];
        d999 = if d999 < 100 {
            if d999 < 10 {
                2
            } else {
                1
            }
        } else {
            0
        };
        zeroes -= d999;
        let digits = rfs.len() * 3 - d999;
        total += (zeroes as f64) / (digits as f64);
        let ratio = total / (f as f64);
        if ratio >= 0.16 {
            first = 0;
        } else if first == 0 {
            first = f;
        }
        if f == 100 || f == 1000 || f == 10000 {
            let duration = t0.elapsed();
            println!(
                "Mean proportion of zero digits in factorials to {} is {:.10}. ({}ms)",
                f,
                ratio,
                duration.as_millis()
            );
        }
    }
    let duration = t0.elapsed();
    println!(
        "The mean proportion dips permanently below 0.16 at {}. ({}ms)",
        first,
        duration.as_millis()
    );
}

Output:

Mean proportion of zero digits in factorials to 100 is 0.2467531862. (0ms)
Mean proportion of zero digits in factorials to 1000 is 0.2035445511. (1ms)
Mean proportion of zero digits in factorials to 10000 is 0.1730038482. (149ms)
The mean proportion dips permanently below 0.16 at 47332. (4485ms)

Ruby

[100, 1000, 10_000].each do |n|
  v = 1
  total_proportion = (1..n).sum do |k|
    v *= k
    digits =  v.digits
    Rational(digits.count(0),  digits.size)
  end
  puts "The mean proportion of 0 in factorials from 1 to #{n} is #{(total_proportion/n).to_f}."
end

Output:

The mean proportion of 0 in factorials from 1 to 100 is 0.24675318616743222.
The mean proportion of 0 in factorials from 1 to 1000 is 0.20354455110316463.
The mean proportion of 0 in factorials from 1 to 10000 is 0.17300384824186604.

Sidef

func mean_factorial_digits(n, d = 0) {

    var v = 1
    var total = 0.float

    for k in (1..n) {
        v *= k
        total += v.digits.count(d)/v.len
    }

    total / n
}

say mean_factorial_digits(100)
say mean_factorial_digits(1000)
say mean_factorial_digits(10000)

Output:

0.246753186167432217778415887197352699112940703327
0.203544551103164635640043803171145530298574116789
0.173003848241866053180036642893070615681027880906

Wren

Brute force

Library: Wren-big

Library: Wren-fmt

Very slow indeed, 10.75 minutes to reach N = 10,000.

import "./big" for BigInt
import "./fmt" for Fmt

var fact = BigInt.one
var sum = 0
System.print("The mean proportion of zero digits in factorials up to the following are:")
for (n in 1..10000) {
    fact = fact * n
    var bytes = fact.toString.bytes
    var digits = bytes.count
    var zeros  = bytes.count { |b| b == 48 }
    sum = sum + zeros / digits
    if (n == 100 || n == 1000 || n == 10000) {
        Fmt.print("$,6d = $12.10f", n, sum / n)
    }
}

Output:

The mean proportion of zero digits in factorials up to the following are:
   100 = 0.2467531862
 1,000 = 0.2035445511
10,000 = 0.1730038482

'String math' and base 1000

Translation of: Phix

Around 60 times faster than before with 10,000 now being reached in about 10.5 seconds. Even the stretch goal is now viable and comes in at 5 minutes 41 seconds.

import "./fmt" for Fmt

var rfs = [1]  // reverse factorial(1) in base 1000

var init = Fn.new { |zc|
    for (x in 1..9) {
        zc[x-1] = 2         // 00x
        zc[10*x - 1] = 2    // 0x0
        zc[100*x - 1] = 2   // x00
        var y = 10
        while (y <= 90) {
            zc[y + x - 1] = 1       // 0yx
            zc[10*y + x - 1] = 1    // y0x
            zc[10*(y + x) - 1] = 1  // yx0
            y = y + 10
        }
    }
}

var zc = List.filled(999, 0)
init.call(zc)
var total = 0
var trail = 1
var first = 0
var firstRatio = 0
System.print("The mean proportion of zero digits in factorials up to the following are:")
for (f in 2..50000) {
    var carry = 0
    var d999 = 0
    var zeros = (trail-1) * 3
    var j = trail
    var l = rfs.count
    while (j <= l || carry != 0) {
        if (j <= l) carry = rfs[j-1]*f + carry
        d999 = carry % 1000
        if (j <= l) {
            rfs[j-1] = d999
        } else {
            rfs.add(d999)
        }
        zeros = zeros + ((d999 == 0) ? 3 : zc[d999-1])
        carry = (carry/1000).floor
        j = j + 1
    }
    while (rfs[trail-1] == 0) trail = trail + 1
    // d999 = quick correction for length and zeros
    d999 = rfs[-1]
    d999 = (d999 < 100) ? ((d999 < 10) ? 2 : 1) : 0
    zeros = zeros - d999
    var digits = rfs.count * 3 - d999
    total = total + zeros/digits
    var ratio =  total / f
    if (ratio >= 0.16) {
        first = 0
        firstRatio = 0
    } else if (first == 0) {
        first = f
        firstRatio = ratio
    }
    if (f == 100 || f == 1000 || f == 10000) {
        Fmt.print("$,6d = $12.10f", f, ratio)
    }
}
Fmt.write("$,6d = $12.10f", first, firstRatio)
System.print(" (stays below 0.16 after this)")
Fmt.print("$,6d = $12.10f", 50000, total/50000)

Output:

The mean proportion of zero digits in factorials up to the following are:
   100 = 0.2467531862
 1,000 = 0.2035445511
10,000 = 0.1730038482
47,332 = 0.1599999958 (stays below 0.16 after this)
50,000 = 0.1596200546