Distribution of 0 digits in factorial series: Difference between revisions
No edit summary |
(added pascal) |
||
Line 189: | Line 189: | ||
<pre>Permanently below 0.16 at n = 47332 |
<pre>Permanently below 0.16 at n = 47332 |
||
Execution time: (seconds: 190, nanosecond: 215845101)</pre> |
Execution time: (seconds: 190, nanosecond: 215845101)</pre> |
||
=={{header|Pascal}}== |
|||
Doing the calculation in Base 1,000,000,000 like in [[Primorial_numbers#alternative]].<BR> |
|||
The most time consuming is converting to string and search for zeros.<BR> |
|||
Therefor I do not convert to string.I divide the base in sections of 3 digits with counting zeros in a lookup table. |
|||
<lang pascal>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; |
|||
CoZ_Fac : array[0..LIMIT] of Uint32; |
|||
CoD_Fac : array[0..LIMIT] of Uint32; |
|||
procedure InitCoZ; |
|||
var |
|||
i,j : integer; |
|||
begin |
|||
fillchar(CountOfZero,SizeOf(CountOfZero),#0); |
|||
CountOfZero[0] := 3; |
|||
For i := 1 to 9 do |
|||
Begin |
|||
CountOfZero[i] := 2; |
|||
CountOfZero[10*i] := 2; |
|||
CountOfZero[100*i] := 2; |
|||
j := 10; |
|||
repeat |
|||
CountOfZero[j+i] := 1; |
|||
CountOfZero[10*j+i] := 1; |
|||
CountOfZero[10*(j+i)] := 1; |
|||
inc(j,10) |
|||
until j > 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; |
|||
procedure OutMul(pMul:tpMul;Lmt :NativeInt); |
|||
//for checking purposes |
|||
Begin |
|||
write(pMul[lmt]); |
|||
For lmt := lmt-1 downto 0 do |
|||
write(Format('%.9d',[pMul[lmt]])); |
|||
writeln; |
|||
end; |
|||
function CheckForZero(pMul:tpMul;Lmt :NativeInt):NativeUint; |
|||
// Numbers in base 1,000,000,000 divide in 3 sections |
|||
var |
|||
s : string[15]; |
|||
q,r : LongWord; |
|||
i : NativeInt; |
|||
begin |
|||
result := 0; |
|||
For i := 0 to Lmt-1 do |
|||
Begin |
|||
q := pMul[i]; |
|||
r := q DIV 1000; |
|||
result +=CountOfZero[q-1000*r]; |
|||
q := r; |
|||
r := q DIV 1000; |
|||
result +=CountOfZero[q-1000*r]; |
|||
q := r; |
|||
r := q DIV 1000; |
|||
result +=CountOfZero[q-1000*r]; |
|||
end; |
|||
q := pMul[lmt]; |
|||
while q >= 1000 do |
|||
begin |
|||
r := q DIV 1000; |
|||
result +=CountOfZero[q-1000*r]; |
|||
q := r; |
|||
end; |
|||
if q > 0 then |
|||
Begin |
|||
str(q,s); |
|||
For i := Length(s) downto 1 do |
|||
IF s[i] ='0' then |
|||
inc(result); |
|||
end; |
|||
end; |
|||
function GetCoD(pMul:tpMul;Lmt :NativeInt):NativeUint; |
|||
//calculate used 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 getFactorialExact(n:NativeInt); |
|||
var |
|||
MulArr : tMul; |
|||
pMul : tpMul; |
|||
prod,carry : Uint64; |
|||
i,j,ul : NativeInt; |
|||
begin |
|||
i := getFactorialDecDigits(n) DIV 9 +10; |
|||
Setlength(MulArr,i); |
|||
pMul := @MulArr[0]; |
|||
Ul := 0; |
|||
pMul[Ul]:= 1; |
|||
i := 1; |
|||
repeat |
|||
carry := 0; |
|||
For j := 0 to UL do |
|||
Begin |
|||
prod := i*pMul[j]+Carry; |
|||
Carry := prod Div LongWordDec; |
|||
pMul[j] := Prod - Carry*LongWordDec; |
|||
end; |
|||
IF Carry <> 0 then |
|||
Begin |
|||
inc(Ul); |
|||
pMul[UL]:= Carry; |
|||
End; |
|||
CoD_Fac[i] := GetCoD(pMul,UL); |
|||
CoZ_Fac[i] := CheckForZero(pMul,UL); |
|||
SumOfRatio[i] := SumOfRatio[i-1] + CoZ_Fac[i]/CoD_Fac[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; |
|||
CoD_Fac[0]:= 0; |
|||
CoZ_Fac[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.</lang> |
|||
{{out}} |
|||
<pre>TIO.RUN |
|||
100 0.246753186167432 |
|||
1000 0.203544551103165 |
|||
10000 0.173003848241866 |
|||
50000 0.159620054602269 |
|||
First ratio < 0.16 47332 0.15999999579985665 |
|||
Real time: 5.521 s CPU share: 98.03 % // 2.67s on 2200G freepascal 3.2.2</pre> |
|||
=={{header|Python}}== |
=={{header|Python}}== |
Revision as of 11:48, 12 June 2021
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.
Go
Brute force as I'll be surprised if there is a faster 'exact' algorithm for this task.
However, the combination of a native code compiler and GMP really cuts down the times (2.8 seconds for the basic task and 182.5 seconds for the stretch goal). Expect these times to be reduced further by the fastest languages. <lang go>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.Println(sum/ 50000)
}</lang>
- 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)
Julia
<lang 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)
</lang>
- 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)
Nim
Task
<lang Nim>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</lang>
- Output:
100: 0.2467531861674322 1000: 0.2035445511031646 10000: 0.1730038482418671 (seconds: 2, nanosecond: 857794404)
Stretch task
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%.
<lang Nim>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</lang>
- 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.
<lang pascal>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; CoZ_Fac : array[0..LIMIT] of Uint32; CoD_Fac : array[0..LIMIT] of Uint32;
procedure InitCoZ; var
i,j : integer;
begin
fillchar(CountOfZero,SizeOf(CountOfZero),#0); CountOfZero[0] := 3; For i := 1 to 9 do Begin CountOfZero[i] := 2; CountOfZero[10*i] := 2; CountOfZero[100*i] := 2; j := 10; repeat CountOfZero[j+i] := 1; CountOfZero[10*j+i] := 1; CountOfZero[10*(j+i)] := 1; inc(j,10) until j > 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;
procedure OutMul(pMul:tpMul;Lmt :NativeInt); //for checking purposes Begin
write(pMul[lmt]); For lmt := lmt-1 downto 0 do write(Format('%.9d',[pMul[lmt]])); writeln;
end;
function CheckForZero(pMul:tpMul;Lmt :NativeInt):NativeUint; // Numbers in base 1,000,000,000 divide in 3 sections var
s : string[15]; q,r : LongWord; i : NativeInt;
begin
result := 0; For i := 0 to Lmt-1 do Begin q := pMul[i]; r := q DIV 1000; result +=CountOfZero[q-1000*r]; q := r; r := q DIV 1000; result +=CountOfZero[q-1000*r]; q := r; r := q DIV 1000; result +=CountOfZero[q-1000*r]; end; q := pMul[lmt]; while q >= 1000 do begin r := q DIV 1000; result +=CountOfZero[q-1000*r]; q := r; end; if q > 0 then Begin str(q,s); For i := Length(s) downto 1 do IF s[i] ='0' then inc(result); end;
end;
function GetCoD(pMul:tpMul;Lmt :NativeInt):NativeUint; //calculate used 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 getFactorialExact(n:NativeInt); var
MulArr : tMul; pMul : tpMul; prod,carry : Uint64; i,j,ul : NativeInt;
begin
i := getFactorialDecDigits(n) DIV 9 +10; Setlength(MulArr,i); pMul := @MulArr[0]; Ul := 0; pMul[Ul]:= 1; i := 1; repeat carry := 0; For j := 0 to UL do Begin prod := i*pMul[j]+Carry; Carry := prod Div LongWordDec; pMul[j] := Prod - Carry*LongWordDec; end;
IF Carry <> 0 then Begin inc(Ul); pMul[UL]:= Carry; End;
CoD_Fac[i] := GetCoD(pMul,UL); CoZ_Fac[i] := CheckForZero(pMul,UL); SumOfRatio[i] := SumOfRatio[i-1] + CoZ_Fac[i]/CoD_Fac[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; CoD_Fac[0]:= 0; CoZ_Fac[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.</lang>
- Output:
TIO.RUN 100 0.246753186167432 1000 0.203544551103165 10000 0.173003848241866 50000 0.159620054602269 First ratio < 0.16 47332 0.15999999579985665 Real time: 5.521 s CPU share: 98.03 % // 2.67s on 2200G freepascal 3.2.2
Python
<lang 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))
</lang>
- 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: <lang python>import matplotlib.pyplot as plt plt.plot([i+1 for i in range(len(props))], props) </lang>
Raku
Works, but depressingly slow for 10000.
<lang perl6>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] }
.say for (
100 ,1000 ,10000
).map: -> \n { "{n}: {([+] (^n).map: *.&zs) / n}" }</lang>
- Output:
100: 0.24485445199021696 1000: 0.20336075048011162 10000: 0.17298757510670162
REXX
<lang 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 /*if its in exponential format, 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</lang>
- 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... ───────────┴────────────────────────────────────────────────────────────────────────────────
Wren
Very slow indeed, 10.75 minutes to reach N = 10,000. <lang ecmascript>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) }
}</lang>
- 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