Least m such that n! + m is prime

From Rosetta Code
Least m such that n! + m is prime is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.

Find the minimum positive integer m such that n factorial plus m is prime.


E.G.
   0! = 1. The next prime greater than 1 is 2. 2 - 1 = 1, so a(0) = 1.
   1! = 1. The next prime greater than 1 is 2. 2 - 1 = 1, so a(1) = 1.
   2! = 2. The next prime greater than 2 is 3. 3 - 2 = 1, so a(2) = 1.
   3! = 6. The next prime greater than 6 is 7. 7 - 6 = 1, so a(3) = 1.
   4! = 24. The next prime greater than 24 is 29. 29 - 24 = 5, so a(4) = 5.

and so on...


Task
  • Find and display the first fifty terms in the series. (0! through 49!)
  • Find and display the position and value of the first m greater than 1000.


Stretch
  • Find and display the position and value of each the first m greater than 2000, 3000, 4000 ... 10,000.


See also


ALGOL 68

Works with: ALGOL 68G version Any - tested with release 2.8.3.win32

Uses ALGOL 68G's LONG LONG INT which has programmer specifiable precision, 500 digits is sufficient for the basic task (the Millar Rabin routine needs extra digits, even though 107! has around 170 digits).

BEGIN # find the least m such that n! + m is prime for various n              #
    PR precision 500 PR                  # set the precision of LONG LOMG INT #
    PR read "primes.incl.a68" PR                    # include priee utilities #
    LONG LONG INT factorial n := 1;
    INT m := 0;
    FOR n FROM 0 WHILE m < 1000 DO
        IF n > 0 THEN
            factorial n *:= n
        FI;
        m := 1;
        WHILE NOT is probably prime( factorial n + m ) DO
            m +:= 2
        OD;
        IF n < 50 THEN
            print( ( " ", whole( m, -4 ) ) );
            IF ( n + 1 ) MOD 10 = 0 THEN print( ( newline ) ) FI
        ELIF m > 1000 THEN
            print( ( "First m > 1000: ", whole( m, 0 ), " for ", whole( n, 0 ), "!", newline ) )
        FI
    OD
END
Output:
    1    1    1    1    5    7    7   11   23   17
   11    1   29   67   19   43   23   31   37   89
   29   31   31   97  131   41   59    1   67  223
  107  127   79   37   97   61  131    1   43   97
   53    1   97   71   47  239  101  233   53   83
First m > 1000: 1069 for 107!

C

Translation of: Wren
Library: GMP
#include <stdio.h>
#include <gmp.h>
#include <locale.h>

#define LIMIT 10000

int main() {
    mpz_t fact, p;
    mpz_init_set_ui(fact, 1);
    mpz_init(p);
    int i, diffs[50], t = 1000;
    unsigned long n, m;
    for (n = 0; ; ++n) {
        if (n > 0) mpz_mul_ui(fact, fact, n);
        mpz_nextprime(p, fact);
        mpz_sub(p, p, fact);
        m = mpz_get_ui(p);
        setlocale(LC_NUMERIC, "");
        if (n < 50) diffs[n] = m;
        if (n == 49) {
            printf("Least positive m such that n! + m is prime; first 50:\n");
            for (i = 0; i < 50; ++i) {
                printf("%3d  ", diffs[i]);
                if (!((i+1)%10)) printf("\n");
            }
            printf("\n");
        } else if (m > t) {
            do {
                printf("First m > %'6d is %'6ld at position %ld\n", t, m, n);
                t += 1000;
            } while (m > t);
            if (t > LIMIT) break;
        }
    }
    mpz_clear(fact);
    mpz_clear(p);
    return 0;
}
Output:
Same as Wren example.

J

   (4&p:-])!i.5 10x
  1   1  1  1   5   7   7  11 23  17
 11   1 29 67  19  43  23  31 37  89
 29  31 31 97 131  41  59   1 67 223
107 127 79 37  97  61 131   1 43  97
 53   1 97 71  47 239 101 233 53  83
   1 i.~1000 < (4&p:-])!i.200x
107

Julia

Translation of: Wren
""" rosettacode.orgwiki/Least_m_such_that_n!_%2B_m_is_prime """

using Primes

function least_m_fact_to_prime(number_to_print, delta_limit)
    fact, p, m, n, t = big"1", big"0", big"0", 0, 1000
    diffs = zeros(BigInt, number_to_print)
    while true
        if n > 0
            fact *= n
            p = nextprime(fact + 1)
            m = p - fact
            if n < number_to_print
                diffs[n] = m
            end
            if n == number_to_print - 1
                println("Least positive m such that n! + m is prime; first $number_to_print:")
                for (i, k) in enumerate(diffs)
                    print(lpad(k, 5), i % 10 == 0 ? "\n" : "")
                end
            elseif m > t
                while true
                    print("\nFirst m > $t is $m at position $n.")
                    t += 1000
                    if m <= t
                        break
                    end
                end
                if t > delta_limit
                    return
                end
            end
        end
        n += 1
    end
end

least_m_fact_to_prime(50, 10_000)
Output:
Least positive m such that n! + m is prime; first 50:
    1    1    1    5    7    7   11   23   17   11
    1   29   67   19   43   23   31   37   89   29
   31   31   97  131   41   59    1   67  223  107
  127   79   37   97   61  131    1   43   97   53
    1   97   71   47  239  101  233   53   83    0

First m > 1000 is 1069 at position 107.
First m > 2000 is 3391 at position 192.
First m > 3000 is 3391 at position 192.
First m > 4000 is 4943 at position 284.
First m > 5000 is 5233 at position 384.
First m > 6000 is 6131 at position 388.
First m > 7000 is 9067 at position 445.
First m > 8000 is 9067 at position 445.
First m > 9000 is 9067 at position 445.
First m > 10000 is 12619 at position 599.
First m > 11000 is 12619 at position 599.
First m > 12000 is 12619 at position 599.

Nim

import std/strformat
import integers

const Lim = 10000
const Step = 1000

var n = 0
var lim = 1000
var f = newInteger(1)
echo "Least positive m such that n! + m is prime; first 50:"
while true:
  var m = nextPrime(f) - f
  if n < 50:
    stdout.write &"{m:3}"
    stdout.write if n mod 10 == 9: '\n' else: ' '
    if n == 49: echo()
  else:
    while m > lim:
      echo &"First m > {lim:5} is {m:5} at position {n}."
      inc lim, Step
    if lim > Lim: break
  inc n
  f *= n
Output:
Least positive m such that n! + m is prime; first 50:
  1   1   1   1   5   7   7  11  23  17
 11   1  29  67  19  43  23  31  37  89
 29  31  31  97 131  41  59   1  67 223
107 127  79  37  97  61 131   1  43  97
 53   1  97  71  47 239 101 233  53  83

First m >  1000 is  1069 at position 107.
First m >  2000 is  3391 at position 192.
First m >  3000 is  3391 at position 192.
First m >  4000 is  4943 at position 284.
First m >  5000 is  5233 at position 384.
First m >  6000 is  6131 at position 388.
First m >  7000 is  9067 at position 445.
First m >  8000 is  9067 at position 445.
First m >  9000 is  9067 at position 445.
First m > 10000 is 12619 at position 599.
First m > 11000 is 12619 at position 599.
First m > 12000 is 12619 at position 599.

Perl

Library: ntheory
use strict;
use warnings;
use ntheory qw<next_prime factorial vecfirstidx>;

my($n,@least_m) = 0;
do {
    my $f = factorial($n++);
    push @least_m, next_prime($f) - $f;
} until $least_m[-1] > 10000;

print "Least positive m such that n! + m is prime; first fifty:\n";
print sprintf(('%4d')x50, @least_m[0..49]) =~ s/.{40}\K(?=.)/\n/gr . "\n\n";

for my $n (map { 1000 * $_ } 1..10) {
    my $key = vecfirstidx { $_ > $n } @least_m;
    printf "First m > $n is %d at position %d\n", $least_m[$key], $key;
}
Output:
Least positive m such that n! + m is prime; first fifty:
   1   1   1   1   5   7   7  11  23  17
  11   1  29  67  19  43  23  31  37  89
  29  31  31  97 131  41  59   1  67 223
 107 127  79  37  97  61 131   1  43  97
  53   1  97  71  47 239 101 233  53  83

First m > 1000 is 1069 at position 107
First m > 2000 is 3391 at position 192
First m > 3000 is 3391 at position 192
First m > 4000 is 4943 at position 284
First m > 5000 is 5233 at position 384
First m > 6000 is 6131 at position 388
First m > 7000 is 9067 at position 445
First m > 8000 is 9067 at position 445
First m > 9000 is 9067 at position 445
First m > 10000 is 12619 at position 599

Phix

Translation of: C
with javascript_semantics
atom t0 = time()
requires("1.0.3") -- mpz_nextprime() added 
constant LIMIT = iff(platform()=JS?1000:10000)
include mpfr.e
mpz {fact, p} = mpz_inits(2,1)
sequence diffs = {}
integer n=0, m, t = 1000
while t<=LIMIT do
    progress("position: %d\r",{n})
    if n>0 then mpz_mul_si(fact, fact, n) end if
    mpz_nextprime(p, fact)
    mpz_sub(p, p, fact);
    m = mpz_get_integer(p);
    if length(diffs)<50 then
        diffs &= m
        if length(diffs)=50 then
            printf(1,"Least positive m such that n! + m is prime; first 50:\n%s\n",
                     {join_by(diffs,1,10,"  ", fmt:="%3d")})
        end if
    elsif m>t then
        string e = elapsed(time()-t0,0," (%s)")
        do
            printf(1,"First m > %,6d is %,6d at position %,d%s\n", {t, m, n, e})
            e = ""
            t += 1000
        until t>m
    end if
    n += 1
end while
?elapsed(time()-t0)
Output:
Least positive m such that n! + m is prime; first 50:
  1    1    1    1    5    7    7   11   23   17
 11    1   29   67   19   43   23   31   37   89
 29   31   31   97  131   41   59    1   67  223
107  127   79   37   97   61  131    1   43   97
 53    1   97   71   47  239  101  233   53   83

First m >  1,000 is  1,069 at position 107 (0.2s)
First m >  2,000 is  3,391 at position 192 (3.4s)
First m >  3,000 is  3,391 at position 192
First m >  4,000 is  4,943 at position 284 (27.3s)
First m >  5,000 is  5,233 at position 384 (2 minutes and 15s)
First m >  6,000 is  6,131 at position 388
First m >  7,000 is  9,067 at position 445 (4 minutes and 56s)
First m >  8,000 is  9,067 at position 445
First m >  9,000 is  9,067 at position 445
First m > 10,000 is 12,619 at position 599 (26 minutes and 15s)
First m > 11,000 is 12,619 at position 599
First m > 12,000 is 12,619 at position 599
"26 minutes and 15s"

For comparison, the Julia entry above took 18 mins 38s on the same box, and Perl an even more impressive 10 mins 44s.

Quackery

from, index, and end are defined at Loops/Increment loop index within loop body#Quackery.

prime is defined at Miller–Rabin primality test#Quackery.

  [ 1 swap times [ i^ 1+ *  ] ] is ! ( n --> n )

  [] 50 times
   [ i^ !
     1+ from
       [ index prime if
           [ index i^ ! -
             join end ] ] ]
  [] swap
  witheach
    [ number$ nested join ]
  49 wrap$
Output:
1 1 1 1 5 7 7 11 23 17 11 1 29 67 19 43 23 31 37
89 29 31 31 97 131 41 59 1 67 223 107 127 79 37
97 61 131 1 43 97 53 1 97 71 47 239 101 233 53 83

Raku

my @f = lazy flat 1, [\×] 1..*;

sink @f[700]; # pre-reify for concurrency

my @least-m = lazy (^∞).hyper(:2batch).map: {(1..*).first: -> \n {(@f[$_] + n).is-prime}};

say "Least positive m such that n! + m is prime; first fifty:\n" 
 ~ @least-m[^50].batch(10)».fmt("%3d").join: "\n";

for (1..10).map: * × 1e3 {
    my $key = @least-m.first: * > $_, :k;
    printf "\nFirst m > $_ is %d at position %d\n", @least-m[$key], $key;
}
Output:
Least positive m such that n! + m is prime; first fifty:
  1   1   1   1   5   7   7  11  23  17
 11   1  29  67  19  43  23  31  37  89
 29  31  31  97 131  41  59   1  67 223
107 127  79  37  97  61 131   1  43  97
 53   1  97  71  47 239 101 233  53  83

First m > 1000 is 1069 at position 107

First m > 2000 is 3391 at position 192

First m > 3000 is 3391 at position 192

First m > 4000 is 4943 at position 284

First m > 5000 is 5233 at position 384

First m > 6000 is 6131 at position 388

First m > 7000 is 9067 at position 445

First m > 8000 is 9067 at position 445

First m > 9000 is 9067 at position 445

First m > 10000 is 12619 at position 599

RPL

« 49 0 « ←n FACT DUP NEXTPRIME SWAP - »
  → max ←n m
  « m '←n' 0 max 1 SEQ 
    max '←n' STO
    DO '←n' INCR DROP m EVAL
    UNTIL 1000 > END 
    ←n
» » 'TASK' STO
Output:
2: { 1 1 1 1 5 7 7 11 23 17 11 1 29 67 19 43 23 31 37 89 29 31 31 97 131 41 59 1 67 223 107 127 79 37 97 61 131 1 43 97 53 1 97 71 47 239 101 233 53 83 }
1: 107

Needs over an hour on an iOS emulator (iHP48).

Ruby

require 'openssl'

def next_prime(n) = ((n+1)..).detect{|n| OpenSSL::BN.new(n).prime?}
def fact(n) = (1..n).inject(:*) || 1

enum_diffs = (0..).lazy.map do |n|
  f = fact(n)
  next_prime(f) - f
end

enum_diffs.first(50).each_slice(10){|s| puts "%4d"*s.size % s}
puts "\nFirst m > 1000 is %d at position %d." % enum_diffs.with_index.detect{|d,_id| d>1000}
Output:
   1   1   1   1   5   7   7  11  23  17
  11   1  29  67  19  43  23  31  37  89
  29  31  31  97 131  41  59   1  67 223
 107 127  79  37  97  61 131   1  43  97
  53   1  97  71  47 239 101 233  53  83

First m > 1000 is 1069 at position 107.

Sidef

with (50) {|N|
    say "Least positive m such that n! + m is prime (first #{N}):"
    ^N -> map {|n|
        var f = n!; 1..Inf -> first {|k| f+k -> is_prime }
    }.each_slice(10, {|*s|
        say s.map{ '%3s' % _ }.join(' ')
    })
}

say ''; var prev = 0
for n in (1..5 -> map { 1e3*_ }) {
    var m = (prev..Inf -> lazy.map{|k|
        var f = k!; [k, f.next_prime - f]
    }.first {|k|
        k.tail >= n
    })
    say "First m > #{n} is #{m.tail} at position #{m.head}"
    prev = m.head
}
Output:
Least positive m such that n! + m is prime (first 50):
  1   1   1   1   5   7   7  11  23  17
 11   1  29  67  19  43  23  31  37  89
 29  31  31  97 131  41  59   1  67 223
107 127  79  37  97  61 131   1  43  97
 53   1  97  71  47 239 101 233  53  83

First m > 1000 is 1069 at position 107
First m > 2000 is 3391 at position 192
First m > 3000 is 3391 at position 192
First m > 4000 is 4943 at position 284
First m > 5000 is 5233 at position 384

Wren

Library: Wren-gmp
Library: Wren-fmt
import "./gmp" for Mpz
import "./fmt" for Fmt

var fact = Mpz.one
var p = Mpz.new()
var diffs = List.filled(50, 0)
var n = 0
var t = 1000
var limit = 10000
while (true) {
    if (n > 0) fact.mul(n)
    p.nextPrime(fact)
    var m = p.sub(fact).toNum
    if (n < 50) diffs[n] = m
    if (n == 49) {
        System.print("Least positive m such that n! + m is prime; first 50:")
        Fmt.tprint("$3d ", diffs, 10)
        System.print()
    } else if (m > t) {
        while (true) {
            Fmt.print("First m > $,6d is $,6d at position $d", t, m, n)
            t = t + 1000
            if (m <= t) break
        }
        if (t > limit) return
    }
    n = n + 1
}
Output:
Least positive m such that n! + m is prime; first 50:
  1    1    1    1    5    7    7   11   23   17 
 11    1   29   67   19   43   23   31   37   89 
 29   31   31   97  131   41   59    1   67  223 
107  127   79   37   97   61  131    1   43   97 
 53    1   97   71   47  239  101  233   53   83 

First m >  1,000 is  1,069 at position 107
First m >  2,000 is  3,391 at position 192
First m >  3,000 is  3,391 at position 192
First m >  4,000 is  4,943 at position 284
First m >  5,000 is  5,233 at position 384
First m >  6,000 is  6,131 at position 388
First m >  7,000 is  9,067 at position 445
First m >  8,000 is  9,067 at position 445
First m >  9,000 is  9,067 at position 445
First m > 10,000 is 12,619 at position 599
First m > 11,000 is 12,619 at position 599
First m > 12,000 is 12,619 at position 599