# Yellowstone sequence

Yellowstone sequence 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.

The Yellowstone sequence, also called the Yellowstone permutation, is defined as:

For n <= 3,

a(n) = n.

For n >= 4,

a(n) = the smallest number not already in sequence such that a(n) is relatively prime to a(n-1) and is not relatively prime to a(n-2).

The sequence is a permutation of the natural numbers, and gets its name from what its authors felt was a spiking, geyser like appearance of a plot of the sequence.

Example

a(4) is 4 because 4 is the smallest number following 1, 2, 3 in the sequence that is relatively prime to the entry before it (3), and is not relatively prime to the number two entries before it (2).

Find and show as output the first 30 Yellowstone numbers.

Extra
Demonstrate how to plot, with x = n and y coordinate a(n), the first 100 Yellowstone numbers.

## Factor

Works with: Factor version 0.99 2020-01-23
`USING: accessors assocs colors.constantscombinators.short-circuit io kernel math prettyprint sequencessets ui ui.gadgets ui.gadgets.charts ui.gadgets.charts.lines ; : yellowstone? ( n hs seq -- ? )    {        [ drop in? not ]        [ nip last gcd nip 1 = ]        [ nip dup length 2 - swap nth gcd nip 1 > ]    } 3&& ; : next-yellowstone ( hs seq -- n )    [ 4 ] 2dip [ 3dup yellowstone? ] [ [ 1 + ] 2dip ] until    2drop ; : next ( hs seq -- hs' seq' )    2dup next-yellowstone [ suffix! ] [ pick adjoin ] bi ; : <yellowstone> ( n -- seq )    [ HS{ 1 2 3 } clone dup V{ } set-like ] dip dup 3 <=    [ head nip ] [ 3 - [ next ] times nip ] if ;  ! Show first 30 Yellowstone numbers. "First 30 Yellowstone numbers:" print30 <yellowstone> [ pprint bl ] each nl ! Plot first 100 Yellowstone numbers. chart new { { 0 100 } { 0 175 } } >>axesline new COLOR: blue >>color100 <iota> 100 <yellowstone> zip >>dataadd-gadget "Yellowstone numbers" open-window`
Output:
```First 30 Yellowstone numbers:
1 2 3 4 9 8 15 14 5 6 25 12 35 16 7 10 21 20 27 22 39 11 13 33 26 45 28 51 32 17
```

## Go

This uses Gnuplot-X11 to do the plotting rather than a third party Go plotting library.

`package main import (    "fmt"    "log"    "os/exec") func gcd(x, y int) int {    for y != 0 {        x, y = y, x%y    }    return x} func yellowstone(n int) []int {    m := make(map[int]bool)    a := make([]int, n+1)    for i := 1; i < 4; i++ {        a[i] = i        m[i] = true    }    min := 4    for c := 4; c <= n; c++ {        for i := min; ; i++ {            if !m[i] && gcd(a[c-1], i) == 1 && gcd(a[c-2], i) > 1 {                a[c] = i                m[i] = true                if i == min {                    min++                }                break            }        }    }        return a[1:]} func check(err error) {    if err != nil {        log.Fatal(err)    }} func main() {    x := make([]int, 100)    for i := 0; i < 100; i++ {        x[i] = i + 1    }    y := yellowstone(100)    fmt.Println("The first 30 Yellowstone numbers are:")    fmt.Println(y[:30])    g := exec.Command("gnuplot", "-persist")    w, err := g.StdinPipe()    check(err)    check(g.Start())    fmt.Fprintln(w, "unset key; plot '-'")    for i, xi := range x {        fmt.Fprintf(w, "%d %d\n", xi, y[i])    }    fmt.Fprintln(w, "e")    w.Close()    g.Wait()}`
Output:
```The first 30 Yellowstone numbers are:
[1 2 3 4 9 8 15 14 5 6 25 12 35 16 7 10 21 20 27 22 39 11 13 33 26 45 28 51 32 17]
```

`import Data.List (unfoldr) yellowstone :: [Integer]yellowstone = 1 : 2 : 3 : unfoldr (Just . f) (2,3,[4..]) where  f :: (Integer, Integer, [Integer]) -> (Integer, (Integer, Integer, [Integer]))  f (p2, p1, rest) = (next, (p1, next, rest_)) where    (next, rest_) = select rest    select :: [Integer] -> (Integer, [Integer])    select (x:xs)      | gcd x p1 == 1 && gcd x p2 /= 1 = (x, xs)      | otherwise = (y, x:ys)      where (y, ys) = select xs main :: IO ()main = print \$ take 30 yellowstone`
Output:
```[1,2,3,4,9,8,15,14,5,6,25,12,35,16,7,10,21,20,27,22,39,11,13,33,26,45,28,51,32,17]
```

## Julia

`using Plots function yellowstone(N)    a = [1, 2, 3]    b = Dict(1 => 1, 2 => 1, 3 => 1)    start = 4    while length(a) < N        inseries = true        for i in start:typemax(Int)            if haskey(b, i)                if inseries                    start += 1                end            else                inseries = false            end            if !haskey(b, i) && (gcd(i, a[end]) == 1) && (gcd(i, a[end - 1]) > 1)                push!(a, i)                b[i] = 1                break            end        end    end    return aend println("The first 30 entries of the Yellowstone permutation:\n", yellowstone(30)) x = 1:100y = yellowstone(100)plot(x, y) `
Output:
```The first 30 entries of the Yellowstone permutation:
[1, 2, 3, 4, 9, 8, 15, 14, 5, 6, 25, 12, 35, 16, 7, 10, 21, 20, 27, 22, 39, 11, 13, 33, 26, 45, 28, 51, 32, 17]
```

## Perl

`use strict;use warnings;use feature 'say'; use List::Util qw(first);use GD::Graph::bars; use constant Inf  => 1e5; sub gcd {  my (\$u, \$v) = @_;  while (\$v) {    (\$u, \$v) = (\$v, \$u % \$v);  }  return abs(\$u);} sub yellowstone {    my(\$terms) = @_;    my @s = (1, 2, 3);    my @used = (1) x 4;    my \$min  = 3;    while (1) {        my \$index = first { not defined \$used[\$_] and gcd(\$_,\$s[-2]) != 1 and gcd(\$_,\$s[-1]) == 1 } \$min .. Inf;        \$used[\$index] = 1;        \$min = (first { not defined \$used[\$_] } 0..@used-1) || @used-1;        push @s, \$index;        last if @s == \$terms;    }    @s;} say "The first 30 terms in the Yellowstone sequence:\n" . join ' ', yellowstone(30); my @data = ( [1..500], [yellowstone(500)]);my \$graph = GD::Graph::bars->new(800, 600);\$graph->set(    title          => 'Yellowstone sequence',    y_max_value    => 1400,    x_tick_number  => 5,    r_margin       => 10,    dclrs          => [ 'blue' ],) or die \$graph->error;my \$gd = \$graph->plot(\@data) or die \$graph->error; open my \$fh, '>', 'yellowstone-sequence.png';binmode \$fh;print \$fh \$gd->png();close \$fh;`
Output:
```The first 30 terms in the Yellowstone sequence:
1 2 3 4 9 8 15 14 5 6 25 12 35 16 7 10 21 20 27 22 39 11 13 33 26 45 28 51 32 17```

See graph at off-site PNG image

## Perl 6

Works with: Rakudo version 2020.01

Not really clear whether a line graph or bar graph was desired, so generate both. Also, 100 points don't really give a good feel for the overall shape so do 500.

`my @yellowstone = 1, 2, 3, -> \$q, \$p {    state @used = True xx 4;    state \$min  = 3;    my \index = (\$min .. *).first: { not @used[\$_] and \$_ gcd \$q != 1 and \$_ gcd \$p == 1 };    @used[index] = True;    \$min = @used.first(!*, :k) // +@used - 1;    index} … *; put "The first 30 terms in the Yellowstone sequence:\n", @yellowstone[^30]; use SVG;use SVG::Plot; my @x = ^500; my \$chart = SVG::Plot.new(    background  => 'white',    width       => 1000,    height      => 600,    plot-width  => 950,    plot-height => 550,    x           => @x,    x-tick-step => { 10 },    y-tick-step => { 50 },    min-y-axis  => 0,    values      => [@yellowstone[@x],],    title       => "Yellowstone Sequence - First {[email protected]} values (zero indexed)",); my \$line = './Yellowstone-sequence-line-perl6.svg'.IO;my \$bars = './Yellowstone-sequence-bars-perl6.svg'.IO; \$line.spurt: SVG.serialize: \$chart.plot: :lines;\$bars.spurt: SVG.serialize: \$chart.plot: :bars;`
Output:
```The first 30 terms in the Yellowstone sequence:
1 2 3 4 9 8 15 14 5 6 25 12 35 16 7 10 21 20 27 22 39 11 13 33 26 45 28 51 32 17```

See (offsite SVG images) Line graph or Bar graph

## Phix

Translation of: Julia
`function yellowstone(integer N)    sequence a = {1, 2, 3},             b = repeat(true,3)    integer i = 4    while length(a) < N do        if (i>length(b) or b[i]=false)        and gcd(i,a[\$])=1        and gcd(i,a[\$-1])>1 then            a &= i            if i>length(b) then                b &= repeat(false,i-length(b))            end if            b[i] = true            i = 4        end if        i += 1    end while    return aend function printf(1,"The first 30 entries of the Yellowstone permutation:\n%v\n", {yellowstone(30)})`
Output:
```The first 30 entries of the Yellowstone permutation:
{1,2,3,4,9,8,15,14,5,6,25,12,35,16,7,10,21,20,27,22,39,11,13,33,26,45,28,51,32,17}
```

### a simple plot

Library: pGUI
`include pGUI.eIupOpen()IupControlsOpen()Ihandle plot = IupPlot("MENUITEMPROPERTIES=Yes, SIZE=640x320")IupSetAttribute(plot, "TITLE", "Yellowstone Numbers");IupSetAttribute(plot, "TITLEFONTSIZE", "10");IupSetAttribute(plot, "TITLEFONTSTYLE", "ITALIC");IupSetAttribute(plot, "GRIDLINESTYLE", "DOTTED");IupSetAttribute(plot, "GRID", "YES");IupSetAttribute(plot, "AXS_XLABEL", "n");IupSetAttribute(plot, "AXS_YLABEL", "a(n)");IupSetAttribute(plot, "AXS_XFONTSTYLE", "ITALIC");IupSetAttribute(plot, "AXS_YFONTSTYLE", "ITALIC");IupSetAttribute(plot, "AXS_YTICKSIZEAUTO", "NO");IupSetAttribute(plot, "AXS_YTICKMAJORSIZE", "8");IupSetAttribute(plot, "AXS_YTICKMINORSIZE", "0");IupPlotBegin(plot)sequence y500 = yellowstone(500)for x=1 to 500 do    IupPlotAdd(plot, x, y500[x])end for{} = IupPlotEnd(plot)--IupSetAttribute(plot, "DS_MODE", "BAR") -- (optional)Ihandle dlg = IupDialog(plot)IupCloseOnEscape(dlg)IupSetAttribute(dlg, "TITLE", "Yellowstone Names")IupMap(dlg)IupShowXY(dlg,IUP_CENTER,IUP_CENTER)IupMainLoop()IupClose()`

## Python

Works with: Python version 3.7
`'''Yellowstone permutation OEIS A098550''' from itertools import chain, count, islicefrom operator import itemgetterfrom math import gcd from matplotlib import pyplot  # yellowstone :: [Int]def yellowstone():    '''A non-finite stream of terms from       the Yellowstone permutation.       OEIS A098550.    '''    # relativelyPrime :: Int -> Int -> Bool    def relativelyPrime(a):        return lambda b: 1 == gcd(a, b)     # nextWindow :: (Int, Int, [Int]) -> (Int, Int, [Int])    def nextWindow(triple):        p2, p1, rest = triple        [rp2, rp1] = map(relativelyPrime, [p2, p1])         # match :: [Int] -> (Int, [Int])        def match(xxs):            x, xs = uncons(xxs)['Just']            return (x, xs) if rp1(x) and not rp2(x) else (                second(cons(x))(                    match(xs)                )            )        n, residue = match(rest)        return (p1, n, residue)     return chain(        range(1, 3),        map(            itemgetter(1),            iterate(nextWindow)(                (2, 3, count(4))            )        )    )  # TEST ----------------------------------------------------# main :: IO ()def main():    '''Terms of the Yellowstone permutation.'''     print(showList(        take(30)(yellowstone())    ))    pyplot.plot(        take(100)(yellowstone())    )    pyplot.xlabel(main.__doc__)    pyplot.show()  # GENERIC ------------------------------------------------- # Just :: a -> Maybe adef Just(x):    '''Constructor for an inhabited Maybe (option type) value.       Wrapper containing the result of a computation.    '''    return {'type': 'Maybe', 'Nothing': False, 'Just': x}  # Nothing :: Maybe adef Nothing():    '''Constructor for an empty Maybe (option type) value.       Empty wrapper returned where a computation is not possible.    '''    return {'type': 'Maybe', 'Nothing': True}  # cons :: a -> [a] -> [a]def cons(x):    '''Construction of a list from x as head,       and xs as tail.    '''    return lambda xs: [x] + xs if (        isinstance(xs, list)    ) else x + xs if (        isinstance(xs, str)    ) else chain([x], xs)  # iterate :: (a -> a) -> a -> Gen [a]def iterate(f):    '''An infinite list of repeated       applications of f to x.    '''    def go(x):        v = x        while True:            yield v            v = f(v)    return lambda x: go(x)  # second :: (a -> b) -> ((c, a) -> (c, b))def second(f):    '''A simple function lifted to a function over a tuple,       with f applied only to the second of two values.    '''    return lambda xy: (xy, f(xy))  # showList :: [a] -> Stringdef showList(xs):    '''Stringification of a list.'''    return '[' + ','.join(repr(x) for x in xs) + ']'  # take :: Int -> [a] -> [a]# take :: Int -> String -> Stringdef take(n):    '''The prefix of xs of length n,       or xs itself if n > length xs.    '''    return lambda xs: (        xs[0:n]        if isinstance(xs, (list, tuple))        else list(islice(xs, n))    )  # uncons :: [a] -> Maybe (a, [a])def uncons(xs):    '''The deconstruction of a non-empty list       (or generator stream) into two parts:       a head value, and the remaining values.    '''    if isinstance(xs, list):        return Just((xs, xs[1:])) if xs else Nothing()    else:        nxt = take(1)(xs)        return Just((nxt, xs)) if nxt else Nothing()  # MAIN ---if __name__ == '__main__':    main()`
Output:
`1,2,3,4,9,8,15,14,5,6,25,12,35,16,7,10,21,20,27,22,39,11,13,33,26,45,28,51,32,17]`

## REXX

`/*REXX program calculates any number of terms in the Yellowstone (permutation) sequence.*/parse arg m .                                    /*obtain optional argument from the CL.*/if m=='' | m==","  then m= 30                    /*Not specified?  Then use the default.*/!.= 0                                            /*initialize an array of numbers(used).*/# = 0                                            /*count of Yellowstone numbers in seq. */\$=                                               /*list   "      "         "     "  "   */      do j=1  until #==m;  prev= # - 1      if j<5  then do;  #= #+1;   @.#= j;  !.#= j;  !.j= 1;  \$= strip(\$ j);  iterate;  end          do k=1;   if !.k  then iterate          /*Already used?  Then skip this number.*/         if gcd(k, @.#)\==1  |  gcd(k, @.prev)<2  then iterate   /*not meet requirement?*/         #= #+1;   @.#= k;     !.k= 1;   \$= \$ k  /*bump ctr; assign; mark used; add list*/         leave                                   /*find the next Yellowstone seq. number*/         end   /*k*/      end      /*j*/say \$                                            /*display a list of a Yellowstone seq. */exit                                             /*stick a fork in it,  we're all done. *//*──────────────────────────────────────────────────────────────────────────────────────*/gcd: parse arg x,y;  do until y==0;  parse value  x//y  y   with   y  x;  end;    return x`
output   when using the default input:
```1 2 3 4 9 8 15 14 5 6 25 12 35 16 7 10 21 20 27 22 39 11 13 33 26 45 28 51 32 17
```

## zkl

Translation of: Julia

This sequence is limited to the max size of a Dictionary, 64k

`fcn yellowstoneW{	// --> iterator   Walker.zero().tweak(fcn(a,b){      foreach i in ([1..]){         if(not b.holds(i) and i.gcd(a[-1])==1 and i.gcd(a[-2]) >1){	    a.del(0).append(i);	// only keep last two terms	    b[i]=True;	    return(i);	 }      }   }.fp(List(2,3), Dictionary(1,True, 2,True, 3,True))).push(1,2,3);}`
`println("The first 30 entries of the Yellowstone permutation:");yellowstoneW().walk(30).concat(", ").println();`
Output:
```The first 30 entries of the Yellowstone permutation:
1, 2, 3, 4, 9, 8, 15, 14, 5, 6, 25, 12, 35, 16, 7, 10, 21, 20, 27, 22, 39, 11, 13, 33, 26, 45, 28, 51, 32, 17
```

Plot using Gnuplot

`gnuplot:=System.popen("gnuplot","w");gnuplot.writeln("unset key; plot '-'"); yellowstoneW().pump(1_000, gnuplot.writeln.fp(" "));  // " 1\n", " 2\n", ...gnuplot.writeln("e");gnuplot.flush();ask("Hit return to finish"); gnuplot.close();`

Offsite Image: yellowstone