A naturally occurring prime generator

Question

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.

Leaderboard:

Here is a Stack Snippet to generate both a regular leaderboard and an overview of winners by language. To make sure that your answer shows up, please start your answer with a headline, using the following Markdown template:

# 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>

Answer 1

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.

Answer 2

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:

• Our good friend 1%x, and
• math.gcd, as opposed to fractions.gcd.

Answer 3

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

Answer 4

PHP, 101 96 99 98 77 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.

Answer 5

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

Answer 6

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

Try it here online.

Answer 7

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.

Answer 8

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>

Answer 9

Haskell, 74 71 66 bytes

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

Answer 10

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

Answer 11

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!

Old answer: 46 bytes

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

Try it online!

Answer 12

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
 ß             Recursively apply this link
  S            Sum
   _           Subtract n
    +6         Add 6
      Æf       Prime factorization
        Ḣ      First element

Answer 13

C++, 193 182 180 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);}

Answer 14

Mathematica, 59 bytes

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

Answer 15

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

Answer 16

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.