## 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.

Pages: 1 2

### 4 Responses to “The Last Prime Digit”

1. Zack said

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…

2. Paul said

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
```
3. Paul said

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
```
4. Daniel said

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
```