# 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 '15 at 6:48 • If a(1) = 7, why doesn't the sequence begin with 7? – feersum Aug 25 '15 at 6:54 • @feersum because the sequence we are concerned with is a(n)-a(n-1) – Maltysen Aug 25 '15 at 6:55 • Can n be zero? – Sp3000 Aug 25 '15 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. – Jakube Aug 25 '15 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. – mbomb007 Aug 27 '15 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) – Sp3000 Aug 27 '15 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 – Beta Decay Aug 25 '15 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. – Glen O Aug 25 '15 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. – Ismael Miguel Aug 25 '15 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. – Ismael Miguel Aug 25 '15 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 – primo Aug 26 '15 at 7:45
• @primo thanks for the nice explanation! Didn't know I could do all that. – jrenk Aug 26 '15 at 7:51
• A few more: Combine $e+7 with $e+=$t (-2). Leave $i uninitialized, use ~++$i instead (-3). codepad.org/fDIImajp – primo Aug 26 '15 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 – Sp3000 Aug 25 '15 at 8:16
• @Sp3000 that is a bug in Pyth. I'll put in a pull request. – Maltysen Aug 25 '15 at 8:19
• Bug found and fixed - patch will be implemented once github stops being DDoS'd. – isaacg Aug 25 '15 at 11:49
• Here it is: meta.codegolf.stackexchange.com/questions/5318/…. Personally I'd consider bug fixes in programming languages as an answer – Thomas Weller Aug 25 '15 at 21:28
• @ThomasWeller It kind of achieved the whole language ... – isaacg Aug 26 '15 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. – alephalpha Aug 27 '15 at 7:56
• @alephalpha - I think I see the error. Easily fixed. Even shaved two bytes off in the process. – Glen O Aug 27 '15 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 (11) bytes

’ÑS_+6ÆfḢ
RÇ€


Try it online!

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

Some other answers return the nth term instead of the first n terms. If that is allowed, this can be 11 bytes: ’Rß€S_+6ÆfḢ

## Explanation

’ÑS_+6ÆfḢ   Monadic link that calculates a(n)
’           Decrement
S         Sum
_        Subtract [n]
Æf    Prime factorization
Ḣ   First element
R           Range from 1 [to n]
€         Map over the list...

• Welcome to the site, and impressive first answer! – caird coinheringaahing Oct 27 at 18:09
• @cairdcoinheringaahing Thanks! I've been lurking for a long time and decided to finally join. – xigoi Oct 27 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) – caird coinheringaahing Oct 27 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. – Jakube Aug 27 '15 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.