$MaxExtraPrecision  =  2000000

"You must change the precision in order to compute all of the \
necessary large primes. You must also load the 
NumberTheory package in order to be able to use the count primes \
functions.";

<< NumberTheory`NumberTheoryFunctions`

"Starting value is either 1 or 10^8 for the first prime The \
starting value is what prime you want to start your data at.";
Start =  PrimePi[10^8] + 1 + 1000

"Dataset is the array of all of the data. The different primes are \
the rows and the digit distributions are the columns.
The columns are the frequencies of the digits 1-100 separately, and then the \
frequencies of all digits on the intervals 1-100, 1-1000, 1-50000, 1-100000, \
1-1000000. Dataset merely computes the frequencies of the digits in the first \
50,000 digits of the continued fraction."

Array[Dataset,  {1000,  120}];

\(\(\(
\)\(Start = \ \((5761456 + 1000)\); \ 
    iterations\  = \ 1000; \ Num\  = \ 50000;

    \(For[k = Start, \ 
        k ≤ \ 
          Start + iterations, \ \(k++\), 
For[i\  = \ 1, \ 
          i\  ≤ \ 120, \ \(i++\), \ 
          NumDigit[i\ ]\  = \ 0\ \ ]; \ 
\ p = Prime[k]; \ 
        a = ContinuedFraction[\((p^\((1/3)\))\), \ 
            Num]; 
For[i\  = \ 1, \ 
          i\  ≤ \ 
            Num, \ \(i++\), \ \ \ \ \ \ 
\ \ \ \ \ \(If[\ 
              a[\([i]\)]\  < \ 101, \ 
              NumDigit[a[\([i]\)]\ ]\  = \ 
                NumDigit[a[\([i]\)]\ ]\  + \ 1, 
	\(If[\ 
                  a[\([i]\)]\  < \ 1001, \ 
                  NumDigit[\ 101\ ]\  = \ 
                    NumDigit[\ 101\ ]\  + \ 
                      1, \ 
	\ \ \ \ \ \(If[\ \ a[\([i]\)]\  < \ 10001, \ NumDigit[102\ ]\  = \ 
                        NumDigit[102]\  + \ 
                          1, \ 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \(If[\ a[\([i]\)] < 50001, \ 
                          NumDigit[103]\  = \ 
                            NumDigit[103]\  + \ 
                              1, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \(If[a[\([i]\)] < 100001, \ 
                              NumDigit[104] = 
                                NumDigit[104] + 
                                  1, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ If[a[\([i]\)] < 1000000, \ 
                                NumDigit[105] = NumDigit[105] + 1, \ 
                                NumDigit[106] = 
                                  NumDigit[106] + 
                                    1]\ \ \ \ \ \ \ 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ];\)\ \ \ 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ];\)
\ \ \ \ \ \ \ \ \ \ \ \ \ \ ];\)
\ \ \ \ \ \ \ \ \ \ ];\)
\ \ \ \ \ ];\)
\ \ \ \ ]; 

        For[j\  = \ 1, \ 
          j\  ≤ \ 106, \ \(j++\), \ \(Dataset[k]\)[j] = 
            NumDigit[j]\ ];
];\)\)\)\)

"secondDataset is necessary to compute the \
frequencies of the digits to 100,000 digits. Dataset is subtracted from \
secondDataset in order to find the frequencies of digits from the 50,001 \
digit to the 100,000 digit.";

Array[secondDataset,  {1000,  120}];
Num2  =  100000;
\(\(For[k = Start, \ 
        k ≤ \ 
          Start + iterations, \ \(k++\), 
For[i\  = \ 1, \ 
          i\  ≤ \ 120, \ \(i++\), \ 
          NumDigit2[i\ ]\  = \ 0\ \ ]; \ 
\ p = Prime[k]; \ 
        a = ContinuedFraction[\((p^\((1/3)\))\), \ 
            Num2]; 
For[i\  = \ 1, \ 
          i\  ≤ \ 
            Num2, \ \(i++\), \ \ \ \ \ \ 
\ \ \ \ \ \(If[\ 
              a[\([i]\)]\  < \ 101, \ 
              NumDigit2[a[\([i]\)]\ ]\  = \ 
                NumDigit2[a[\([i]\)]\ ]\  + \ 1, 
	\(If[\ 
                  a[\([i]\)]\  < \ 1001, \ 
                  NumDigit2[\ 101\ ]\  = \ 
                    NumDigit2[\ 101\ ]\  + \ 
                      1, \ 
	\ \ \ \ \ \(If[\ \ a[\([i]\)]\  < \ 10001, \ NumDigit2[102\ ]\  = \ 
                        NumDigit2[102]\  + \ 
                          1, \ 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \(If[\ a[\([i]\)] < 50001, \ 
                          NumDigit2[103]\  = \ 
                            NumDigit2[103]\  + \ 
                              1, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \(If[a[\([i]\)] < 100001, \ 
                              NumDigit2[104] = 
                                NumDigit2[104] + 
                                  1, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ If[a[\([i]\)] < 1000000, \ 
                                NumDigit2[105] = NumDigit2[105] + 1, \ 
                                NumDigit2[106] = 
                                  NumDigit2[106] + 
                                    1]\ \ \ \ \ \ \ 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ];\)\ \ \ 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ];\)
\ \ \ \ \ \ \ \ \ \ \ \ \ \ ];\)
\ \ \ \ \ \ \ \ \ \ ];\)
\ \ \ \ \ ];\)
\ \ \ \ ]; 

        For[j\  = \ 1, \ 
          j\  ≤ \ 106, \ \(j++\), \ \(secondDataset[k]\)[j] = 
            NumDigit2[
              j]\ ];
];\)\)

\(\(For[k\  = Start, \ 
        k\  ≤ \ \ Start\  + \ 
            iterations, \ \(k++\), \ 
		\(For[j = 1, \ 
            j ≤ \ 
              106, \ \(j++\), \ \(completedata[k]\)[
                j] = \(secondDataset[k]\)[j] - \(Dataset[k]\)[
                  j]];\)
];\)\)

\(\("The s table provides the matrix of all of the data. The primes change by the rows and the digit distributions are in the columns.";\)
\)
s = 
        Table[completedata[i][j],  {i, 1000},  {j, 
            8}];
MatrixForm[s]