# A naturally occurring prime generator

There are quite a large number of prime generating functions. Pretty much all of them are constructed and are based on the sieve of Eratosthenes, the Möbius function or the Wilson's theorem and are generally infeasible to compute in practice. But there are also generators, that have a very easy structure and were found by accident.

In 2003 Stephen Wolfram explored a class of nested recurrence equations in a live computer experiment at the NKS Summer School. A group of people around Matthew Frank followed up with additional experiments and discovered an interesting property of the simply recurrence

a(n) = a(n-1) + gcd(n,a(n-1))


with the start value of a(1) = 7. The difference a(n) - a(n-1) = gcd(n,a(n-1)) always seemed to be 1 or a prime. The first few differences are (OEIS A132199):

1, 1, 1, 5, 3, 1, 1, 1, 1, 11, 3, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 23, 3, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 47, 3, 1, 5, 3, ...


If we only omit the 1s we get the following sequence (OEIS A137613):

5, 3, 11, 3, 23, 3, 47, 3, 5, 3, 101, 3, 7, 11, 3, 13, 233, 3, 467, 3, 5, 3,
941, 3, 7, 1889, 3, 3779, 3, 7559, 3, 13, 15131, 3, 53, 3, 7, 30323, 3, ...


Eric S. Rowland proved the primeness of each element in this list a few years later. As you can see, the primes are mixed and some of them appear multiple times. It also has been proven, that the sequence includes infinitely many distinct primes. Furthermore it is conjectured, that all odd primes appear.

Because this prime generator was not constructed but simply found by accident, the prime generator is called "naturally occurring". But notice that in practice this generator is also quite infeasible to compute. As it turns out, a prime p appears only after (p–3)/2 consecutive 1s. Nevertheless implementing this prime generator will be your task.

# Challenge:

Write a function or a program that prints the first n elements of the sequence A137613 (the sequence without the 1s). You can read the input number n >= 0 via STDIN, command-line argument, prompt or function argument. Output the first n elements in any readable format to STDOUT, or return an array or a list with these values.

This is code-golf. Therefore the shortest code wins.

# Language Name, N bytes


where N is the size of your submission. If you improve your score, you can keep old scores in the headline, by striking them through. For instance:

# Ruby, <s>104</s> <s>101</s> 96 bytes


var QUESTION_ID=55272;function answersUrl(e){return"http://api.stackexchange.com/2.2/questions/"+QUESTION_ID+"/answers?page="+e+"&pagesize=100&order=desc&sort=creation&site=codegolf&filter="+ANSWER_FILTER}function getAnswers(){jQuery.ajax({url:answersUrl(page++),method:"get",dataType:"jsonp",crossDomain:!0,success:function(e){answers.push.apply(answers,e.items),e.has_more?getAnswers():process()}})}function shouldHaveHeading(e){var a=!1,r=e.body_markdown.split("\n");try{a|=/^#/.test(e.body_markdown),a|=["-","="].indexOf(r[1][0])>-1,a&=LANGUAGE_REG.test(e.body_markdown)}catch(n){}return a}function shouldHaveScore(e){var a=!1;try{a|=SIZE_REG.test(e.body_markdown.split("\n")[0])}catch(r){}return a}function getAuthorName(e){return e.owner.display_name}function process(){answers=answers.filter(shouldHaveScore).filter(shouldHaveHeading),answers.sort(function(e,a){var r=+(e.body_markdown.split("\n")[0].match(SIZE_REG)||[1/0])[0],n=+(a.body_markdown.split("\n")[0].match(SIZE_REG)||[1/0])[0];return r-n});var e={},a=1,r=null,n=1;answers.forEach(function(s){var t=s.body_markdown.split("\n")[0],o=jQuery("#answer-template").html(),l=(t.match(NUMBER_REG)[0],(t.match(SIZE_REG)||[0])[0]),c=t.match(LANGUAGE_REG)[1],i=getAuthorName(s);l!=r&&(n=a),r=l,++a,o=o.replace("{{PLACE}}",n+".").replace("{{NAME}}",i).replace("{{LANGUAGE}}",c).replace("{{SIZE}}",l).replace("{{LINK}}",s.share_link),o=jQuery(o),jQuery("#answers").append(o),e[c]=e[c]||{lang:c,user:i,size:l,link:s.share_link}});var s=[];for(var t in e)e.hasOwnProperty(t)&&s.push(e[t]);s.sort(function(e,a){return e.lang>a.lang?1:e.lang<a.lang?-1:0});for(var o=0;o<s.length;++o){var l=jQuery("#language-template").html(),t=s[o];l=l.replace("{{LANGUAGE}}",t.lang).replace("{{NAME}}",t.user).replace("{{SIZE}}",t.size).replace("{{LINK}}",t.link),l=jQuery(l),jQuery("#languages").append(l)}}var ANSWER_FILTER="!t)IWYnsLAZle2tQ3KqrVveCRJfxcRLe",answers=[],page=1;getAnswers();var SIZE_REG=/\d+(?=[^\d&]*(?:&lt;(?:s&gt;[^&]*&lt;\/s&gt;|[^&]+&gt;)[^\d&]*)*)/,NUMBER_REG=/\d+/,LANGUAGE_REG=/^#*\s*([^,]+)/; body{text-align:left!important}#answer-list,#language-list{padding:10px;width:290px;float:left}table thead{font-weight:700}table td{padding:5px} <script src="https://ajax.googleapis.com/ajax/libs/jquery/2.1.1/jquery.min.js"></script><link rel="stylesheet" type="text/css" href="//cdn.sstatic.net/codegolf/all.css?v=83c949450c8b"><div id="answer-list"> <h2>Leaderboard</h2> <table class="answer-list"> <thead> <tr><td></td><td>Author</td><td>Language</td><td>Size</td></tr></thead> <tbody id="answers"> </tbody> </table></div><div id="language-list"> <h2>Winners by Language</h2> <table class="language-list"> <thead> <tr><td>Language</td><td>User</td><td>Score</td></tr></thead> <tbody id="languages"> </tbody> </table></div><table style="display: none"> <tbody id="answer-template"> <tr><td>{{PLACE}}</td><td>{{NAME}}</td><td>{{LANGUAGE}}</td><td>{{SIZE}}</td><td><a href="{{LINK}}">Link</a></td></tr></tbody></table><table style="display: none"> <tbody id="language-template"> <tr><td>{{LANGUAGE}}</td><td>{{NAME}}</td><td>{{SIZE}}</td><td><a href="{{LINK}}">Link</a></td></tr></tbody></table> • While the prime generator is not constructed, you're effectively implementing a trial division using the recursion. – orlp Aug 25, 2015 at 6:48 • If a(1) = 7, why doesn't the sequence begin with 7? Aug 25, 2015 at 6:54 • @feersum because the sequence we are concerned with is a(n)-a(n-1) Aug 25, 2015 at 6:55 • Can n be zero? Aug 25, 2015 at 7:50 • @jrenk Not sure. Maybe count it as 2 bytes (since you're removing 2 chars //) and explain it in your submission. If anyone disagrees with you, you can always edit your post. Aug 25, 2015 at 8:22 ## 16 Answers # Pyth, 14 13 bytes meaYhP+sY-5dQ  Uses a(n) = Lpf(6 - n + sum(a(i) for i in range(1, n)) where Lpf means least prime factor. Try it here online. # Python 3.5.0b1+, 95 93 bytes Link to Python 3.5.0b1+ release import math def f(k,n=2,a=7,L=[]):x=math.gcd(n,a);return k and f(k-1%x,n+1,a+x,L+1%x*[x])or L  A direct implementation of the recurrence, featuring: • What does 1%x do? Side question: where do I find documentation of the revision history of Python that includes betas? Edit: Nevermind, found it at the bottom of the revision history. Aug 27, 2015 at 4:33 • @mbomb007 Since x >= 1, 1%x returns 0 if x == 1, 1 otherwise (used to decide whether to add x to the list) Aug 27, 2015 at 4:42 # Julia, 110 bytes n->(a(n)=(n≥1&&(n==1?7:a(n-1)+gcd(n,a(n-1))));i=2;j=0;while j<n x=a(i)-a(i-1);x>1&&(j+=1;println(x));i+=1end)  Ungolfed: function a(n::Int) n ≥ 1 && (n == 1 ? 7 : a(n-1) + gcd(n, a(n-1))) end function f(n::Int) i = 2; j = 0; while j < n x = a(i) - a(i-1) if x > 1 j += 1 println(x) end i += 1 end end  • Wow, a perfect 8k, nice :D Aug 25, 2015 at 9:52 • Use n<2 instead of n==1. Also, if you look forwards instead of backwards, you can use i=1 and x=a(i)-a(i+=1), and then println(-x) and -x>1 to correct for the negativeness, thereby avoiding the need for a separate increment of i. And ≥ is three bytes, while >= is two... but then, you can use n<1||() rather than n>=1&&()... and yet, it's not even necessary in the first place (drop the conditional, n will never be less than 1). You also don't need the outermost brackets when defining a(n). With these changes, you should at least get down to 97 bytes. Aug 25, 2015 at 11:11 # PHP, 10196999877 72 bytes <?for(;2>t=gmp_strval(gmp_gcd(~++$i,7+$e+=$t))or$argv[1]-=print"$t ";);  Usage: Call the Script with an argument: php -d error_reporting=0 script.php 30 If you want to test it you need to uncomment ;extension=php_gmp.dll in your php.ini --> extension=php_gmp.dll Should I add the extension to my byte count? Any thoughts? Log: Saved 3 bytes thanks to Ismael Miguel. Saved 26 bytes thanks to primo. • You can shorten your opening tag to <? and remove the definition of $j. Aug 25, 2015 at 19:26
• Yes, it counts. But you can remove that newline. Which will save 1-2 bytes, depending on how you counted your code size. Aug 25, 2015 at 19:36
• Minor improvements: Use < in $j<=$argv[1] (prints one too many) (-1). Leave $e uninitialized, use $e+7 instead (-3). Use for(;;) instead of while(), making use of the pre- and post-expressions (-2). Replace echo$t.' ';$j++ with $j+=print"$t ", drop the brackets (-3). Replace if($t>1) with 2>$t|| (-2). Combine the assignment to $t with the conditional, switch || for or, drop the brackets (-5). Move $argv[1] to the $j increment, move the entire expression to the for conditional (-2). Change >=$j+=print to -=print (-3). Step by step: codepad.org/s6LNSPSM Aug 26, 2015 at 7:45
• @primo thanks for the nice explanation! Didn't know I could do all that. Aug 26, 2015 at 7:51
• A few more: Combine $e+7 with $e+=$t (-2). Leave $i uninitialized, use ~++$i instead (-3). codepad.org/fDIImajp Aug 26, 2015 at 8:03 # Haskell, 51 bytes d=zipWith gcd[2..]$scanl(+)7d
f=($filter(>1)d).take  Note that f is a function that will return the first n elements. Rather than computing the a(n) and then working out the differences, we compute the differences d(n) and sum them together to get a(n). (Those unfamiliar with Haskell may protest that we need a(n) first in order to get d(n), but of course lazy evaluation gets us around this problem!) Ungolfed: a = scanl (+) 7 d -- yielding a(n) = 7 + d(1) + d(2) + ... + d(n-1) d = zipWith gcd [2..] a -- yielding d(n) = gcd(n+1, a(n)) f n = take n$ filter (> 1) d -- get rid of 1s and take the first n


# Pyth, 30 bytes

Very badly golfed, can be considerably reduced. Defines recursive function at front, filters .first-n, and then maps the difference.

L?tb+KytbibK7m-yhdyd.ft-yhZyZQ

• This gives the wrong output for n = 0 Aug 25, 2015 at 8:16
• @Sp3000 that is a bug in Pyth. I'll put in a pull request. Aug 25, 2015 at 8:19
• Bug found and fixed - patch will be implemented once github stops being DDoS'd. Aug 25, 2015 at 11:49
• Here it is: meta.codegolf.stackexchange.com/questions/5318/…. Personally I'd consider bug fixes in programming languages as an answer Aug 25, 2015 at 21:28
• @ThomasWeller It kind of achieved the whole language ... Aug 26, 2015 at 2:42

# Julia, 69 67 bytes

n->(i=1;a=7;while n>0 x=gcd(i+=1,a);a+=x;x>1&&(n-=1;println(x))end)


This is a simple iterative solution to the problem. x is the difference (which is the gcd), and then I update a by adding x.

• I think it prints A231900. Aug 27, 2015 at 7:56
• @alephalpha - I think I see the error. Easily fixed. Even shaved two bytes off in the process. Aug 27, 2015 at 12:07

# JavaScript (ES6), 91

Recursive gcd, iterative main function. Not so fast.

Usual note: test running the snippet on any EcmaScript 6 compliant browser (notably not Chrome not MSIE. I tested on Firefox, Safari 9 could go)

F=m=>{
for(G=(a,b)=>b?G(b,a%b):a,o=[],p=7,n=1;m;d>1&&(o.push(d),--m))
p+=d=G(++n,p);
return o
}

O.innerHTML=F(+I.value)
<input id=I value=10><button onclick='O.innerHTML=F(+I.value)'>-></button>
<pre id=O></pre>

f=($filter(>1)$tail>>=zipWith(-)$scanl(\x->(x+).gcd x)7[2..]).take  Used the trick here: https://codegolf.stackexchange.com/a/39730/43318, and made point-free. (Previous: 71 bytes) a=scanl(\x->(x+).gcd x)7[2..] f m=take m$filter(>1)$zipWith(-)(tail a)a  First make the sequence of a's, and then take the differences. (Previous: 74 bytes) f m=take m$filter(>1)$map snd$scanl(\(x,d)->(\y->(x+y,y)).gcd x)(7,1)[2..]


Standard list functions, plus clever use of lambda function. Note this is 1 byte shorter than the more obvious

g m=take m$filter(>1)$map snd$scanl(\(x,d)n->(x+gcd x n,gcd x n))(7,1)[2..]  If we don't count imports, I can get this down to 66. import Data.List h m=take m$filter(>1)$snd$mapAccumL(\x->(\y->(x+y,y)).gcd x)7[2..]


## PARI/GP, 60 bytes

a(n)=a=7;i=1;while(n,if(1<d=gcd(i+=1,a),n-=1;print(d));a+=d)


Taken more or less straight from the definition a(n) - a(n-1) = gcd(n, a(n-1))

Output for a(20):

5
3
11
3
23
3
47
3
5
3
101
3
7
11
3
13
233
3
467
3


# APL (Dyalog Extended), 24 bytes

{⍵<2:5⋄⊃⍭6-⍵-+/∇¨⍳⍵-1}¨⍳


{⍵<2:5⋄⊃⍭6-⍵-+/∇¨⍳⍵-1} dfn to calculate the nth element

⍵<2:5 if ⍵<2 return 5

⋄ else

∇¨⍳⍵-1 recursively find the first ⍵-1 terms

+/ sum them

6-⍵-s in APL is 6-(⍵-s) = 6-⍵+s, as APL is evaluated right to left

⊃⍭ take the first (least) prime factor

¨⍳ apply the function to the first n numbers, to find the first n elements

Try it online!

{l←⍵⋄{l≡≢r←1~⍨2-⍨/⍵:r⋄∇⍵,(+/(0,1+≢)∨1↑⌽)⍵}7 8}


Try it online!

# Jelly, 13 10 bytes

’ßS_+6ÆfḢ)


Try it online!

Edit: I decided to revisit my oldest answer and found out that it can be very trivially golfed by 3 bytes. I sure have improved…

Uses the formula a(n) = Lpf(6-n+sum{i=1,...,n-1}a(i)).

## Explanation

’ßS_+6ÆfḢ)   Main monadic link
)   Map over the range [1..n]:
’              Decrement
S            Sum
_           Subtract n
Æf       Prime factorization
Ḣ      First element

• Welcome to the site, and impressive first answer! Oct 27, 2020 at 18:09
• @cairdcoinheringaahing Thanks! I've been lurking for a long time and decided to finally join. Oct 27, 2020 at 18:11
• If you're interested in improving your Jelly skills/talking to other Jelly users, be sure to check out the Jelly Hypertraining chat room as well! (You'll need 20 rep to do so, but I have no doubt that you can gather a couple of upvotes for your answer) Oct 27, 2020 at 18:12

# C++, 193182180 172 bytes

Thanks @Jakube - saved 8 bytes on output.

int g(int a,int b){return a==b?a:a>b?g(b,a-b):g(a,b-a);}void f(int *r,int n){int m=1,i=0,a=7,b;while(i<n){b=g(a,++m);if(b>1){r[i]=b;++i;}a+=b;}}int main(){int r[6];f(r,6);}

• You can probably save a few bytes by defining a function f, that returns an array with the results. This way you can remove the include, the scanf and the print. Aug 27, 2015 at 7:24

# Mathematica, 59 bytes

For[i=j=1;a=7,i<=#,,d=GCD[++j,a];If[d>1,Print@d;i++];a+=d]&


# Husk, 13 bytes

↑¡ö▼p+6§-LΣ;5


Try it online!

Uses Shevelev's formula from A137613: a(n) = Lpf(6-n+sum{i=1,...,n-1}a(i)).

↑               # input-th element of
¡              # list of repeated application of
ö             # 4 functions:
▼            # minimum of
p           # prime factors of
+6         # 6+
Σ     # sum of list so far
§-L      # minus length of list so far
;5   # starting with 5


# Husk, 18 bytes

↑f←Ẋ-¡λS+ȯ⌋→L¹→);7


Try it online!

Alternative derivation (following the story of the question): difference between adjacent elements of a(n) = a(n-1) + gcd(n,a(n-1)) with a(1)=7, ignoring 1s.

↑                   # input-th element of
f←                 # (without 1s)
Ẋ-               # differences between adjacent elements of
¡              # list of repeated application of
λ        )    # this function:
S+     →     # sum of last element of list so far, plus
ȯ⌋         # gcd of
→        # last element of list so far, and
L¹      # length of list so far
;7  # starting with 7


# Japt-g, 12 bytes

n6+UÆßXÃÅx)k


Try it

Uses the LPF formula from the OEIS entry with a little bit of Japt trickery: when reducing a 2D-array by addition Japt just grabs the first element of each sub-array*.

n6+UÆßXÃÅx)k     :Implicit input of integer U
n                :Subtract U from
6+              :  6 plus
UÆ            :  Map each X in the range [0,U)
ßX          :    Recursive call with argument U=X
Ã         :  End map
Å        :  Slice off the first element
x       :  Reduce by addition (of first elements)
)      :End subtraction
k     :Prime factors of result
:Implicit output of first element of final array


*Or, for those interested, what actually happens is that each sub-array is run through JavaScript's parseFloat function coercing them to comma delimited strings and then parseFloat parses that string to an integer until it hits an illegal character, i.e. the first comma.