# of 4x4x4 centers correctly placed in a scramble

$ perl -le 'sub _ { my $out; for $n (0..5) { $out += grep {$ _ == $n} @{$_[$n]} } $out } for (1..1e6) { @centers = map {($_)x4} 0..5; my @cube; while (@centers) { push @{$cube[@centers % 6]}, splice @centers, (rand @centers), 1 } $c += _ @cube } print "$c / 1e6 = ", $c / 1e6'
3999087 / 1e6 = 3.999087

$ perl -mList::Util=max -le 'sub _ { my $out; for $n (0..5) { $out += grep {$ _ == $n} @{$_[$n]} } $out } for (1..1e5) { my @t; for (1..24) { @centers = map {($_)x4} 0..5; my @cube; while (@centers) { push @{$cube[@centers % 6]}, splice @centers, (rand @centers), 1 } push @t,  _ @cube } $c += max @t } print "$c / 1e5 = ", $c / 1e5'
792854 / 1e5 = 7.92854

$ perl -mList::Util=max,sum -le 'sub _ { my $out; for $n (0..5) { $out += grep {$ _ == $n} @{$_[$n]} } $out } for (1..1e7) { my @t; for (1..24) { @centers = map {($_)x4} 0..5; my @cube; while (@centers) { push @{$cube[@centers % 6]}, splice @centers, (rand @centers), 1 } push @t,  _ @cube; $ca[$t[-1]]++ } $cb[max @t]++ } ($a,$b) = ((sum @ca),sum @cb); print "\ncorrect centers | times | distribution\n\na)\ntotal: $a\nmean: ",  (sum map $_*$ca[$_], grep $ca[$_], 0..24) / $a, "\n"; $ca[$_] and printf "%2d: %8d  %.3f\n", $_, $ca[$_], $ca[$_]/$a for 0..24; print "\nb)\ntotal: $b\nmean: ", (sum map $_*$cb[$_], grep $cb[$_], 0..24) / $b, "\n"; $cb[$_] and printf "%2d: %8d  %.3f\n", $_, $cb[$_], $cb[$_]/$b for 0..24'

correct centers | times | distribution

a)
total: 240000000
mean: 4.00009665416667

 0:  3338073  0.014
 1: 15230287  0.063
 2: 33765699  0.141
 3: 48334329  0.201
 4: 50206280  0.209
 5: 40249724  0.168
 6: 25874279  0.108
 7: 13676951  0.057
 8:  6054183  0.025
 9:  2274718  0.009
10:   731267  0.003
11:   203462  0.001
12:    48524  0.000
13:    10035  0.000
14:     1850  0.000
15:      295  0.000
16:       40  0.000
17:        4  0.000

b)
total: 10000000
mean: 7.9294182

 4:      161  0.000
 5:    42235  0.004
 6:   848691  0.085
 7:  2971835  0.297
 8:  3331731  0.333
 9:  1856234  0.186
10:   688208  0.069
11:   200317  0.020
12:    48369  0.005
13:    10030  0.001
14:     1850  0.000
15:      295  0.000
16:       40  0.000
17:        4  0.000

<< DiscreteMath`Combinatorica`
a = Table[Floor[0.75 + n/4], {n, 1, 24}];
b = Table[0, {25}];
Do[(++b[[#]]) & /@ Table[Total[If[# == 0, 1,0] & /@(RandomPermutation[a] - a)] + 1, {1000}], {1000}]
(MapIndexed[(#*(#2[[1]] - 1)) &, b] // Total)/1000000 // N

a={139216, 635534, 1407486, 2013693, 2091757, 1674656, 1079544, 569305, 252506, 94666, 30534, 8533, 2076, 408, 66, 17, 2, 1, 0, 0, 0, 0, 0, 0, 0};
a=a/Total[a];
c=1;
x=(Table[d=c*((1-(Total[Drop[b, n]])/c)^24);c-=d;d,{n,Length[b]}]/.Indeterminate->0)// N;
(MapIndexed[(#*(#2[[1]]-1))&, x]//Total)/Total[x]//N