Yellowstone sequence
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).
- Task
- 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.
- Related tasks
- See also
-
- The OEIS entry: A098550 The Yellowstone permutation.
- Applegate et al, 2015: The Yellowstone Permutation [1].
Factor
<lang factor>USING: accessors assocs colors.constants combinators.short-circuit io kernel math prettyprint sequences sets 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:" print 30 <yellowstone> [ pprint bl ] each nl
! Plot first 100 Yellowstone numbers.
chart new { { 0 100 } { 0 175 } } >>axes line new COLOR: blue >>color 100 <iota> 100 <yellowstone> zip >>data add-gadget "Yellowstone numbers" open-window</lang>
- 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. <lang go>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()
}</lang>
- 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]
Haskell
<lang haskell>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</lang>
- 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]
A variation on the definition of the Yellowstone permutation (in terms of iterate, as an alternative to unfoldr), and a basic chart of the first 100 terms. (Plotting tested only on Haskell for Mac).
<lang haskell>import qualified Graphics.SVGFonts.ReadFont as F import Graphics.Rendering.Chart.Backend.Diagrams import Graphics.Rendering.Chart.Easy import Diagrams.Backend.Rasterific import Diagrams.Prelude import Codec.Picture import Control.Arrow (second)
YELLOWSTONE PERMUTATION------------------
yellowstone :: [Integer] yellowstone = 1 : 2 : (active <$> iterate nextWindow (2, 3, [4 ..]))
where
nextWindow :: (Integer, Integer, [Integer]) -> (Integer, Integer, [Integer])
nextWindow (p2, p1, rest) = (p1, n, residue)
where
[rp2, rp1] = relativelyPrime <$> [p2, p1]
go (x:xs)
| rp1 x && not (rp2 x) = (x, xs)
| otherwise = second ((:) x) (go xs)
(n, residue) = go rest
active (_, x, _) = x
relativelyPrime :: Integer -> Integer -> Bool relativelyPrime a b = 1 == gcd a b
30 FIRST TERMS, AND CHART OF FIRST 100-----------
main :: IO (Image PixelRGBA8) main = do
print $ take 30 yellowstone env <- chartEnv return $ chartRender env $ plot (line "Yellowstone terms" [zip [1 ..] (take 100 yellowstone)])
CHART GENERATION----------------------
chartRender
:: (Default r, ToRenderable r) => DEnv Double -> EC r () -> Image PixelRGBA8
chartRender env ec =
let (width, _) = envOutputSize env
in renderDia Rasterific (RasterificOptions (mkWidth width)) $
fst $ runBackendR env (toRenderable (execEC ec))
LOCAL FONT------------------------
chartEnv :: IO (DEnv Double) chartEnv = do
sansR <- F.loadFont "SourceSansPro_R.svg"
sansRB <- F.loadFont "SourceSansPro_RB.svg"
let fontChosen fs =
case (_font_name fs, _font_slant fs, _font_weight fs) of
("sans-serif", FontSlantNormal, FontWeightNormal) -> sansR
("sans-serif", FontSlantNormal, FontWeightBold) -> sansRB
return $ createEnv vectorAlignmentFns 640 400 fontChosen</lang>
- 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]
J
<lang J> GCD=: +. relatively_prime=: 1 = GCD
yellowstone=: monad define
start=. #\ i. 4 + y NB. prepare minimal starting values
s=. 3 {. start NB. the sequence vector
start=. 3 }. start
while. y > # s do.
z=. {. start NB. z is the lowest number not in the sequence
while.do.
if. 0 1 -: (_2 {. s) relatively_prime z do.
if. z -.@e. s do.
break.
end.
end.
z =. >: z
end.
start=. start -. z NB. remove z from the list of starting values
s=. s , z
end.
s
) </lang>
yellowstone 30 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 load'plot' 'marker'plot yellowstone 100
Java
<lang java> import java.util.ArrayList; import java.util.List;
public class YellowstoneSequence {
public static void main(String[] args) {
System.out.printf("First 30 values in the yellowstone sequence:%n%s%n", yellowstoneSequence(30));
}
private static List<Integer> yellowstoneSequence(int sequenceCount) {
List<Integer> yellowstoneList = new ArrayList<Integer>();
yellowstoneList.add(1);
yellowstoneList.add(2);
yellowstoneList.add(3);
int num = 4;
List<Integer> notYellowstoneList = new ArrayList<Integer>();
int yellowSize = 3;
while ( yellowSize < sequenceCount ) {
int found = -1;
for ( int index = 0 ; index < notYellowstoneList.size() ; index++ ) {
int test = notYellowstoneList.get(index);
if ( gcd(yellowstoneList.get(yellowSize-2), test) > 1 && gcd(yellowstoneList.get(yellowSize-1), test) == 1 ) {
found = index;
break;
}
}
if ( found >= 0 ) {
yellowstoneList.add(notYellowstoneList.remove(found));
yellowSize++;
}
else {
while ( true ) {
if ( gcd(yellowstoneList.get(yellowSize-2), num) > 1 && gcd(yellowstoneList.get(yellowSize-1), num) == 1 ) {
yellowstoneList.add(num);
yellowSize++;
num++;
break;
}
notYellowstoneList.add(num);
num++;
}
}
}
return yellowstoneList;
}
private static final int gcd(int a, int b) {
if ( b == 0 ) {
return a;
}
return gcd(b, a%b);
}
} </lang>
- Output:
First 30 values 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]
JavaScript
<lang javascript>(() => {
'use strict';
// yellowstone :: Generator [Int]
function* yellowstone() {
// A non finite stream of terms in the
// Yellowstone permutation of the natural numbers.
// OEIS A098550
const nextWindow = ([p2, p1, rest]) => {
const [rp2, rp1] = [p2, p1].map(
relativelyPrime
);
const go = xxs => {
const [x, xs] = Array.from(
uncons(xxs).Just
);
return rp1(x) && !rp2(x) ? (
Tuple(x)(xs)
) : secondArrow(cons(x))(
go(xs)
);
};
return [p1, ...Array.from(go(rest))];
};
const A098550 = fmapGen(x => x[1])(
iterate(nextWindow)(
[2, 3, enumFrom(4)]
)
);
yield 1
yield 2
while (true)(
yield A098550.next().value
)
};
// relativelyPrime :: Int -> Int -> Bool
const relativelyPrime = a =>
// True if a is relatively prime to b.
b => 1 === gcd(a)(b);
// ------------------------TEST------------------------
const main = () => console.log(
take(30)(
yellowstone()
)
);
// -----------------GENERIC FUNCTIONS------------------
// Just :: a -> Maybe a
const Just = x => ({
type: 'Maybe',
Nothing: false,
Just: x
});
// Nothing :: Maybe a
const Nothing = () => ({
type: 'Maybe',
Nothing: true,
});
// Tuple (,) :: a -> b -> (a, b)
const Tuple = a =>
b => ({
type: 'Tuple',
'0': a,
'1': b,
length: 2
});
// abs :: Num -> Num
const abs =
// Absolute value of a given number - without the sign.
Math.abs;
// cons :: a -> [a] -> [a]
const cons = x =>
xs => Array.isArray(xs) ? (
[x].concat(xs)
) : 'GeneratorFunction' !== xs
.constructor.constructor.name ? (
x + xs
) : ( // cons(x)(Generator)
function*() {
yield x;
let nxt = xs.next()
while (!nxt.done) {
yield nxt.value;
nxt = xs.next();
}
}
)();
// enumFrom :: Enum a => a -> [a]
function* enumFrom(x) {
// A non-finite succession of enumerable
// values, starting with the value x.
let v = x;
while (true) {
yield v;
v = 1 + v;
}
}
// fmapGen <$> :: (a -> b) -> Gen [a] -> Gen [b]
const fmapGen = f =>
function*(gen) {
let v = take(1)(gen);
while (0 < v.length) {
yield(f(v[0]))
v = take(1)(gen)
}
};
// gcd :: Int -> Int -> Int
const gcd = x => y => {
const
_gcd = (a, b) => (0 === b ? a : _gcd(b, a % b)),
abs = Math.abs;
return _gcd(abs(x), abs(y));
};
// iterate :: (a -> a) -> a -> Gen [a]
const iterate = f =>
function*(x) {
let v = x;
while (true) {
yield(v);
v = f(v);
}
};
// length :: [a] -> Int
const length = xs =>
// Returns Infinity over objects without finite
// length. This enables zip and zipWith to choose
// the shorter argument when one is non-finite,
// like cycle, repeat etc
(Array.isArray(xs) || 'string' === typeof xs) ? (
xs.length
) : Infinity;
// secondArrow :: (a -> b) -> ((c, a) -> (c, b))
const secondArrow = f => xy =>
// A function over a simple value lifted
// to a function over a tuple.
// f (a, b) -> (a, f(b))
Tuple(xy[0])(
f(xy[1])
);
// take :: Int -> [a] -> [a]
// take :: Int -> String -> String
const take = n =>
// The first n elements of a list,
// string of characters, or stream.
xs => 'GeneratorFunction' !== xs
.constructor.constructor.name ? (
xs.slice(0, n)
) : [].concat.apply([], Array.from({
length: n
}, () => {
const x = xs.next();
return x.done ? [] : [x.value];
}));
// uncons :: [a] -> Maybe (a, [a])
const uncons = xs => {
// Just a tuple of the head of xs and its tail,
// Or Nothing if xs is an empty list.
const lng = length(xs);
return (0 < lng) ? (
Infinity > lng ? (
Just(Tuple(xs[0])(xs.slice(1))) // Finite list
) : (() => {
const nxt = take(1)(xs);
return 0 < nxt.length ? (
Just(Tuple(nxt[0])(xs))
) : Nothing();
})() // Lazy generator
) : Nothing();
};
// MAIN --- return main();
})();</lang>
- 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
<lang 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 a
end
println("The first 30 entries of the Yellowstone permutation:\n", yellowstone(30))
x = 1:100 y = yellowstone(100) plot(x, y)
</lang>
- 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
<lang 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;</lang>
- 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
Phix
<lang Phix>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 a
end function
printf(1,"The first 30 entries of the Yellowstone permutation:\n%v\n", {yellowstone(30)})</lang>
- 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
<lang Phix>include pGUI.e IupOpen() 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()</lang>
PicoLisp
<lang PicoLisp>(load "@lib/frac.l") (de yellow (N)
(let (L (list 3 2 1) I 4 C 3 D)
(while (> N C)
(when
(and
(not (idx 'D I))
(=1 (gcd I (get L 1)))
(> (gcd I (get L 2)) 1) )
(push 'L I)
(idx 'D I T)
(setq I 4)
(inc 'C) )
(inc 'I) )
(flip L) ) )
(println (yellow 30))</lang>
- 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)
Python
<lang python>Yellowstone permutation OEIS A098550
from itertools import chain, count, islice from operator import itemgetter from 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 a
def 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 a
def 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 go
- 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[0], f(xy[1]))
- showList :: [a] -> String
def showList(xs):
Stringification of a list. return '[' + ','.join(repr(x) for x in xs) + ']'
- take :: Int -> [a] -> [a]
- take :: Int -> String -> String
def 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[0], xs[1:])) if xs else Nothing()
else:
nxt = take(1)(xs)
return Just((nxt[0], xs)) if nxt else Nothing()
- MAIN ---
if __name__ == '__main__':
main()</lang>
- 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]
Raku
(formerly Perl 6)
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.
<lang perl6>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 {+@x} 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;</lang>
- 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
REXX
horizontal list of numbers
<lang 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</lang>
- output when using the default input: 30
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
vertical histogram plot
A horizontal histogram could also be shown, but it would require a taller (higher) plot with more vertical screen real estate. <lang 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*/
call $histo $ '(vertical)' /*invoke a REXX vertical histogram plot*/ 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</lang>
- output when using the input: 100
The terminal screen size has a height of 55 characters (including values), the plot is shown at three quarter scale.
■
│
│
│ ■
■ │ │
│ │ ■ │
│ │ │ │
│ │ │ │
│ │ ■ │ │
│ │ │ │ │
│ │ │ │ │
│ │ ■ │ │ │
■ │ │ │ │ │ │
│ │ │ ■ │ │ │ │
│ │ │ │ │ ■ │ │ │
│ │ │ ■ │ │ │ │ │ │
■ │ │ │ ■ │ │ │ │ │ │ │
■ │ │ │ │ │ │ │ │ │ │ │ │
│ │ │ │ │ ■ │ │ │ │ │ │ │ │
│ │ │ │ │ │ │ │ │ │ │ │ │ │
│ │ │ │ │ ■ │ │ │ │ │ │ │ │ │
│ ■ │ │ │ │ │ │ │ │ │ │ │ │ │ │
│ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │
■ │ ■ │ │ │ │ │ │ │ │ │ │ │■│ │ │■│
■ │ │ │ │ │ │ │ │ │ │ │ │ │ │││ │ │││
■ │ │ │ ■ │ │ │ │ │ │ │ │ │ │■ │■│││ │ │││
│ │ │ │ │ │ │ │ │ │ │ │ │ │ ││ ■│││││ │■│││
│ │ │ │ ■ │ │ │ │ │ │ │ │ │ │■││ ││││││■│││││
│ │ │ │ ■ │ │ │ │ ■ │ │ │ │ │ │■││││ ││││││││││││
│ │ │ │ │ │ │ │ │ │ │ │ │■│■│ ││││││ ││││││││││││
│ │ │ ■ │ │ │ │ │ │■│ │ │ │││││ ■││││││ ││││││││││││
│ ■ │ │ │ │ │ │ │ │■│││ │ │■│││││ │││││││ ││││││││││││
│ │ │ ■ │ │ │ │ │ │■ │││││ │ │││││││ │││││││ ││││││││││││
│ │ │ │ │ ■ │ │ │ │ ││ ■│││││ │■│││││││ │││││││ ││││││││││││
│ ■ │ │ │ │ │ │ │ │ │■││ ││││││■│││││││││ │││││││ ││││││││││││
■ │ │ │ │ │ │ │ │■ │■│ ││││ ││││││││││││││││ │││││││ ││││││││││││
│ │ │ │ │ ■ │ │ │ ││ ■│││ ││││ ││││││││││││││││ │││││││ ││││││││││││
■ │ │ │ │ │ │ │ │■│ ││ ││││■││││ ││││││││││││││││ │││││││ ■││││││││││││
│ │ │ │ │ │ │ │ ■│││■││ │││││││││ ││││││││││││││││ │││││││■│││││││││││││
■ │ │ │ │ │■│ │■│ │││││││ │││││││││ ││││││││││││││││■│││││││││││││││││││││
■ │ ■ │ │ │ │■│││ │││ │││││││ │││││││││ ││││││││││││││││││││││││││││││││││││││
│ │ │ │ │■ │ │││││■│││ │││││││ │││││││││■■││││││││││││││││││││││││││││││││││││││
│ ■ │ │■│■││ │ │││││││││ │││││││ │││││││││││││││││││││││││││││││││││││││││││││││││
■ │ │■│ ││││││ │■│││││││││ │││││││■│││││││││││││││││││││││││││││││││││││││││││││││││
│ │ ■■│││ ││││││ │││││││││││■│││││││││││││││││││││││││││││││││││││││││││││││││││││││││
■ │ │■ │││││ ││││││■■│││││││││││││││││││││││││││││││││││││││││││││││││││││││││││││││││││││
│■ │■││ │││││ ■│││││││││││││││││││││││││││││││││││││││││││││││││││││││││││││││││││││││││││││
■■││ ││││ ■│││││■││││││││││││││││││││││││││││││││││││││││││││││││││││││││││││││││││││││││││││││
────┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴┴
1234981156213171222231132425311825363934361494546524157485862359696751616171713717181844818191718191
54 5256 0107291336581278545456109839254709038125616729183498741036627407016102613228322844546
5 9 3 1 5 5 1 7 3 9 5 5 5 5 5
Ruby
<lang ruby>def yellow(n)
a = [1, 2, 3]
b = { 1 => true, 2 => true, 3 => true }
i = 4
while n > a.length
if !b[i] && i.gcd(a[-1]) == 1 && i.gcd(a[-2]) > 1
a << i
b[i] = true
i = 4
end
i += 1
end
a
end
p yellow(30)</lang>
- 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]
Tcl
<lang Tcl>proc gcd {a b} {
while {$b} {
lassign [list $b [expr {$a % $b}]] a b
}
return $a
}
proc gen_yellowstones Template:MaxN 30 {
set r {}
for {set n 1} {$n <= $maxN} {incr n} {
if {$n <= 3} {
lappend r $n
} else {
## NB: list indices start at 0, not 1.
set pred [lindex $r end ] ;# a(n-1): coprime
set prepred [lindex $r end-1] ;# a(n-2): not coprime
for {set k 4} {1} {incr k} {
if {[lsearch -exact $r $k] >= 0} { continue }
if {1 != [gcd $k $pred ]} { continue }
if {1 == [gcd $k $prepred]} { continue }
## candidate k survived all tests...
break
}
lappend r $k
}
}
return $r
} puts "The first 30 Yellowstone numbers are:" puts [gen_yellowstones]</lang>
- 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
zkl
This sequence is limited to the max size of a Dictionary, 64k <lang zkl>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);
}</lang> <lang zkl>println("The first 30 entries of the Yellowstone permutation:"); yellowstoneW().walk(30).concat(", ").println();</lang>
- 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 <lang zkl>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();</lang> Offsite Image: yellowstone
