Leverage FunctionCompile to dramatically accelerate the algorithm for generating terms of A383732 up to 1 billion digits of Pi, where a(n) denotes the smallest k such that each digit from 0 to 9 appears at least n times within the first k decimal digits of Pi.

Explain OEIS A383732

The OEIS sequence goes like
{
a
1
,
a
2
,
a
3
,
a
4
,
a
5
,...}={32,50,54,65,71,77,...}
. Here
a
1
=32
is the least number that each of
{0,1,2,3...,9}
appears at least once among the first 32 digits of π:
In[]:=
Count[#][piDigits[[;;31]]]&/@Range[0,9,1](*zeroismissinginthefirst31digits*)
Out[]=
{0,2,4,6,3,4,3,2,3,4}
Count[#][piDigits[[;;32]]]&/@Range[0,9,1](*Inthefirst32digitsofPi,eachdigitappearsatleastonce*)
Out[]=
{1,2,4,6,3,4,3,2,3,4}
Count[#][piDigits[[;;33]]]&/@Range[0,9,1](*Inthefirst33digitsofPi,eachdigitappearsatleastonce.But33>32*)
Out[]=
{1,2,5,6,3,4,3,2,3,4}
Here
a
2
=50
is the least number that each of
{0,1,2,3...,9}
appears at least twice among the first 50 digits of π:
In[]:=
Count[#][piDigits[[;;49]]]&/@Range[0,9,1](*zeroappearsonlyoncefirst49digits*)
Out[]=
{1,5,5,8,4,5,4,4,5,8}
In[]:=
Count[#][piDigits[[;;50]]]&/@Range[0,9,1](*Inthefirst50digitsofPi,eachdigitappearsatleasttwice*)
Out[]=
{2,5,5,8,4,5,4,4,5,8}
In[]:=
Count[#][piDigits[[;;51]]]&/@Range[0,9,1](*Inthefirst51digitsofPi,eachdigitappearsatleasttwicebut51>50*)
Out[]=
{2,5,5,8,4,6,4,4,5,8}
This pattern continues consistently for subsequent ordered terms in the sequence.

Import 1 billion digit of Pi

Download the Pi.zip from https://archive.org/download/Math_Constants and unzip the file. Put this notebook in the same directory with the “Pi - Dec.txt”.
SetDirectory[NotebookDirectory[]];
In[]:=
piDigit1B=ReadString["Pi - Dec.txt"];
Get rid of the leading 3 and the decimal point:
In[]:=
piDigit1B=StringTake[piDigit1B,{3,-1}];
Check the number of digits after decimal point:
In[]:=
StringLength[piDigit1B]
Out[]=
1000000000
Compare the imported value with the build-in value for the first 20 digits:
In[]:=
{StringPart[piDigit1B,1;;20],Rest[RealDigits[π,10,21][[1]]]}
Out[]=
{{1,4,1,5,9,2,6,5,3,5,8,9,7,9,3,2,3,8,4,6},{1,4,1,5,9,2,6,5,3,5,8,9,7,9,3,2,3,8,4,6}}
Use a compiled function to convert string value in to a list vector with decimal digits:
In[]:=
toDigit=FunctionCompile[Function[Typed[arg,"String"],ToCharacterCode[arg]-48]]
In[]:=
piDigits=toDigit[piDigit1B];

Compiled version of WL to generate the sequence

This is the compiled version of the existing algorithm on the OEIS page:
In[]:=
piDigitCt=FunctionCompile[Function[Typed[pd,"PackedArray"::["MachineInteger",1]],​​Module[{resq,len,count,tempSum,n=1},​​count=TypeHint[CreateDataStructure["FixedArray",10],"FixedArray"::["MachineInteger"]];​​resq=TypeHint[CreateDataStructure["Queue"],"Queue"::["MachineInteger"]];​​len=Length[pd];​​Do[Do[​​If[pd[[i]]==j-1,count["SetPart",j,1+count["Part",j]]]​​,{j,10}];​​tempSum=0;​​Do[tempSum=tempSum+If[count["Part",i]>=n,1,0],{i,10}];​​If[tempSum==10,resq["Push",i];n++]​​,{i,len}];​​resq["Elements"]]]]
Generate OEIS A383732 using 10 million digits of π and it takes less than 0.4 seconds:
In[]:=
AbsoluteTiming[piDigitCt[piDigits[[;;10^7]]];]
Out[]=
{0.381366,Null}
Generate OEIS A383732 using 1 billion digits of π and it takes around 35 seconds on modern computer:
In[]:=
AbsoluteTiming[res=piDigitCt[piDigits];]
Out[]=
{36.7261,Null}
Exam The vector length of result:
In[]:=
Length[res]
Out[]=
99986911
Compare our result with the sequence on OEIS A383732:
In[]:=
res[[;;59]]
Out[]=
{32,50,54,65,71,77,96,99,120,139,156,166,209,224,232,235,242,288,299,301,306,320,343,351,405,407,412,429,439,452,458,463,468,475,478,486,506,538,540,544,548,556,559,560,567,569,575,577,584,591,609,621,622,625,626,631,633,634,641}
In[]:=
ResourceFunction["OEISSequenceData"]["A383732"]
Out[]=
{32,50,54,65,71,77,96,99,120,139,156,166,209,224,232,235,242,288,299,301,306,320,343,351,405,407,412,429,439,452,458,463,468,475,478,486,506,538,540,544,548,556,559,560,567,569,575,577,584,591,609,621,622,625,626,631,633,634,641}
Display the last 10 entries in our result:
In[]:=
res[[-10;;]]
Out[]=
{999999917,999999922,999999929,999999932,999999944,999999953,999999961,999999965,999999967,999999986}
Use Histogram function to visualize the distribution of the gaps between adjacent terms:
In[]:=
gap=Differences[res];
In[]:=
Histogram[Differences[data],PlotLabel->"Distribution of gaps in A383732"]
Out[]=
The average gap is about 10 and median is only 7:
In[]:=
{Mean[gap],Median[gap]}//N
Out[]=
{10.0013,7.}

CITE THIS NOTEBOOK

High-performance computation of OEIS A383732 using FunctionCompile and 1 billion digits of pi​
by Shenghui Yang​
Wolfram Community, STAFF PICKS, May 21, 2025
​https://community.wolfram.com/groups/-/m/t/3464929