As many near-repdigit primes as possible

Question

A near-repdigit number is a positive integer where all the digits are the same, except one. For example 101 and 227 are near-repdigits. A near-repdigit prime is a near-repdigit that is also prime. For example:

101, 113, 131, 151, 181, 191, 199, 211, 223, 227, 229, 233, 277, 311, 313, 331, 337, 353, 373, 383, 433, 443, 449, 499, 557, 577, 599, 661, 677, 727, 733, 757, 773, 787, 797, 811, 877, 881, 883, 887, 911, 919, 929, 977, 991, 997, 1117, 1151, 1171, 1181, 1511

These are all near-repdigit primes.

The smallest near-repdigit prime has two digits, but it is an open mathematical question whether there exists a near-repdigit prime for every possible number of digits larger than 2. It is conjectured, by me, that the answer is in fact yes.

Task

For each integer n where n >= 2, compute a near-repdigit prime with n digits.

Your code can either output a near-repdigit prime with 2, 3, ... digits or, as they are easily compressible, output a compressed representation of each number. For example, any near-repdigit can be represented by four smaller numbers. The first is the number of digits, the second the majority digit, the third the minority digit and the fourth the location of the minority digit. You can choose whichever representation you prefer.

Primality testing

There are many different ways to test if a number is prime. You can choose any method subject to the following conditions.

• You can use any primality test that is guaranteed never to make a mistake.
• You can use any well-known primality test which hasn't been proved to be correct but for which there is no known example number for which it gives the wrong answer.
• You can use a probabilistic primality test if the probability of giving the wrong answer is less than 1/1000.

I will test your code for 5 minutes on my ubuntu 22.04 machine, but please quote how high you get on your machine in your answer.

This challenge is judged per language.

In the very unlikely event that you find a number of digits n for which there is no near-repdigit prime, I will award a 500-rep bounty and you will get mathematical immortality.

Results so far

• n=1291 by Kirill L. in Julia
• n=1291 by c-- in C with gmp
• n=1232 by jdt in C++ with gmp
• n=972 by Kirill L. in Julia
• n=851 by alephalpha in Pari/GP
• n=770 by ZaMoC in Wolfram Language (not tested on my PC)
• n=722 by gsitcia in Pypy/Python
• n=721 by jdt in C++
• n=665 by c-- in C++
• n=575 by Seggan in Kotlin
• n=403 by Arnauld in nodejs
• n=9 by py3_and_c_programmer in Python

Answer 1

Julia n = 972 (single-threaded)

using Primes, ThreadPools

function nrdprime(n) 
  a = fill(0, n)
  for j = 1:9
    a[1:n] .= j
    for k = n:-1:1, l = 0:9
      if j == l continue end
      if l == 0 && k == 1 continue end
      a[k] = l
      if iseven(a[n]) continue end
      x = parse(BigInt, join(a))
      if isprime(x)
        print("$n $j $l $k\n")
        return x
      end
      a[k] = j
    end
  end
  println(n, " has no answer")
end

@qthreads for i = 2:2000 nrdprime(i) end

Uses Primes.jl package (import Pkg; Pkg.add("Primes")). According to the source code, it relies on Miller-Rabin algorithm implementation from libgmp library.

UPDATE: Added multithreading with ThreadPools.jl (import Pkg; Pkg.add("ThreadPools")). Compared to Base.Threads module, using a pool allows scheduling calculations in a more convienient way for our use case. In order to take advantage of concurrency, run the program with julia --threads N, setting the number N in accordance with your CPU specification.

The output was changed to compressed form (n major minor position) and I also tweaked the code a bit to reduce the number of array allocations, but it looks like this isn't really a bottleneck.

Try it online! for n up to about 550 in 1 min, using an ancient version of Julia that still has primes-related functionality integrated in base, and thus available on TIO.

Answer 2

C++ (gcc) with boost, n = 595

~354 digits in 60 seconds on TIO.

~595 digits in 5 minutes on my Windows machine, compiled with MSVC.

#include <atomic>
#include <chrono>
#include <iostream>
#include <mutex>
#include <thread>
#include <unordered_map>
#include <boost/multiprecision/miller_rabin.hpp>

namespace mp = boost::multiprecision;

std::atomic<int> digits{ 2 };
std::mutex mtx;
std::unordered_map<int, mp::cpp_int> results;

mp::cpp_int calcPrimes(int n)
{
    char buf[10000] = { 0 };
    for (int i = '1'; i <= '9'; i+=2)
    {
        std::memset(buf, i, n);
        for (int j = '0'; j <= '9'; j++)
        {
            if (j != i)
            {
                for (int k = 0; k < n; k++)
                {
                    if (k > 0 || j > '0')
                    {
                        buf[k] = j;
                        mp::cpp_int num(buf);
                        if (miller_rabin_test(num, 10))
                        {
                            return num;
                        }
                        buf[k] = i;
                    }
                }
            }
        }
    }
    return 0;
}

void threadJob()
{
    for (;;)
    {
        int n = digits.fetch_add(1);
        auto result = calcPrimes(n);
        mtx.lock();
        results[n] = result;
        mtx.unlock();
    }
}

int main() {
    int count = 1;
    const auto start = std::chrono::steady_clock::now();
    const auto tc = std::thread::hardware_concurrency();
    std::vector<std::thread> threads;
    for (int j = 0; j < tc; j++)
        threads.emplace_back(threadJob);
    for (;;)
    {
        std::this_thread::sleep_for(std::chrono::milliseconds(50));
        mtx.lock();
        auto seconds = std::chrono::duration_cast<std::chrono::seconds>(std::chrono::steady_clock::now() - start);
        for (const auto& result : results) {
            std::cout << std::setfill(' ') << std::setw(3) << result.first << ": " << result.second;
            std::cout << "  time: " << seconds.count() << " sec, count: " << count++ << "\n";
            if (result.second == 0)
            {
                std::cout << "Mathematical Immortality!!!";
                exit(1);
            }
        }
        results.clear();
        mtx.unlock();
    }
    for (int j = 0; j < tc; j++)
        threads[j].join();
}

Try it online!

Answer 3

Vyxal, at least n=116

ƛkdn1-v*9ʀẊ'ƒc¬;9ʀẊƛf÷$Ṁ;U'⌊=$⌊æ∧;h

Try it Online!

A simple modification of my answer to a CMC for near-rep numbers of length 8.

It can do up to n=116 in a minute on the online interpreter. Maybe more depending on where you run it.

The score here is cumulative total. If you want individual limits of n, it can handle at least n=200 in less than 10 seconds with a few modifications.

I'll explain this properly a bit later, but basically:

kd7v*9ʀẊ gets the Cartesian product of each digit from 0 to 9 repeated n times and the range [0, 9].

'ƒc¬; keeps items where the digit from the range isn't in the string (because this Cartesian product is [repeated digit, number to insert at every position])

9ʀẊ is once again Cartesian product with range [0, 9]. It gives a list of [[repeated digit, number to insert], location to insert]

ƛf÷$Ṁ; flattens each list, swaps the order of the number to insert and the position and performs the insertion

Finally, '⌊=$⌊æ∧; uniquifies the result and keeps only those invariant under conversion to int and which are also prime.

Answer 4

JavaScript (Node.js), n = 351

Reached \$n=351\$ in 5 minutes on my laptop. Reaches \$n\approx 250\$ in 1 minute on TIO.

Uses a Miller-Rabin primality test with 10 rounds.

NB: Because the challenge allows a probability of false positive as high as 1/1000, we could actually perform fewer rounds -- at least until \$n\$ becomes really large. This would however not improve the performance drastically, as most of the computation is used for processing non-primes for which the test is most likely to fail on the first round anyway.

const ROUNDS = 10;

function modpow(k, e, m) {
  let r = 1n;

  while(e) {
    if(e & 1n) {
      r = (r * k) % m;
    }
    e >>= 1n;
    k = (k * k) % m;
  }
  return r;
}

function millerRabin(sz, n, k) {
  let d = n - 1n, s = 0n;
  do {
    s++;
    d >>= 1n;
  } while(!(d & 1n));

  while(k--) {
    let a;

    do {
      a = BigInt([...Array(sz)].map(_ => Math.random() * 10 | 0).join(''));
    } while(a < 2n || a > n - 1n);

    if(millerRabinIteration(n, d, s, a)) {
      return false;
    }
  }
  return true;
}

function millerRabinIteration(n, d, s, a) {
  let x = modpow(a, d, n);

  if(x == 1n || x == n - 1n) {
    return false;
  }

  for(let i = 0n; i < s - 1n; i++) {
    x = modpow(x, 2n, n);
    if(x == n - 1n) {
      return false;
    }
  }
  return true;
}

let ts = Date.now(), time;

for(let sz = 2; (time = Math.round((Date.now() - ts) / 1000)) < 300; sz++) {
  console.log(time + 's', 'n=' + sz, search(sz));
}

function search(sz) {
  for(let i = 1; i <= 9; i++) {
    let a = Array(sz).fill(i);

    for(let j = 0; j <= 9; j++) {
      for(let k = j ? 0 : 1; i !=j && k < sz; k++) {
        a[k] = j;
        let n = BigInt(a.join(''));
        if(
          n % 2n && n % 3n && n % 5n && n % 7n &&
          millerRabin(sz, n, ROUNDS)
        ) {
          return n;
        }
        a[k] = i;
      }
    }
  }
  return false;
}

Try it online!

Answer 5

Python 3 (PyPy), n = 808

~432 digits on TIO in one minute (no multiprocessing)

~808 digits on my computer in 5 minutes

~697 digits on my computer in 5 minutes on 1 core

Uses Baillie-PSW primality testing.

Avoids testing any multiple of 2, 3 or 5.

# aab -- a in [1,2,3,4,5,6,7,8,9]
#     -- b in [1,3,7,9]
# aba -- a in [1,3,7,9]
#     -- b in [0,1,2,3,4,5,6,7,8,9]
# baa -- a in [1,3,7,9]
#     -- b in [1,2,3,4,5,6,7,8,9]
# b + a*(n-1) != 0 mod 3

aab = [[(a,b) for a in [1,2,3,4,5,6,7,8,9] for b in [1,3,7,9] if (b + a*(n-1))%3] for n in range(3)]
aba = [[(a,b) for a in [1,3,7,9] for b in [0,1,2,3,4,5,6,7,8,9] if (b + a*(n-1))%3] for n in range(3)]
baa = [[(a,b) for a in [1,3,7,9] for b in [1,2,3,4,5,6,7,8,9] if (b + a*(n-1))%3] for n in range(3)]

def f(n):
    A = (10**n-1)//9
    # aab
    for a,b in aab[n%3]:
        if baillie_psw(a*A + (b-a)): return (n,a,b,0)
    # baa
    B = 10**(n-1)
    for a,b in baa[n%3]:
        if baillie_psw(a*A + (b-a)*B): return (n,a,b,n-1)
    # aba
    for k in range(1, n-1):
        B = 10**k
        for a,b in aba[n%3]:
            if baillie_psw(a*A + (b-a)*B): return (n,a,b,k)
    sike

def main():
    import multiprocessing
    import time
    t = time.time()
    for i in range(3, 2500, 2):
        if baillie_psw(i): small_primes.append(i)
    with multiprocessing.Pool(multiprocessing.cpu_count()-1 or 1) as pool:
        for i in pool.imap(f, range(2, 10000)):
            print(i, time.time() - t)

def baillie_psw(n):
    for p in small_primes:
        if n == p: return True
        if n%p == 0: return False
    return strong_pseudoprime(n) and lucas_prime(n)

small_primes = [2] # calculated out to 3000 in main

def ds(d):
    s = 0
    while d%2 == 0:
        s += 1
        d //= 2
    return d,s

def strong_pseudoprime(n):
    n1 = n-1
    d,s = ds(n1)
    x = pow(2,d,n)
    if x == 1 or x == n1: return True
    for _ in range(s):
        x = (x*x)%n
        if x == n1: return True
    return False

# code from https://en.wikipedia.org/wiki/Jacobi_symbol
def jacobi(D, n):
    D %= n
    t = 1
    while D:
        while D%2 == 0:
            D //= 2
            r = n%8
            if r == 3 or r == 5: t = -t
        n,D = D,n
        if D%4 == 3 == n%4: t = -t
        D %= n
    if n == 1: return t
    return 0

def sign(n):
    if n > 0: return 1
    if n < 0: return -1
    return 0

def isqrt(n):
    if n < 2: return n
    k = n//2
    k1 = n
    while k < k1:
        k1 = k
        k = (k+n//k)//2
    return k

def is_square(n):
    np = isqrt(n)
    return np*np == n

def UVQ(P, Q, D, k, m):
    u,v,q = 0,0,-1
    B = f'{k:b}'
    B1 = B.rstrip('0')
    for i in B1:
        u,v,q = (u*v)%m,(v*v-2*q)%m,(q*q)%m
        if i == '1':
            a = P*u + v
            b = D*u + P*v
            if a%2: a += m
            if b%2: b += m
            u,v,q = (a//2)%m, (b//2)%m, (Q*q)%m
    if u == 0 or v == 0: return True
    for _ in range(len(B)-len(B1)-1):
        u,v,q = (u*v)%m,(v*v-2*q)%m,(q*q)%m
        if v == 0: return True
    return False

def lucas_prime(n):
    for D in [5, -7, 9, -11, 13, -15, 17, -19, 21, -23, 25, -27, 29, -31, 33, -35, 37, -39, 41, -43]:
        J = jacobi(D, n)
        if J == -1: break
        if J == 0 and D%n: return False
    else:
        if is_square(n): return False
        else:
            while J != -1:
                D = -D - 2*sign(D)
                J = jacobi(D, n)
                if J == 0 and D%n: return False
    P = 1
    Q = (1-D)//4
    if Q == -1: P = Q = 5
    return UVQ(P, Q, D, n+1, n)

if __name__ == '__main__':
    main()

Try it online!

In profiling I found that about 97% of the time was spent in strong_pseudoprime, specifically the line x = pow(2,d,n).

I wanted to avoid using external libraries in my answer, but it is possible to improve the speed significantly by using gmpy2 (or gmpy_cffi with pypy), and replacing the line with x = pow(mpz(2),d,n). This achieved ~1050 digits in 5 minutes.

Answer 6

Wolfram Language (Mathematica), n=770

(n = #;
  r = RandomChoice[{1,3,7,9}];
  While[! 
     PrimeQ[p = 
       FromDigits@
        RandomSample[
         Join[Table[r, 
           n - 2], {r1 = 
            RandomChoice@Complement[Range[0, 9], {r}]}, {r}]]] || 
    r1 == r];
  i = IntegerDigits@p;
  f = First /@ Reverse@SortBy[Tally@i, Last];
  Flatten[{n, f, Position[i, f[[2]]]}]) &

Try it online!

n=500 in 1min
n=770 in 5min

Answer 7

Kotlin, n = 575

fun main() {
    var n = 2
    while (true) {
        println(primesForN(n))
        n++
    }
}

fun primesForN(n: Int): String {
    for (major in '1'..'9') {
        // Java's String.repeat is much faster than Kotlin's
        @Suppress("PLATFORM_CLASS_MAPPED_TO_KOTLIN")
        val digits = (major.toString() as java.lang.String).repeat(n - 1)
        for (minor in '0'..'9') {
            if (major == minor) continue
            if (major.isEven()) {
                if (minor.isEven()) continue
                val number = digits + minor
                if (number.toBigInteger().isProbablePrime(10)) {
                    return number
                }
            } else {
                for (i in 0 until n - 1) {
                    if ((i == n - 1 && minor.isEven()) || (i == 0 && minor == '0')) continue
                    val number = digits.substring(0, i) + minor + digits.substring(i)
                    if (number.toBigInteger().isProbablePrime(10)) {
                        return number
                    }
                }
            }
        }
    }
    return "No near-repdigit primes found for n = $n"
}

fun Char.isEven() = this == '0' || this == '2' || this == '4' || this == '6' || this == '8'

Try it online!

Gotta love Kotlin's CharRange. The isProbablePrime test is a Miller-Rabin test that has a \$\frac{1}{2^{10}}\$ chance of failing. For ease, here is the jar file.

Answer 8

C++ (clang) with GMP n = 1164

To be compiled with the following flags: -O3 -lgmp -lpthread.

#include <atomic>
#include <chrono>
#include <cstring>
#include <iostream>
#include <mutex>
#include <thread>
#include <map>
#include <vector>
#include <gmp.h>

struct result
{
    char major;
    char minor;
    int pos;
};

std::atomic<int> digits{ 2 };
std::atomic<bool> timeOut{ false };
std::mutex mtx;
std::map<int, result> results;

result calcPrimes(int n)
{
    mpz_t x;
    char buf[10000] = { 0 };
    for (char i = '1'; i <= '9'; i+=2)
    {
        std::memset(buf, i, n);
        for (char j = '0'; j <= '9'; j++)
        {
            if (j != i)
            {
                for (int k = 0; k < n; k++)
                {
                    if (k > 0 || j > '0')
                    {
                        buf[k] = j;
                        mpz_init_set_str(x, buf, 10);
                        if (mpz_probab_prime_p(x, 10))
                        {
                            return { i, j, k };
                        }
                        buf[k] = i;
                    }
                }
            }
        }
    }
    return {};
}

void threadJob()
{
    while (!timeOut)
    {
        int n = digits.fetch_add(1);
        auto result = calcPrimes(n);
        mtx.lock();
        results[n] = result;
        mtx.unlock();
    }
}

int main() {
    int count = 0;
    int maxTime = 5 * 60;
    const auto start = std::chrono::steady_clock::now();
    const auto tc = std::thread::hardware_concurrency();
    std::vector<std::thread> threads;
    for (int j = 0; j < tc; j++)
        threads.emplace_back(threadJob);
    for (;;)
    {
        std::this_thread::sleep_for(std::chrono::milliseconds(50));
        mtx.lock();
        auto seconds = std::chrono::duration_cast<std::chrono::seconds>(std::chrono::steady_clock::now() - start);
        if (seconds.count() >= maxTime)
        {
            timeOut = true;
            mtx.unlock();
            break;
        }
        for (const auto& result : results) {
            std::cout << result.first << ": " << "major: " << result.second.major << ", minor: " << result.second.minor << ", pos: " << result.second.pos;
            std::cout << ", time: " << seconds.count() << " sec, count: " << ++count << "\n";
        }
        results.clear();
        mtx.unlock();
    }
    for (int j = 0; j < tc; j++)
        threads[j].join();
    std::cout << "Time: " << maxTime << " seconds\n";
    std::cout << "Count: " << count << "\n";
}

Try it online!

Answer 9

Python 3, computed till n=8

"Takes 51 seconds for 8 digits"

Not a very interesting answer, I fear, and it's very slow.

from itertools import product as pt

def isprime(n):
    # This was borrowed from https://en.wikipedia.org/wiki/Primality_test
    # Modified to save bytes unnecessarily
    if n <= 3:return n > 1
    if not n%2 or not n%3:return False
    for i in range(5, int(n**.5)+1, 6):
        if not n%i or not n%(i+2):return False
    return True

def isnrep(n):
    if (not isprime(n)) or not (len(set(str(n)))==2):return False
    return list(str(n)).count(list(set(str(n)))[0])in[len(str(n))-1, 1]
    
def getnrepfor(n):
    x=[list(i)for i in eval("list(pt((1,2,3,4,5,6,7,8,9),"+("(0,1,2,3,4,5,6,7,8,9), "*(n-2))+"(1,3,7,9)))")]
    for i in range(len(x)):
        for j in range(len(x[i])):
            x[i][j]=str(x[i][j])
    x=[int("".join(i)) if isnrep(int("".join(i))) else '' for i in x]
    for i in x:
        if i!='':return i
    return 'It has been proven that for %d, there does not exist a near-repdigit' % n

Back to the bounty... I've found out that the amount of near-repdigit primes actually increases when the number of digits is increased. So I guess nobody will get that bounty.

Uses the 6k+-1 computation of primes.

I've updated this with an approach similar to lyxal's answer, but it hasn't produced much of a speedup. Oh well, guess esoteric languages are faster than those around 1991.

Answer 10

C (gcc) + GMP, n = 1015 (w/ -j2)

• Compile with same flags as in TIO.
• Specify number of threads as -j jobs depending on your hardware.

This uses mpz_probab_prime_p() from GMP which uses a combination of a Baille-PSW test, and then Miller-Rabin tests.

Output: digits majority minority location(1-indexed).

#define _XOPEN_SOURCE 700
#include <aio.h>
#include <gmp.h>
#include <pthread.h>
#include <stdatomic.h>
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>

typedef struct RepDigit {
    unsigned long maj;
    unsigned long min;
    unsigned long pos;
} RepDigit;

static RepDigit calc_primes(unsigned long digits);
static void init_ones(mpz_t rop, unsigned long digits);
static void *job(void *);

RepDigit
calc_primes(unsigned long digits)
{
    mpz_t rep, digit, num;
    unsigned long i, j, k;

    mpz_inits(rep, digit, NULL);
    init_ones(rep, digits);
    mpz_init_set(num, rep);
    mpz_add(rep, rep, rep);

    for (i = 1; i <= 9; mpz_add(num, num, rep), i += 2ul * (1ul + (i == 3))) {
        mpz_set_ui(digit, 10);
        for (k = 1; k < digits; mpz_mul_ui(digit, digit,  10), ++k) {
            mpz_submul_ui(num, digit, i);
            for (j = k == digits - 1; j <= 9; ++j) {
                if (j == i)
                    continue;
                mpz_addmul_ui(num, digit, j);
                if (mpz_probab_prime_p(num, 5))
                    goto end;
                mpz_submul_ui(num, digit, j);
            }
            mpz_addmul_ui(num, digit, i);
        }
    }

end:
    mpz_clears(rep, digit, num, NULL);
    return (RepDigit) {i, j, digits - k};
}

void
init_ones(mpz_t rop, unsigned long digits)
{
    mpz_ui_pow_ui(rop, 10, digits);
    mpz_sub_ui(rop, rop, 1);
    mpz_divexact_ui(rop, rop, 9);
}

void *
job(void *arg)
{
    RepDigit p;
    unsigned long n;
    char buf[2048][16];
    struct aiocb aiocbs[2048], *out;

    static atomic_ulong digits = 2;

    (void) arg;
    for (out = &aiocbs[0];; ++out) {
        n = atomic_fetch_add_explicit(&digits, 1, memory_order_relaxed);
        p = calc_primes(n);
        /* digits majority minority location(1-indexed) */
        out->aio_fildes = STDOUT_FILENO;
        out->aio_offset = 0;
        out->aio_buf = buf[out - aiocbs];
        out->aio_nbytes = snprintf(buf[out - aiocbs], sizeof buf[0], "%lu %lu %lu %lu\n", n, p.maj, p.min, p.pos);
        out->aio_reqprio = 0;
        out->aio_sigevent.sigev_notify = SIGEV_NONE;
        aio_write(out);
    }

    /* unreachable */
    return NULL;
}

int
main(int argc, char *argv[])
{
    int opt, jobs;
    pthread_t t;
    extern char *optarg;

    jobs = 1;
    while ((opt = getopt(argc, argv, "j:")) > 0) {
        switch (opt) {
        case 'j':
            jobs = atoi(optarg);
            if (jobs > 0)
                break;
            /* fallthrough */
        default:
            fprintf(stderr, "usage: %s [-j jobs]\n", argv[0]);
            exit(EXIT_FAILURE);
        }
    }

    while (jobs --> 1) {
        pthread_create(&t, NULL, &job, NULL);
        pthread_detach(t);
    }

    job(NULL);

    /* unreachable */
    return 0;
}

Try it online! n ~ 700 (w/ -j2).

Answer 11

Python, with gmpy2, n=920

~920 digits in 5 minutes on my Windows machine

import gmpy2
import time

def nrdprime(n):
    for j in range(1,10):
        for k in range(n, 0, -1):
            for l in range(10):
                if j == l: continue
                if l == 0 and k == 1: continue
                a = [j] * n
                a[k-1] = l
                if a[n-1] % 2 == 0: continue
                x = int(''.join(map(str, a)))
                if gmpy2.is_prime(x):
                    print(n, j, l, k)
                    return x
    print(n, "has no answer")

start_time = time.time()
i = 2

while time.time() - start_time < 300:  # Run for 300 seconds (5 minutes)
    nrdprime(i)
    i += 1

print(f"Processed up to n = {i-1} in 5 minutes.")

Answer 12

PARI/GP, 164 bytes, n = 839

f(n)=d=(10^n-1)/9;for(i=1,9,for(j=0,n-1,for(k=0,9,if(i!=k&&(k>0||j<n-1)&&((i*n+k-i)%3)&&ispseudoprime(p=d*i+(k-i)*10^j),return(p)))))
for(n=2,oo,print(n" -> "f(n)))

Attempt This Online!

The ispseudoprime function uses the BPSW primality test. Currently no counterexample is known.

The package name for PARI/GP is pari-gp on Ubuntu, Debian, Arch Linux and some other distros.

Runs as gp -qf ./file_name.gp.

Answer 13

C++ (gcc), n = 595

• Compile with -O3 -pthread -lfmt

#include <atomic>
#include <boost/multiprecision/cpp_int.hpp>
#include <boost/multiprecision/miller_rabin.hpp>
#include <fmt/core.h>
#include <thread>

namespace mp = boost::multiprecision;

struct RepDigit {
    int pos;
    char maj, min;
};

static std::atomic<int> digits{2};

static RepDigit
calcPrimes(int digits)
{
    static const mp::cpp_int ten{10};
    char buf[2048] = {};
    for (char i = '1'; i <= '9'; i += 2) {
        std::memset(buf, i, digits);
        mp::cpp_int num{buf};
        mp::cpp_int n{ten};
        for (int k = 1; k < digits; n *= ten, ++k) {
            for (char j = '9'; j >= '0'; --j) {
                if (j == i)
                    continue;
                num += n * (j - i);
                if (mp::miller_rabin_test(num, 10))
                    return {digits - k, i, j};
                num -= n * (j - i);
            }
        }
    }
    return {};
}

static void
job()
{
    for (;;) {
        auto n = digits.fetch_add(1);
        auto p = calcPrimes(n);
        /* digits majority minority location(1-indexed) */
        fmt::print("{} {} {} {}\n", n, p.maj, p.min, p.pos);
    }
}

int
main()
{
    for (auto hc = std::thread::hardware_concurrency(); hc > 1; --hc)
        std::thread{job}.detach();

    job();
}

No TIO link because it doesn't have {fmt}.