Erdős-primes: Difference between revisions
Line 594: | Line 594: | ||
Implementation: |
Implementation: |
||
<lang J>NB. erdos primes greater than !k and less than or equal to !k+1 |
<lang J>NB. erdos primes greater than !k and less than or equal to !k+1 (where !k is the factorial of k) |
||
erdospr=: {{ k=. y |
erdospr=: {{ k=. y |
||
f=. !k+0 1 |
f=. !k+0 1 |
Revision as of 04:44, 4 July 2022
- Definitions
In mathematics, Erdős primes are prime numbers such that all p-k! for 1<=k!<p are composite.
- Task
Write a program to determine (and show here) all Erdős primes less than 2500.
Optionally, show the number of Erdős primes.
- Stretch goal
Show that the 7,875th Erdős prime is 999,721 (the highest below 1,000,000)
- Also see
-
- the OEIS entry: A064152 Erdos primes.
11l
<lang 11l>F primes_upto(limit)
V is_prime = [0B] * 2 [+] [1B] * (limit - 1) L(n) 0 .< Int(limit ^ 0.5 + 1.5) I is_prime[n] L(i) (n * n .. limit).step(n) is_prime[i] = 0B R is_prime
V is_prime = primes_upto(1'000'000) V primeList = enumerate(is_prime).filter((i, prime) -> prime).map((i, prime) -> i)
[Int] factorials L(n) 1..
I factorial(n) >= 1'000'000 L.break factorials.append(factorial(n))
F isErdosPrime(p)
L(f) :factorials I f >= p L.break I :is_prime[p - f] R 0B R 1B
[Int] erdosList2500 L(p) primeList
I p >= 2500 L.break I isErdosPrime(p) erdosList2500.append(p)
print(‘Found ’erdosList2500.len‘ Erdos primes less than 2500:’) L(prime) erdosList2500
print(‘#5’.format(prime), end' I (L.index + 1) % 10 == 0 {"\n"} E ‘ ’)
print()
V count = 0 L(p) primeList
I isErdosPrime(p) count++ I count == 7875 print("\nThe 7875th Erdos prime is "p‘.’) L.break</lang>
- Output:
Found 25 Erdos primes less than 2500: 2 101 211 367 409 419 461 557 673 709 769 937 967 1009 1201 1259 1709 1831 1889 2141 2221 2309 2351 2411 2437 The 7875th Erdos prime is 999721.
Action!
<lang Action!>INCLUDE "H6:SIEVE.ACT"
BYTE Func IsErdosPrime(INT x BYTE ARRAY primes)
INT k,f
IF primes(x)=0 THEN RETURN (0) FI
k=1 f=1 WHILE f<x DO IF primes(x-f) THEN RETURN (0) FI k==+1 f==*k OD
RETURN (1)
PROC Main()
DEFINE MAX="2499" BYTE ARRAY primes(MAX+1) INT i,count=[0]
Put(125) PutE() ;clear the screen Sieve(primes,MAX+1) FOR i=2 TO MAX DO IF IsErdosPrime(i,primes) THEN PrintI(i) Put(32) count==+1 FI OD PrintF("%E%EThere are %I Erdos primes",count)
RETURN</lang>
- Output:
Screenshot from Atari 8-bit computer
2 101 211 367 409 419 461 557 673 709 769 937 967 1009 1201 1259 1709 1831 1889 2141 2221 2309 2351 2411 2437 There are 25 Erdos primes
ALGOL 68
<lang algol68>BEGIN # find Erdős primes - primes p such that p-k! is composite for all 1<=k!<p #
# returns TRUE if p is an Erdős prime # PROC is erdos prime = ( INT p )BOOL: IF NOT prime[ p ] THEN FALSE ELSE BOOL result := TRUE; FOR k WHILE factorial[ k ] < p AND result DO result := NOT prime[ p - factorial[ k ] ] OD; result FI # is erdos prime # ; INT max prime = 2500; # maximum number we will consider # # construct a table of factorials large enough for max prime # [ 1 : 12 ]INT factorial; factorial[ 1 ] := 1; FOR f FROM 2 TO UPB factorial DO factorial[ f ] := factorial[ f - 1 ] * f OD; # sieve the primes to max prime # PR read "primes.incl.a68" PR []BOOL prime = PRIMESIEVE max prime; # find the Erdős primes # INT e count := 0; IF is erdos prime( 2 ) THEN print( ( " 2" ) ); e count +:= 1 FI; FOR p FROM 3 BY 2 TO UPB prime DO IF is erdos prime( p ) THEN print( ( " ", whole( p, 0 ) ) ); e count +:= 1 FI OD; print( ( newline, "Found ", whole( e count, 0 ), " Erdos primes" ) )
END</lang>
- Output:
2 101 211 367 409 419 461 557 673 709 769 937 967 1009 1201 1259 1709 1831 1889 2141 2221 2309 2351 2411 2437 Found 25 Erdos primes
Arturo
<lang rebol>factorials: map 1..20 => [product 1..&] erdos?: function [x][
if not? prime? x -> return false
loop factorials 'f [ if f >= x -> break if prime? x - f -> return false ] return true
]
loop split.every:10 select 2..2500 => erdos? 'a ->
print map a => [pad to :string & 5]</lang>
- Output:
2 101 211 367 409 419 461 557 673 709 769 937 967 1009 1201 1259 1709 1831 1889 2141 2221 2309 2351 2411 2437
AWK
<lang AWK>
- syntax: GAWK -f ERDOS-PRIMES.AWK
- converted from FreeBASIC
BEGIN {
while (++i) { if (is_erdos_prime(i)) { if (i < 2500) { printf("%d ",i) count1++ } if (++count2 == 7875) { printf("\nErdos primes 1-2500: %d\nErdos prime %d: %d\n",count1,count2,i) break } } } exit(0)
} function is_erdos_prime(p, kf,m) {
if (!is_prime(p)) { return(0) } kf = m = 1 while (kf < p) { kf *= m++ if (is_prime(p-kf)) { return(0) } } return(1)
} function is_prime(x, i) {
if (x <= 1) { return(0) } for (i=2; i<=int(sqrt(x)); i++) { if (x % i == 0) { return(0) } } return(1)
} </lang>
- Output:
2 101 211 367 409 419 461 557 673 709 769 937 967 1009 1201 1259 1709 1831 1889 2141 2221 2309 2351 2411 2437 Erdos primes 1-2500: 25 Erdos prime 7875: 999721
C#
<lang csharp>using System; using static System.Console; class Program {
const int lmt = (int)1e6, first = 2500; static int[] f = new int[10]; static void Main(string[] args) { f[0] = 1; for (int a = 0, b = 1; b < f.Length; a = b++) f[b] = f[a] * (b + 1); int pc = 0, nth = 0, lv = 0; for (int i = 2; i < lmt; i++) if (is_erdos_prime(i)) { if (i < first) Write("{0,5:n0}{1}", i, pc++ % 5 == 4 ? "\n" : " "); nth++; lv = i; } Write("\nCount of Erdős primes between 1 and {0:n0}: {1}\n{2} Erdős prime (the last one under {3:n0}): {4:n0}", first, pc, ord(nth), lmt, lv); }
static string ord(int n) { return string.Format("{0:n0}", n) + new string[]{"th", "st", "nd", "rd", "th", "th", "th", "th", "th", "th"}[n % 10]; }
static bool is_erdos_prime(int p) { if (!is_pr(p)) return false; int m = 0, t; while ((t = p - f[m++]) > 0) if (is_pr(t)) return false; return true; bool is_pr(int x) { if (x < 4) return x > 1; if ((x & 1) == 0) return false; for (int i = 3; i * i <= x; i += 2) if (x % i == 0) return false; return true; } } }</lang>
- Output:
2 101 211 367 409 419 461 557 673 709 769 937 967 1,009 1,201 1,259 1,709 1,831 1,889 2,141 2,221 2,309 2,351 2,411 2,437 Count of Erdős primes between 1 and 2,500: 25 7,875th Erdős prime (the last one under 1,000,000): 999,721
C++
<lang cpp>#include <cstdint>
- include <iomanip>
- include <iostream>
- include <set>
- include <primesieve.hpp>
class erdos_prime_generator { public:
erdos_prime_generator() {} uint64_t next();
private:
bool erdos(uint64_t p) const; primesieve::iterator iter_; std::set<uint64_t> primes_;
};
uint64_t erdos_prime_generator::next() {
uint64_t prime; for (;;) { prime = iter_.next_prime(); primes_.insert(prime); if (erdos(prime)) break; } return prime;
}
bool erdos_prime_generator::erdos(uint64_t p) const {
for (uint64_t k = 1, f = 1; f < p; ++k, f *= k) { if (primes_.find(p - f) != primes_.end()) return false; } return true;
}
int main() {
std::wcout.imbue(std::locale("")); erdos_prime_generator epgen; const int max_print = 2500; const int max_count = 7875; uint64_t p; std::wcout << L"Erd\x151s primes less than " << max_print << L":\n"; for (int count = 1; count <= max_count; ++count) { p = epgen.next(); if (p < max_print) std::wcout << std::setw(6) << p << (count % 10 == 0 ? '\n' : ' '); } std::wcout << L"\n\nThe " << max_count << L"th Erd\x151s prime is " << p << L".\n"; return 0;
}</lang>
- Output:
Erdős primes less than 2,500: 2 101 211 367 409 419 461 557 673 709 769 937 967 1,009 1,201 1,259 1,709 1,831 1,889 2,141 2,221 2,309 2,351 2,411 2,437 The 7,875th Erdős prime is 999,721.
F#
This task uses Extensible Prime Generator (F#) <lang fsharp> // Erdős Primes. Nigel Galloway: March 22nd., 2021 let rec fN g=function 1->g |n->fN(g*n)(n-1) let rec fG n g=seq{let i=fN 1 n in if i<g then yield (isPrime>>not)(g-i); yield! fG(n+1) g} let eP()=primes32()|>Seq.filter(fG 1>>Seq.forall id) eP()|>Seq.takeWhile((>)2500)|>Seq.iter(printf "%d "); printfn "\n\n7875th Erdős prime is %d" (eP()|>Seq.item 7874) </lang>
- Output:
2 101 211 367 409 419 461 557 673 709 769 937 967 1009 1201 1259 1709 1831 1889 2141 2221 2309 2351 2411 2437 7875th Erdos prime is 999721
Factor
<lang>USING: formatting io kernel lists lists.lazy math math.factorials math.primes math.primes.lists math.vectors prettyprint sequences tools.memory.private ;
- facts ( -- list ) 1 lfrom [ n! ] lmap-lazy ;
- n!< ( p -- seq ) [ facts ] dip [ < ] curry lwhile list>array ;
- erdős? ( p -- ? ) dup n!< n-v [ prime? ] none? ;
- erdős ( -- list ) lprimes [ erdős? ] lfilter ;
erdős [ 2500 < ] lwhile list>array dup length "Found %d Erdős primes < 2500:\n" printf [ bl ] [ pprint ] interleave nl nl
7874 erdős lnth commas "The 7,875th Erdős prime is %s.\n" printf</lang>
- Output:
Found 25 Erdős primes < 2500: 2 101 211 367 409 419 461 557 673 709 769 937 967 1009 1201 1259 1709 1831 1889 2141 2221 2309 2351 2411 2437 The 7,875th Erdős prime is 999,721.
Forth
<lang forth>: prime? ( n -- ? ) here + c@ 0= ;
- notprime! ( n -- ) here + 1 swap c! ;
- prime_sieve { n -- }
here n erase 0 notprime! 1 notprime! n 4 > if n 4 do i notprime! 2 +loop then 3 begin dup dup * n < while dup prime? if n over dup * do i notprime! dup 2* +loop then 2 + repeat drop ;
- erdos_prime? { p -- ? }
p prime? if 1 1 begin dup p < while p over - prime? if 2drop false exit then swap 1+ swap over * repeat 2drop true else false then ;
- print_erdos_primes { n -- }
." Erdos primes < " n 1 .r ." :" cr n prime_sieve 0 n 0 do i erdos_prime? if i 5 .r 1+ dup 10 mod 0= if cr then then loop cr ." Count: " . cr ;
2500 print_erdos_primes bye</lang>
- Output:
Erdos primes < 2500: 2 101 211 367 409 419 461 557 673 709 769 937 967 1009 1201 1259 1709 1831 1889 2141 2221 2309 2351 2411 2437 Count: 25
FreeBASIC
I won't bother reproducing a primality-testing function; use the one from Primality_by_trial_division#FreeBASIC. <lang freebasic>#include "isprime.bas"
function is_erdos_prime( p as uinteger ) as boolean
if not isprime(p) then return false dim as uinteger kf=1, m=1 while kf < p kf*=m : m+=1 if isprime(p - kf) then return false wend return true
end function
dim as integer c = 0, i = 1 while c<7875
i+=1 if is_erdos_prime(i) then c+=1 if i<2500 or c=7875 then print c, i end if
wend</lang>
- Output:
1 22 101 3 211 4 367 5 409 6 419 7 461 8 557 9 673 10 709 11 769 12 937 13 967 14 1009 15 1201 16 1259 17 1709 18 1831 19 1889 20 2141 21 2221 22 2309 23 2351 24 2411 25 2437 7875 999721
Go
<lang go>package main
import "fmt"
func sieve(limit int) []bool {
limit++ // True denotes composite, false denotes prime. c := make([]bool, limit) // all false by default c[0] = true c[1] = true for i := 4; i < limit; i += 2 { c[i] = true } p := 3 // Start from 3. for { p2 := p * p if p2 >= limit { break } for i := p2; i < limit; i += 2 * p { c[i] = true } for { p += 2 if !c[p] { break } } } return c
}
func commatize(n int) string {
s := fmt.Sprintf("%d", n) if n < 0 { s = s[1:] } le := len(s) for i := le - 3; i >= 1; i -= 3 { s = s[0:i] + "," + s[i:] } if n >= 0 { return s } return "-" + s
}
func main() {
limit := int(1e6) c := sieve(limit - 1) var erdos []int for i := 2; i < limit; { if !c[i] { found := true for j, fact := 1, 1; fact < i; { if !c[i-fact] { found = false break } j++ fact = fact * j } if found { erdos = append(erdos, i) } } if i > 2 { i += 2 } else { i += 1 } } lowerLimit := 2500 var erdosLower []int for _, e := range erdos { if e < lowerLimit { erdosLower = append(erdosLower, e) } else { break } } fmt.Printf("The %d Erdős primes under %s are\n", len(erdosLower), commatize(lowerLimit)) for i, e := range erdosLower { fmt.Printf("%6d", e) if (i+1)%10 == 0 { fmt.Println() } } show := 7875 fmt.Printf("\n\nThe %s Erdős prime is %s.\n", commatize(show), commatize(erdos[show-1]))
}</lang>
- Output:
The 25 Erdős primes under 2,500 are 2 101 211 367 409 419 461 557 673 709 769 937 967 1009 1201 1259 1709 1831 1889 2141 2221 2309 2351 2411 2437 The 7,875 Erdős prime is 999,721.
J
Implementation:
<lang J>NB. erdos primes greater than !k and less than or equal to !k+1 (where !k is the factorial of k) erdospr=: {{ k=. y
f=. !k+0 1 p=. (#~ 1= f&I.) p:(+i.)/0 1+p:inv f p#~*/|:0=1 p:p-/!i.1+k
}}
NB. erdos primes less than j erdosprs=: Template:(</lang>
Task examples:
<lang J> erdosprs 2500 2 101 211 367 409 419 461 557 673 709 769 937 967 1009 1201 1259 1709 1831 1889 2141 2221 2309 2351 2411 2437
(#,{:) erdosprs 1e6
7875 999721</lang>
Java
<lang java>import java.util.*;
public class ErdosPrimes {
public static void main(String[] args) { boolean[] sieve = primeSieve(1000000); int maxPrint = 2500; int maxCount = 7875; System.out.printf("Erd\u0151s primes less than %d:\n", maxPrint); for (int count = 0, prime = 1; count < maxCount; ++prime) { if (erdos(sieve, prime)) { ++count; if (prime < maxPrint) { System.out.printf("%6d", prime); if (count % 10 == 0) System.out.println(); } if (count == maxCount) System.out.printf("\n\nThe %dth Erd\u0151s prime is %d.\n", maxCount, prime); } } }
private static boolean erdos(boolean[] sieve, int p) { if (!sieve[p]) return false; for (int k = 1, f = 1; f < p; ++k, f *= k) { if (sieve[p - f]) return false; } return true; }
private static boolean[] primeSieve(int limit) { boolean[] sieve = new boolean[limit]; Arrays.fill(sieve, true); if (limit > 0) sieve[0] = false; if (limit > 1) sieve[1] = false; for (int i = 4; i < limit; i += 2) sieve[i] = false; for (int p = 3; ; p += 2) { int q = p * p; if (q >= limit) break; if (sieve[p]) { int inc = 2 * p; for (; q < limit; q += inc) sieve[q] = false; } } return sieve; }
}</lang>
- Output:
Erdős primes less than 2500: 2 101 211 367 409 419 461 557 673 709 769 937 967 1009 1201 1259 1709 1831 1889 2141 2221 2309 2351 2411 2437 The 7875th Erdős prime is 999721.
jq
Works with gojq, the Go implementation of jq (but the second task requires an unreasonable amount of memory)
Preliminaries
<lang jq>def emit_until(cond; stream): label $out | stream | if cond then break $out else . end;
def is_prime:
. as $n | if ($n < 2) then false elif ($n % 2 == 0) then $n == 2 elif ($n % 3 == 0) then $n == 3 elif ($n % 5 == 0) then $n == 5 elif ($n % 7 == 0) then $n == 7 elif ($n % 11 == 0) then $n == 11 elif ($n % 13 == 0) then $n == 13 elif ($n % 17 == 0) then $n == 17 elif ($n % 19 == 0) then $n == 19 else 23 | until( (. * .) > $n or ($n % . == 0); . + 2) | . * . > $n end;
</lang> Erdős-primes <lang jq> def is_Erdos:
. as $p | if is_prime|not then false else label $out | foreach range(1; .+1) as $k (1; . * $k; if . >= $p then true, break $out elif ($p - .) | is_prime then 0, break $out
else empty end) // true
| . == true end ;
- emit the Erdos primes
def Erdos: range(2; infinite) | select(is_Erdos); </lang> The tasks <lang jq>"The Erdős primes less than 2500 are:", emit_until(. >= 2500; Erdos),
"\nThe 7875th Erdős prime is \(nth(7874; Erdos))." </lang>
- Output:
The Erdős primes less than 2500 are: 2 101 211 367 409 419 461 557 673 709 769 937 967 1009 1201 1259 1709 1831 1889 2141 2221 2309 2351 2411 2437 The 7875th Erdős prime is 999721.
Julia
<lang julia>using Primes, Formatting
function isErdős(p::Integer)
isprime(p) || return false for i in 1:100 kfac = factorial(i) kfac >= p && break isprime(p - kfac) && return false end return true
end
const Erdőslist = filter(isErdős, 1:1000000) const E2500 = filter(x -> x < 2500, Erdőslist)
println(length(E2500), " Erdős primes < 2500: ", E2500) println("The 7875th Erdős prime is ", format(Erdőslist[7875], commas=true))
</lang>
- Output:
25 Erdős primes < 2500: [2, 101, 211, 367, 409, 419, 461, 557, 673, 709, 769, 937, 967, 1009, 1201, 1259, 1709, 1831, 1889, 2141, 2221, 2309, 2351, 2411, 2437] The 7875th Erdős prime is 999,721
Mathematica /Wolfram Language
<lang Mathematica>ClearAll[ErdosPrimeQ] ErdosPrimeQ[p_Integer] := Module[{k},
If[PrimeQ[p], k = 1; While[k! < p, If[PrimeQ[p - k!], Return[False]]; k++; ]; True , False ] ]
sel = Select[Range[2500], ErdosPrimeQ] Length[sel] sel = Select[Range[999999], ErdosPrimeQ]; {Length[sel], Last[sel]}</lang>
- Output:
{2, 101, 211, 367, 409, 419, 461, 557, 673, 709, 769, 937, 967, 1009, 1201, 1259, 1709, 1831, 1889, 2141, 2221, 2309, 2351, 2411, 2437} 25 {7875, 999721}
Nim
<lang Nim>import math, sets, strutils, sugar
const N = 1_000_000
- Sieve of Erathostenes.
var isComposite: array[2..N, bool] for n in 2..N:
let n2 = n * n if n2 > N: break if not isComposite[n]: for k in countup(n2, N, n): isComposite[k] = true
template isPrime(n: int): bool = n > 1 and not isComposite[n]
let primeList = collect(newSeq):
for n in 2..N: if n.isPrime: n
const Factorials = collect(newSeq):
for n in 1..20: if fac(n) >= N: break fac(n)
proc isErdösPrime(p: int): bool =
## Check if prime "p" is an Erdös prime. for f in Factorials: if f >= p: break if (p - f).isPrime: return false result = true
let erdösList2500 = collect(newSeq):
for p in primeList: if p >= 2500: break if p.isErdösPrime: p
echo "Found $# Erdös primes less than 2500:".format(erdösList2500.len) for i, prime in erdösList2500:
stdout.write ($prime).align(5) stdout.write if (i+1) mod 10 == 0: '\n' else: ' '
echo()
var erdös7875: int var count = 0 for p in primeList:
if p.isErdösPrime: inc count if count == 7875: erdös7875 = p break
echo "\nThe 7875th Erdös prime is $#.".format(erdös7875)</lang>
- Output:
Found 25 Erdös primes less than 2500: 2 101 211 367 409 419 461 557 673 709 769 937 967 1009 1201 1259 1709 1831 1889 2141 2221 2309 2351 2411 2437 The 7875th Erdös prime is 999721.
Perl
<lang perl>use strict; use warnings; use feature 'say'; use utf8; binmode(STDOUT, ':utf8'); use List::AllUtils 'before'; use ntheory qw<is_prime factorial>;
sub is_erdos {
my($n) = @_; my $k; return unless is_prime($n); while ($n > factorial($k++)) { return if is_prime $n-factorial($k) } 'True'
}
my(@Erdős,$n); do { push @Erdős, $n if is_erdos(++$n) } until $n >= 1e6;
say 'Erdős primes < ' . (my $max = 2500) . ':'; say join ' ', before { $_ > 2500 } @Erdős; say "\nErdős prime #" . @Erdős . ' is ' . $Erdős[-1];</lang>
- Output:
Erdős primes < 2500: 2 101 211 367 409 419 461 557 673 709 769 937 967 1009 1201 1259 1709 1831 1889 2141 2221 2309 2351 2411 2437 Erdős prime #7875 is 999721
Phix
atom t0 = time() sequence facts = {1} function erdos(integer p) while facts[$]<p do facts &= facts[$]*(length(facts)+1) end while for i=length(facts) to 1 by -1 do integer pmk = p-facts[i] if pmk>0 then if is_prime(pmk) then return false end if end if end for return true end function sequence res = filter(get_primes_le(2500),erdos) printf(1,"Found %d Erdos primes < 2,500:\n%s\n\n",{length(res),join(apply(res,sprint))}) res = filter(get_primes_le(1000000),erdos) integer l = length(res) printf(1,"The %,d%s Erdos prime is %,d (%s)\n",{l,ord(l),res[$],elapsed(time()-t0)})
- Output:
Found 25 Erdos primes < 2,500: 2 101 211 367 409 419 461 557 673 709 769 937 967 1009 1201 1259 1709 1831 1889 2141 2221 2309 2351 2411 2437 The 7,875th Erdos prime is 999,721 (1.2s)
Raku
<lang perl6>use Lingua::EN::Numbers;
my @factorial = 1, |[\*] 1..*; my @Erdős = ^Inf .grep: { .is-prime and none($_ «-« @factorial[^(@factorial.first: * > $_, :k)]).is-prime }
put 'Erdős primes < 2500:'; put @Erdős[^(@Erdős.first: * > 2500, :k)]»., put "\nThe 7,875th Erdős prime is: " ~ @Erdős[7874].,</lang>
- Output:
Erdős primes < 2500: 2 101 211 367 409 419 461 557 673 709 769 937 967 1,009 1,201 1,259 1,709 1,831 1,889 2,141 2,221 2,309 2,351 2,411 2,437 The 7,875th Erdős prime is: 999,721
REXX
<lang rexx>/*REXX program counts/displays the number of Erdos primes under a specified number N. */ parse arg n cols . /*get optional number of primes to find*/ if n== | n=="," then n= 2500 /*Not specified? Then assume default.*/ if cols== | cols=="," then cols= 10 /* " " " " " */ nn= n; n= abs(n) /*N<0: shows highest Erdos prime< │N│ */ call genP n /*generate all primes under N. */ w= 10 /*width of a number in any column. */ if cols>0 then say ' index │'center(" Erdos primes that are < " n, 1 + cols*(w+1) ) if cols>0 then say '───────┼'center("" , 1 + cols*(w+1), '─') call facts /*generate a table of needed factorials*/ Eprimes= 0; idx= 1 /*initialize # of additive primes & idx*/ $= /*a list of additive primes (so far). */
do j=1 for #; prime= @.j /* */ do k=1 until fact.k>j /*verify: J-K! for 1≤K!<J are composite*/ z= prime - fact.k /*subtract some factorial from a prime.*/ if !.z then iterate j /*Is Z is a prime? Then skip it. */ end /*j*/
Eprimes= Eprimes + 1; EprimeL= j /*bump the count of Erdos primes. */ if cols<0 then iterate /*Build the list (to be shown later)? */ c= commas(j) /*maybe add some commas to the number. */ $= $ right(c, max(w, length(c) ) ) /*add Erdos prime to list, allow big #.*/ if Eprimes//cols\==0 then iterate /*have we populated a line of output? */ say center(idx, 7)'│' substr($, 2); $= /*display what we have so far (cols). */ idx= idx + cols /*bump the index count for the output*/ end /*j*/
if $\== then say center(idx, 7)"│" substr($, 2) /*possible display residual output.*/ if cols>0 then say '───────┴'center("" , 1 + cols*(w+1), '─') say say 'found ' commas(Eprimes) " Erdos primes < " commas(n) say if nn<0 then say commas(EprimeL) ' is the ' commas(Eprimes)th(Eprimes) " Erdos prime." exit 0 /*stick a fork in it, we're all done. */ /*──────────────────────────────────────────────────────────────────────────────────────*/ commas: parse arg ?; do jc=length(?)-3 to 1 by -3; ?=insert(',', ?, jc); end; return ? facts: arg x; fact.=1; do x=2 until fact.x>1e9; p= x-1; fact.x= x*fact.p; end; return th: parse arg th; return word('th st nd rd', 1+(th//10) *(th//100%10\==1) *(th//10<4)) /*──────────────────────────────────────────────────────────────────────────────────────*/ genP: parse arg n; @.=.; @.1=2; @.2=3; @.3=5; @.4=7; @.5=11; @.6=13; @.7=17; #= 7
w= length(n); !.=0; !.2=1; !.3=1; !.5=1; !.7=1; !.11=1; !.13=1; !.17=1 do j=@.7+2 by 2 while j<n /*continue on with the next odd prime. */ parse var j -1 _ /*obtain the last digit of the J var.*/ if _ ==5 then iterate /*is this integer a multiple of five? */ if j // 3 ==0 then iterate /* " " " " " " three? */ /* [↓] divide by the primes. ___ */ do k=4 to # while k*k<=j /*divide J by other primes ≤ √ J */ if j//@.k == 0 then iterate j /*÷ by prev. prime? ¬prime ___ */ end /*k*/ /* [↑] only divide up to √ J */ #= # + 1; @.#= j; !.j= 1 /*bump prime count; assign prime & flag*/ end /*j*/; return</lang>
- output when using the default inputs:
index │ Erdos primes that are < 2500 ───────┼─────────────────────────────────────────────────────────────────────────────────────────────────────────────── 1 │ 2 101 211 367 409 419 461 557 673 709 11 │ 769 937 967 1,009 1,201 1,259 1,709 1,831 1,889 2,141 21 │ 2,221 2,309 2,351 2,411 2,437 ───────┴─────────────────────────────────────────────────────────────────────────────────────────────────────────────── found 25 Erdos primes < 2500
- output when using the inputs of: 1000000 0
found 7,875 Erdos primes < 1,000,000 999,721 is the 7,875th Erdos prime.
Rust
<lang rust>// [dependencies] // primal = "0.3"
use std::collections::HashSet;
fn erdos_primes() -> impl std::iter::Iterator<Item = usize> {
let mut primes = HashSet::new(); let mut all_primes = primal::Primes::all(); std::iter::from_fn(move || { 'all_primes: for p in all_primes.by_ref() { primes.insert(p); let mut k = 1; let mut f = 1; while f < p { if primes.contains(&(p - f)) { continue 'all_primes; } k += 1; f *= k; } return Some(p); } None })
}
fn main() {
let mut count = 0; println!("Erd\u{151}s primes less than 2500:"); for p in erdos_primes().take_while(|x| *x < 2500) { count += 1; if count % 10 == 0 { println!("{:4}", p); } else { print!("{:4} ", p); } } println!(); if let Some(p) = erdos_primes().nth(7874) { println!("\nThe 7875th Erd\u{151}s prime is {}.", p); }
}</lang>
- Output:
Erdős primes less than 2500: 2 101 211 367 409 419 461 557 673 709 769 937 967 1009 1201 1259 1709 1831 1889 2141 2221 2309 2351 2411 2437 The 7875th Erdős prime is 999721.
Sidef
<lang ruby>func is_erdos_prime(p) {
return true if p==2 return false if !p.is_prime
var f = 1
for (var k = 2; f < p; k++) { p - f -> is_composite || return false f *= k }
return true
}
say ("Erdős primes <= 2500: ", 1..2500 -> grep(is_erdos_prime)) say ("The 7875th Erdős prime is: ", is_erdos_prime.nth(7875))</lang>
- Output:
Erdős primes <= 2500: [2, 101, 211, 367, 409, 419, 461, 557, 673, 709, 769, 937, 967, 1009, 1201, 1259, 1709, 1831, 1889, 2141, 2221, 2309, 2351, 2411, 2437] The 7875th Erdős prime is: 999721
Tiny BASIC
Can't manage the stretch goal because integers are limited to signed 16 bit.
<lang tinybasic> LET P = 1
10 IF P > 2 THEN LET P = P + 2 IF P < 3 THEN LET P = P + 1 LET Z = P GOSUB 1000 IF A = 0 THEN GOTO 10 LET K = 0 20 LET K = K + 1 GOSUB 2000 LET Z = P - F IF Z < 0 THEN GOTO 30 GOSUB 1000 IF A = 1 THEN LET E = 0 IF A = 1 THEN GOTO 10 GOTO 20 30 LET C = C + 1 IF P < 2500 THEN PRINT C," ",P IF P > 2500 THEN END GOTO 10
1000 REM primality of Z by trial division, result is in A
LET Y = 1 LET A = 0 IF Z = 2 THEN LET A = 1 IF Z < 3 THEN RETURN
1010 LET Y = Y + 2
IF (Z/Y)*Y = Z THEN RETURN IF Y*Y < Z THEN GOTO 1010 LET A = 1 RETURN
2000 REM factorial of K, result is in F
LET A = 1 LET F = 1
2010 LET F = F*A
IF A=K THEN RETURN LET A = A + 1 GOTO 2010</lang>
- Output:
1 22 101 3 211 4 367 5 409 6 419 7 461 8 673 9 709 10 769 11 937 12 967 13 1009 14 1201 15 1259 16 1709 17 1831 18 2141 19 2221 20 2351 21 2411 22 2437
Wren
<lang ecmascript>import "/math" for Int import "/seq" for Lst import "/fmt" for Fmt
var limit = 1e6 var primes = Int.primeSieve(limit - 1, true) var erdos = [] for (p in primes) {
var i = 1 var fact = 1 var found = true while (fact < p) { if (Int.isPrime(p - fact)) { found = false break } i = i + 1 fact = fact * i } if (found) erdos.add(p)
} var lowerLimit = 2500 var erdosLower = erdos.where { |e| e < lowerLimit}.toList Fmt.print("The $,d Erdős primes under $,d are:", erdosLower.count, lowerLimit) for (chunk in Lst.chunks(erdosLower, 10)) Fmt.print("$6d", chunk) var show = 7875 Fmt.print("\nThe $,r Erdős prime is $,d.", show, erdos[show-1])</lang>
- Output:
The 25 Erdős primes under 2,500 are: 2 101 211 367 409 419 461 557 673 709 769 937 967 1009 1201 1259 1709 1831 1889 2141 2221 2309 2351 2411 2437 The 7,875th Erdős prime is 999,721.
XPL0
<lang XPL0>func IsPrime(N); \Return 'true' if N is prime int N, I; [if N <= 2 then return N = 2; if (N&1) = 0 then \even >2\ return false; for I:= 3 to sqrt(N) do
[if rem(N/I) = 0 then return false; I:= I+1; ];
return true; ];
func Erdos(N); \Return 'true' if N is an Erdos prime int N, F, K; [if not IsPrime(N) then return false; F:= 1; K:= 1; while F < N do
[if IsPrime(N-F) then return false; K:= K+1; F:= F*K; ];
return true; ];
int Cnt, N, SN; [Format(5, 0); Cnt:= 0; N:= 2; repeat if Erdos(N) then
[RlOut(0, float(N)); Cnt:= Cnt+1; if rem(Cnt/10) = 0 then CrLf(0); ]; N:= N+1;
until N >= 2500; CrLf(0); Text(0, "Found "); IntOut(0, Cnt); Text(0, " Erdos primes^m^j"); Cnt:= 1; N:= 3; repeat if Erdos(N) then [Cnt:= Cnt+1; SN:= N];
N:= N+2;
until N >= 1_000_000; Text(0, "The "); IntOut(0, Cnt); Text(0, "th Erdos prime is indeed "); IntOut(0, SN); CrLf(0); ]</lang>
- Output:
2 101 211 367 409 419 461 557 673 709 769 937 967 1009 1201 1259 1709 1831 1889 2141 2221 2309 2351 2411 2437 Found 25 Erdos primes The 7875th Erdos prime is indeed 999721