19
\$\begingroup\$

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
\$\endgroup\$
14
  • 1
    \$\begingroup\$ Relevant oeis, also related. \$\endgroup\$ Nov 10, 2022 at 21:27
  • \$\begingroup\$ I'd expect there not to be a near-repdigit prime of length n for every n - there are only 81n near-repdigits of length n, and the density of primes is 1/log(n). Asymptotically, we should expect almost none. \$\endgroup\$ Nov 10, 2022 at 21:48
  • \$\begingroup\$ @Spitemaster It will be interesting (and expensive for me) if anyone can find an example. \$\endgroup\$
    – user108721
    Nov 10, 2022 at 21:58
  • 3
    \$\begingroup\$ @Spitemaster With 81n near-repdigits of decimal length n, which have a roughly 1 / ln(10^n) density of primes, don't we expect on the order of 81 / ln (10) ≈ 35 near-repdigit primes of length n? \$\endgroup\$
    – xnor
    Nov 10, 2022 at 23:34
  • 1
    \$\begingroup\$ @Noodle9 I defined it as "near-repdigit number is a positive integer where all the digits are the same, except one." So it's two under that definition. \$\endgroup\$
    – user108721
    Nov 11, 2022 at 13:19

13 Answers 13

5
\$\begingroup\$

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.

\$\endgroup\$
6
  • \$\begingroup\$ Your answer is by far the fastest currently! It's also single threaded so this is even more impressive. \$\endgroup\$
    – user108721
    Nov 16, 2022 at 11:24
  • \$\begingroup\$ Sadly overtaken by C++. Can Julia catch up again? \$\endgroup\$
    – user108721
    Nov 16, 2022 at 17:44
  • \$\begingroup\$ @graffe, I doubt, as it uses the same GMP library under the hood, and Julia has some overhead as a higher-level language. But it's still interesting to see how close we can get, so I also added multithreading. \$\endgroup\$
    – Kirill L.
    Nov 16, 2022 at 19:44
  • \$\begingroup\$ The output from your new code is confusing. If I run it with 8 threads for 5 seconds I get e.g. bpa.st/QOQA What do the last two lines mean? \$\endgroup\$
    – user108721
    Nov 18, 2022 at 11:18
  • 1
    \$\begingroup\$ @graffe, it looks like a race condition with overlapping print statements from different threads. Which is strange, because I don't observe it neither on Win11, nor on Ubuntu. Maybe try it with an updated print statement (see code) that at least prints the entire line in one go. \$\endgroup\$
    – Kirill L.
    Nov 18, 2022 at 14:09
2
\$\begingroup\$

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!

\$\endgroup\$
10
  • 1
    \$\begingroup\$ suggest using a thread pool or any other mechanism which allows you to not have to wait for every thread to finish before starting the next batch \$\endgroup\$
    – c--
    Nov 13, 2022 at 5:27
  • \$\begingroup\$ @c-- thanks for the good suggestion! My implementation could be better but I have to mow the lawn. Feel free to suggest some improvements :-) \$\endgroup\$
    – jdt
    Nov 13, 2022 at 14:40
  • 1
    \$\begingroup\$ There's no need for results to be an unordered_map, it could be a vector or maybe a list (unlikely), trying to be mindful of unnecessary moves and copies since mp::cpp_int isn't trivially destructible. The thread pool is kind of a hack but it gets the job done. \$\endgroup\$
    – c--
    Nov 13, 2022 at 17:04
  • 1
    \$\begingroup\$ AMD Ryzen 5 3400G with g++ 11.3.0 \$\endgroup\$
    – user108721
    Nov 13, 2022 at 19:30
  • 1
    \$\begingroup\$ It was 722 which I hadn't noticed until I started using awk 'NR > 1 && $0 > old + 1 { print old; exit } { old = $0 }' to check :) \$\endgroup\$
    – user108721
    Nov 13, 2022 at 20:20
1
\$\begingroup\$

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.

\$\endgroup\$
1
  • \$\begingroup\$ Ok, that's mindblowing \$\endgroup\$ Nov 12, 2022 at 14:14
1
\$\begingroup\$

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!

\$\endgroup\$
2
  • \$\begingroup\$ Your code outputs "0s n=19 1111111111111111111n" which is not a near-repdigit. \$\endgroup\$
    – user108721
    Nov 13, 2022 at 19:45
  • \$\begingroup\$ @graffe Oh, of course. Now fixed. \$\endgroup\$
    – Arnauld
    Nov 13, 2022 at 20:36
1
\$\begingroup\$

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.

\$\endgroup\$
1
\$\begingroup\$

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

\$\endgroup\$
1
\$\begingroup\$

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.

\$\endgroup\$
5
  • 1
    \$\begingroup\$ Would it help to test for a last digit of 5 as well? \$\endgroup\$
    – Neil
    Nov 12, 2022 at 19:18
  • \$\begingroup\$ Can you add a TIO link? \$\endgroup\$
    – user108721
    Nov 12, 2022 at 19:34
  • \$\begingroup\$ @Neil it probably would. I'll add that in my next revision \$\endgroup\$
    – Seggan
    Nov 14, 2022 at 2:16
  • \$\begingroup\$ @Neil it actually makes it worse for some reason \$\endgroup\$
    – Seggan
    Nov 15, 2022 at 1:20
  • \$\begingroup\$ your profile statement is very true for this question :-) \$\endgroup\$
    – jdt
    Nov 17, 2022 at 0:35
1
\$\begingroup\$

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!

\$\endgroup\$
9
  • 1
    \$\begingroup\$ You are winning now! \$\endgroup\$
    – user108721
    Nov 16, 2022 at 16:38
  • \$\begingroup\$ @graffe Cool. What do you get on your machine if you run it with const auto tc = 1;? \$\endgroup\$
    – jdt
    Nov 16, 2022 at 17:32
  • 1
    \$\begingroup\$ Gets to 906 which is pretty good already. \$\endgroup\$
    – user108721
    Nov 16, 2022 at 17:42
  • 1
    \$\begingroup\$ that's probably because the Baille-PSW test filters most of the non-primes, so it won't even reach the Miller-Rabin test \$\endgroup\$
    – c--
    Nov 16, 2022 at 21:57
  • 2
    \$\begingroup\$ @c-- Say what? Wasting valuable CPU time freeing data ;-) I suppose it could affect the usage of the CPU cache and eventually, it will start paging data from the disk or run out of memory completely. \$\endgroup\$
    – jdt
    Nov 17, 2022 at 19:22
1
\$\begingroup\$

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.

\$\endgroup\$
5
  • \$\begingroup\$ This code is currently broken. It outputs numbers which are not near repdigits. \$\endgroup\$
    – user108721
    Nov 12, 2022 at 12:02
  • \$\begingroup\$ It should be fixed now. \$\endgroup\$ Nov 12, 2022 at 12:19
  • \$\begingroup\$ To be clear, the time I gave is for doing all n up to and including 9. \$\endgroup\$
    – user108721
    Nov 12, 2022 at 14:01
  • 3
    \$\begingroup\$ this is broken, the i += 6 should be replaced with the step argument to the range \$\endgroup\$
    – c--
    Nov 12, 2022 at 15:18
  • 2
    \$\begingroup\$ i was overwritten with the next value in the range after every iteration so the line i += 6 had no effect. Which meant this was testing divisibility by every number in the range(n, int(n**.5)), Try it Online! \$\endgroup\$
    – c--
    Nov 13, 2022 at 15:04
1
\$\begingroup\$

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

\$\endgroup\$
6
  • \$\begingroup\$ Made it to 1027. Any idea why jdt's code is so much faster? \$\endgroup\$
    – user108721
    Nov 20, 2022 at 12:03
  • \$\begingroup\$ @graffe, not really. I get slightly better performace with this version (1013 vs. 1015). Is the score consistent across runs? \$\endgroup\$
    – c--
    Nov 20, 2022 at 15:06
  • \$\begingroup\$ I always add -march=native too. People often seem to ignore that important option. \$\endgroup\$
    – user108721
    Nov 20, 2022 at 21:41
  • \$\begingroup\$ @graffe I guess it could be because he's using async IO, could you try the new version to see if there's a difference? \$\endgroup\$
    – c--
    Dec 6, 2022 at 15:33
  • \$\begingroup\$ I can do that on Friday. Thanks \$\endgroup\$
    – user108721
    Dec 6, 2022 at 19:35
1
\$\begingroup\$

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.")
\$\endgroup\$
0
\$\begingroup\$

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.

\$\endgroup\$
0
0
\$\begingroup\$

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

\$\endgroup\$
6
  • 1
    \$\begingroup\$ i love std::format. The whole std::cout thing was crappy from the start. \$\endgroup\$
    – jdt
    Nov 14, 2022 at 3:20
  • \$\begingroup\$ yeah, it's too bad compilers don't support std::print yet, but {fmt} is really good. \$\endgroup\$
    – c--
    Nov 14, 2022 at 13:37
  • \$\begingroup\$ MSVC 2019 does support it with /std:c++latest. My laptop is pretty much useless for benchmarking, do you get a noticeable improvement with your code vs mine? \$\endgroup\$
    – jdt
    Nov 14, 2022 at 13:52
  • \$\begingroup\$ @jdt not really, I just got 593 with yours vs 588 with mine in 5 minutes. I think it could be that mp::pow is not a good way of generating ones or the fact that I'm searching for the non-repdigit from the end (which could affect the rate at which they're found if for these values of n they're easier to find from the start). \$\endgroup\$
    – c--
    Nov 14, 2022 at 14:35
  • \$\begingroup\$ I think it's more likely the temperature of the CPU. On my laptop, it starts at 120% cpu utilization but drops to about 80% after a while. Turning down the AC or putting your computer in a freezer might make the biggest difference :-) \$\endgroup\$
    – jdt
    Nov 14, 2022 at 14:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.