The Last Prime Digit

April 16, 2019

We don’t want to store lots of prime numbers, only to throw them away, so we generate the primes and calculate the last digits on the fly, trading time for space:

(define (grimes n)
  (let ((ps (primegen)) (counts (make-matrix 4 4 0)))
    (ps) (ps) (ps) ; discard 2,3,5
    (let loop ((prev (ps)) (n (- n 4)))
      (if (zero? n) counts
        (let ((curr (modulo (ps) 10)))
          (matrix-set! counts (quotient prev 3) (quotient curr 3)
            (+ 1 (matrix-ref counts (quotient prev 3) (quotient curr 3))))
          (loop curr (- n 1)))))))

This calculation took a few minutes:

> (grimes 10000000)
#(#(446808 756071 769923 526953)
  #(593195 422302 714795 769915)
  #(639384 681759 422289 756851)
  #(820368 640076 593275 446032))

You can run the program at https://ideone.com/IcT9ib.

4 Responses to “The Last Prime Digit”

Zack said
April 16, 2019 at 10:36 AM
Here is my take in Julia (1.0.2):

lastdigit(n::Int64) = n % 10

for i = 2:N
ld1 = lastdigit(P[i-1])
ld2 = lastdigit(P[i])
a = b = 0
```
if ld1 == 1
    a = 1
elseif ld1 == 3
    a = 2
elseif ld1 == 7
    a = 3
elseif ld1 == 9
    a = 4
end

if ld2 == 1
    b = 1
elseif ld2 == 3
    b = 2
elseif ld2 == 7
    b = 3
elseif ld2 == 9
    b = 4
end

if a*b != 0
    Z[a,b] += 1
end
```
end

Z /= N

4×4 Array{Float64,2}:
0.0462304 0.0742944 0.0750461 0.0544234
0.0601098 0.0444256 0.0704369 0.075029
0.0637398 0.0675519 0.0443935 0.0743187
0.0799143 0.0637294 0.0601274 0.0462292

On another note, I wonder if Dr. Grimes changed his surname so that it rhymes with “primes” because he’s so intrigued by them…

Paul said

April 18, 2019 at 6:04 AM

Sourcecode in Python is on ideone.
Here are some results:

primes < 1_000_000
 0.041   0.080   0.083   0.045
 0.056   0.036   0.074   0.085
 0.065   0.069   0.037   0.079
 0.089   0.065   0.056   0.040

10_000_000
 0.043   0.078   0.080   0.049
 0.058   0.039   0.073   0.080
 0.064   0.069   0.039   0.078
 0.085   0.065   0.058   0.042

100_000_000
 0.044   0.076   0.078   0.052
 0.059   0.042   0.072   0.078
 0.064   0.068   0.042   0.076
 0.083   0.064   0.059   0.044

1_000_000_000
 0.046   0.075   0.076   0.054
 0.060   0.044   0.071   0.076
 0.064   0.068   0.044   0.075
 0.080   0.064   0.060   0.046

10_000_000_000
 0.047   0.074   0.074   0.055
 0.060   0.046   0.070   0.074
 0.064   0.067   0.046   0.074
 0.079   0.064   0.060   0.047

 350_000_000_000 15045.708864189168
 0.049   0.072   0.072   0.057
 0.061   0.048   0.069   0.072
 0.063   0.067   0.048   0.072
 0.077   0.063   0.061   0.049

1_140_000_000_000 32903.831111384796
 0.049   0.072   0.072   0.057
 0.061   0.048   0.069   0.072
 0.063   0.066   0.048   0.072
 0.076   0.063   0.061   0.049

 3_130_000_000_000 88038.33787216355
 0.050   0.072   0.071   0.057
 0.061   0.049   0.069   0.071
 0.063   0.066   0.049   0.072
 0.076   0.063   0.061   0.050

 4_240_000_000_000 118131.35975553434
 0.050   0.072   0.071   0.057
 0.061   0.049   0.068   0.071
 0.063   0.066   0.049   0.072
 0.076   0.063   0.061   0.050

Paul said

April 20, 2019 at 7:33 AM

I let it run longer and with primes upto 10_900_000_000_000 the result only changes a little.

 10_900_000_000_000 298241 seconds
 0.050   0.071   0.071   0.058
 0.061   0.049   0.068   0.071
 0.063   0.066   0.049   0.071
 0.075   0.063   0.061   0.050

Daniel said

April 28, 2019 at 11:46 PM

Here’s a solution in C++. It uses the non-segmented Sieve of Eratosthenes, limiting how large the primes can get before the system runs out of memory.

#include <cstdlib>
#include <cstring>
#include <iomanip>
#include <iostream>
#include <string>
#include <vector>

using namespace std;

typedef unsigned long long ull;

int main(int argc, char* argv[]) {
  if (argc != 2) return EXIT_FAILURE;
  ull n = stoull(argv[1], NULL, 10);
  vector<bool> sieve(n + 1, true);
  sieve[0] = false;
  sieve[1] = false;
  for (ull i = 2; i * i <= n; ++i) {
    for (ull j = i * 2; j <= n; j += i) {
      sieve[j] = false;
    }
  }
  int idxs[10] = {-1, 0, -1, 1, -1, -1, -1, 2, -1, 3};
  ull total = 0;
  ull counts[4][4];
  memset(counts, 0, sizeof(counts));
  // ignore 2 and 5 by initializing 'last' to 3 and
  // starting iteration at 7.
  int last = 3;
  for (ull i = 7; i <= n; ++i) {
    if (!sieve[i]) continue;
    int current = i % 10;
    ++counts[idxs[last]][idxs[current]];
    ++total;
    last = current;
  }
  cout << fixed;
  cout << setprecision(3);
  for (int i = 0; i < 4; ++i) {
    for (int j = 0; j < 4; ++j) {
      if (j > 0) cout << " ";
      cout << (float)counts[i][j] / total;
    }
    cout << endl;
  }
  return EXIT_SUCCESS;
}

Example Usage:

$ ./praxis 1000000
0.041 0.080 0.083 0.045
0.056 0.036 0.074 0.085
0.065 0.069 0.037 0.079
0.089 0.065 0.056 0.040

Programming Praxis