Yellowstone sequence
The Yellowstone sequence, also called the Yellowstone permutation, is defined as:
You are encouraged to solve this task according to the task description, using any language you may know.
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].
T YellowstoneGenerator
min_ = 1
n_ = 0
n1_ = 0
n2_ = 0
Set[Int] sequence_
F next()
.n2_ = .n1_
.n1_ = .n_
I .n_ < 3
.n_++
E
.n_ = .min_
L !(.n_ !C .sequence_ & gcd(.n1_, .n_) == 1 & gcd(.n2_, .n_) > 1)
.n_++
.sequence_.add(.n_)
L
I .min_ !C .sequence_
L.break
.sequence_.remove(.min_)
.min_++
R .n_
print(‘First 30 Yellowstone numbers:’)
V ygen = YellowstoneGenerator()
print(ygen.next(), end' ‘’)
L(i) 1 .< 30
print(‘ ’ygen.next(), end' ‘’)
print()- 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
F yellow(n)
V a = [1, 2, 3]
V b = Set([1, 2, 3])
V i = 4
L n > a.len
I i !C b & gcd(i, a.last) == 1 & gcd(i, a[(len)-2]) > 1
a.append(i)
b.add(i)
i = 4
i++
R a
print(yellow(30))- 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]
BYTE FUNC Gcd(BYTE a,b)
BYTE tmp
IF a<b THEN
tmp=a a=b b=tmp
FI
WHILE b#0
DO
tmp=a MOD b
a=b b=tmp
OD
RETURN (a)
BYTE FUNC Contains(BYTE ARRAY a BYTE len,value)
BYTE i
FOR i=0 TO len-1
DO
IF a(i)=value THEN
RETURN (1)
FI
OD
RETURN (0)
PROC Generate(BYTE ARRAY seq BYTE count)
BYTE i,x
seq(0)=1 seq(1)=2 seq(2)=3
FOR i=3 TO COUNT-1
DO
x=1
DO
IF Contains(seq,i,x)=0 AND
Gcd(x,seq(i-1))=1 AND Gcd(x,seq(i-2))>1 THEN
EXIT
FI
x==+1
OD
seq(i)=x
OD
RETURN
PROC Main()
DEFINE COUNT="30"
BYTE ARRAY seq(COUNT)
BYTE i
Generate(seq,COUNT)
PrintF("First %B Yellowstone numbers:%E",COUNT)
FOR i=0 TO COUNT-1
DO
PrintB(seq(i)) Put(32)
OD
RETURN- Output:
Screenshot from Atari 8-bit computer
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
with Ada.Text_IO;
with Ada.Containers.Ordered_Sets;
procedure Yellowstone_Sequence is
generic -- Allow more than one generator, but must be instantiated
package Yellowstones is
function Next return Integer;
function GCD (Left, Right : Integer) return Integer;
end Yellowstones;
package body Yellowstones
is
package Sequences is
new Ada.Containers.Ordered_Sets (Integer);
-- Internal package state
N_0 : Integer := 0;
N_1 : Integer := 0;
N_2 : Integer := 0;
Seq : Sequences.Set;
Min : Integer := 1;
function GCD (Left, Right : Integer) return Integer
is (if Right = 0
then Left
else GCD (Right, Left mod Right));
function Next return Integer is
begin
N_2 := N_1;
N_1 := N_0;
if N_0 < 3 then
N_0 := N_0 + 1;
else
N_0 := Min;
while
not (not Seq.Contains (N_0)
and then GCD (N_1, N_0) = 1
and then GCD (N_2, N_0) > 1)
loop
N_0 := N_0 + 1;
end loop;
end if;
Seq.Insert (N_0);
while Seq.Contains (Min) loop
Seq.Delete (Min);
Min := Min + 1;
end loop;
return N_0;
end Next;
end Yellowstones;
procedure First_30 is
package Yellowstone is new Yellowstones; -- New generator instance
use Ada.Text_IO;
begin
Put_Line ("First 30 Yellowstone numbers:");
for A in 1 .. 30 loop
Put (Yellowstone.Next'Image); Put (" ");
end loop;
New_Line;
end First_30;
begin
First_30;
end Yellowstone_Sequence;
- 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
scope # Yellowstone sequence - translation of the Pluto sample
local proc yellowstone( n :: number ) :: table
local a, m := [ 1, 2, 3 ], [ true, true, true ];
for x from 4 to n do
a[ x ], m[ x ] := 0, false
od;
local minV := 4;
for c from 4 to n do
local more := true
for i from minV while more do
if not m[ i ] and numtheory.gcd( a[ c - 1 ], i ) = 1 and numtheory.gcd( a[ c - 2 ], i ) > 1
then
a[ c ], m[ i ] := i, true;
if i = minV then minV +:= 1 fi;
more := false
fi
od
od;
return a
end;
local constant ySize, constant perLine := 30, 10;
local y := yellowstone( ySize );
printf( "The first %d Yellowstone numbers are:\n", ySize );
for yPos to ySize do printf( " %2d", y[ yPos ] ); if yPos mod perLine = 0 then print() fi od
end- 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
BEGIN # find members of the yellowstone sequence: starting from 1, 2, 3 the #
# subsequent members are the lowest number coprime to the previous one #
# and not coprime to the one before that, that haven't appeared in the #
# sequence yet #
# iterative Greatest Common Divisor routine, returns the gcd of m and n #
PROC gcd = ( INT m, n )INT:
BEGIN
INT a := ABS m, b := ABS n;
WHILE b /= 0 DO
INT new a = b;
b := a MOD b;
a := new a
OD;
a
END # gcd # ;
# returns an array of the Yellowstone seuence up to n #
OP YELLOWSTONE = ( INT n )[]INT:
BEGIN
[ 1 : n ]INT result;
IF n > 0 THEN
result[ 1 ] := 1;
IF n > 1 THEN
result[ 2 ] := 2;
IF n > 2 THEN
result[ 3 ] := 3;
# guess the maximum element will be n, if it is larger, used will be enlarged #
HEAP[ 1 : n ]BOOL initial used;
REF[]BOOL used := initial used;
used[ 1 ] := used[ 2 ] := used[ 3 ] := TRUE;
FOR i FROM 4 TO UPB used DO used[ i ] := FALSE OD;
FOR i FROM 4 TO UPB result DO
INT p1 = result[ i - 1 ];
INT p2 = result[ i - 2 ];
BOOL found := FALSE;
FOR j WHILE NOT found DO
IF j > UPB used THEN
# not enough elements in used - enlarge it #
HEAP[ 1 : 2 * UPB used ]BOOL new used;
new used[ 1 : UPB used ] := used;
FOR k FROM UPB used + 1 TO UPB new used DO new used[ k ] := FALSE OD;
used := new used
FI;
IF NOT used[ j ] THEN
IF found := gcd( j, p1 ) = 1 AND gcd( j, p2 ) /= 1
THEN
result[ i ] := j;
used[ j ] := TRUE
FI
FI
OD
OD
FI
FI
FI;
result
END # YELLOWSTONE # ;
[]INT ys = YELLOWSTONE 30;
FOR i TO UPB ys DO
print( ( " ", whole( ys[ i ], 0 ) ) )
OD;
print( ( newline ) )
END- 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
yellowstone: function [n][
result: [1 2 3]
present: [1 2 3]
start: 4
while [n > size result][
candidate: new start
while ø [
if all? @[
not? contains? present candidate
1 = gcd @[candidate last result]
1 <> gcd @[candidate get result (size result)-2]
][
'result ++ candidate
'present ++ candidate
while [contains? present start] -> inc 'start
break
]
inc 'candidate
]
]
return result
]
print yellowstone 30
- 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 := [], in_seq := []
loop 30 {
n := A_Index
if n <=3
A[n] := n, in_seq[n] := true
else while true
{
s := A_Index
if !in_seq[s] && relatively_prime(s, A[n-1]) && !relatively_prime(s, A[n-2])
{
A[n] := s
in_seq[s] := true
break
}
}
}
for i, v in A
result .= v ","
MsgBox % result := "[" Trim(result, ",") "]"
return
;--------------------------------------
relatively_prime(a, b){
return (GCD(a, b) = 1)
}
;--------------------------------------
GCD(a, b) {
while b
b := Mod(a | 0x0, a := b)
return a
}
- 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]
#include <stdbool.h>
#include <stdio.h>
#include <stdlib.h>
typedef struct lnode_t {
struct lnode_t *prev;
struct lnode_t *next;
int v;
} Lnode;
Lnode *make_list_node(int v) {
Lnode *node = malloc(sizeof(Lnode));
if (node == NULL) {
return NULL;
}
node->v = v;
node->prev = NULL;
node->next = NULL;
return node;
}
void free_lnode(Lnode *node) {
if (node == NULL) {
return;
}
node->v = 0;
node->prev = NULL;
free_lnode(node->next);
node->next = NULL;
}
typedef struct list_t {
Lnode *front;
Lnode *back;
size_t len;
} List;
List *make_list() {
List *list = malloc(sizeof(List));
if (list == NULL) {
return NULL;
}
list->front = NULL;
list->back = NULL;
list->len = 0;
return list;
}
void free_list(List *list) {
if (list == NULL) {
return;
}
list->len = 0;
list->back = NULL;
free_lnode(list->front);
list->front = NULL;
}
void list_insert(List *list, int v) {
Lnode *node;
if (list == NULL) {
return;
}
node = make_list_node(v);
if (list->front == NULL) {
list->front = node;
list->back = node;
list->len = 1;
} else {
node->prev = list->back;
list->back->next = node;
list->back = node;
list->len++;
}
}
void list_print(List *list) {
Lnode *it;
if (list == NULL) {
return;
}
for (it = list->front; it != NULL; it = it->next) {
printf("%d ", it->v);
}
}
int list_get(List *list, int idx) {
Lnode *it = NULL;
if (list != NULL && list->front != NULL) {
int i;
if (idx < 0) {
it = list->back;
i = -1;
while (it != NULL && i > idx) {
it = it->prev;
i--;
}
} else {
it = list->front;
i = 0;
while (it != NULL && i < idx) {
it = it->next;
i++;
}
}
}
if (it == NULL) {
return INT_MIN;
}
return it->v;
}
///////////////////////////////////////
typedef struct mnode_t {
int k;
bool v;
struct mnode_t *next;
} Mnode;
Mnode *make_map_node(int k, bool v) {
Mnode *node = malloc(sizeof(Mnode));
if (node == NULL) {
return node;
}
node->k = k;
node->v = v;
node->next = NULL;
return node;
}
void free_mnode(Mnode *node) {
if (node == NULL) {
return;
}
node->k = 0;
node->v = false;
free_mnode(node->next);
node->next = NULL;
}
typedef struct map_t {
Mnode *front;
} Map;
Map *make_map() {
Map *map = malloc(sizeof(Map));
if (map == NULL) {
return NULL;
}
map->front = NULL;
return map;
}
void free_map(Map *map) {
if (map == NULL) {
return;
}
free_mnode(map->front);
map->front = NULL;
}
void map_insert(Map *map, int k, bool v) {
if (map == NULL) {
return;
}
if (map->front == NULL) {
map->front = make_map_node(k, v);
} else {
Mnode *it = map->front;
while (it->next != NULL) {
it = it->next;
}
it->next = make_map_node(k, v);
}
}
bool map_get(Map *map, int k) {
if (map != NULL) {
Mnode *it = map->front;
while (it != NULL && it->k != k) {
it = it->next;
}
if (it != NULL) {
return it->v;
}
}
return false;
}
///////////////////////////////////////
int gcd(int u, int v) {
if (u < 0) u = -u;
if (v < 0) v = -v;
if (v) {
while ((u %= v) && (v %= u));
}
return u + v;
}
List *yellow(size_t n) {
List *a;
Map *b;
int i;
a = make_list();
list_insert(a, 1);
list_insert(a, 2);
list_insert(a, 3);
b = make_map();
map_insert(b, 1, true);
map_insert(b, 2, true);
map_insert(b, 3, true);
i = 4;
while (n > a->len) {
if (!map_get(b, i) && gcd(i, list_get(a, -1)) == 1 && gcd(i, list_get(a, -2)) > 1) {
list_insert(a, i);
map_insert(b, i, true);
i = 4;
}
i++;
}
free_map(b);
return a;
}
int main() {
List *a = yellow(30);
list_print(a);
free_list(a);
putc('\n', stdout);
return 0;
}
- 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
#include <iostream>
#include <numeric>
#include <set>
template <typename integer>
class yellowstone_generator {
public:
integer next() {
n2_ = n1_;
n1_ = n_;
if (n_ < 3) {
++n_;
} else {
for (n_ = min_; !(sequence_.count(n_) == 0
&& std::gcd(n1_, n_) == 1
&& std::gcd(n2_, n_) > 1); ++n_) {}
}
sequence_.insert(n_);
for (;;) {
auto it = sequence_.find(min_);
if (it == sequence_.end())
break;
sequence_.erase(it);
++min_;
}
return n_;
}
private:
std::set<integer> sequence_;
integer min_ = 1;
integer n_ = 0;
integer n1_ = 0;
integer n2_ = 0;
};
int main() {
std::cout << "First 30 Yellowstone numbers:\n";
yellowstone_generator<unsigned int> ygen;
std::cout << ygen.next();
for (int i = 1; i < 30; ++i)
std::cout << ' ' << ygen.next();
std::cout << '\n';
return 0;
}
- 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
internal class Program
{
private static void Main(string[] args)
{
List<int> yellowStoneList = YellowstoneSequence(30);
Console.WriteLine("Yellowstone 30");
foreach (int i in yellowStoneList)
{
Console.Write(i + " ");
}
}
private static List<int> YellowstoneSequence(int sequenceCount)
{
List<int> yellowstoneList = [1, 2, 3];
int num = 4;
List<int> notYellowstoneList = new();
int yellowSize = 3;
while (yellowSize < sequenceCount)
{
int found = -1;
for (int index = 0; index < notYellowstoneList.Count; index++)
{
int test = notYellowstoneList[index];
if (GCD(yellowstoneList[yellowSize - 2], test) > 1 && GCD(yellowstoneList[yellowSize - 1], test) == 1)
{
found = index;
break;
}
}
if (found >= 0)
{
yellowstoneList.Add(notYellowstoneList[found]);
notYellowstoneList.RemoveAt(found);
yellowSize++;
}
else
{
while (true)
{
if (GCD(yellowstoneList[yellowSize - 2], num) > 1 && GCD(yellowstoneList[yellowSize - 1], num) == 1)
{
yellowstoneList.Add(num);
yellowSize++;
num++;
break;
}
notYellowstoneList.Add(num);
num++;
}
}
}
return yellowstoneList;
}
private static int GCD(int a, int b)
{
if (b == 0)
{
return a;
}
return GCD(b, a % b);
}
}
- Output:
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
The bar graph looks awful with a lot of monospace fonts (like Consolas or Courier New) since they lack some glyphs which are substituted from another font with different metrics. If it doesn't look good to you, I recommend using fonts with a better Unicode coverage, like Iosevka.
def yellowstone
used = Set(Int32).new
penult, last = 2, 3
start = 4
iter = Iterator.of do
while start.in? used
used.delete start
start += 1
end
(start..).each do |n|
if !used.includes?(n) && n.gcd(last) == 1 && n.gcd(penult) > 1
used << n
penult, last = last, n
break n
end
end
end
(1..3).each.chain iter
end
p yellowstone.first(30).to_a
values = yellowstone.first(100).to_a
height = (values.max - 1) // 8
(0..height).reverse_each do |i|
values.each do |v|
d = v - i*8
print d <= 0 ? ' ' : d >= 8 ? '█' : '\u2580' + d
end
puts
end
- 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]
▇
▁ █ ▅
█ █ ▃ █
█ █ ▁ █ █
█ █ █ █ █
▅ █ █ ▁ ▇ █ █ █
█ █ █ ▃ █ █ ▅ █ █ █
▃ ▇ █ █ █ ▅ █ █ █ █ █ █ █
█ █ █ █ █ ▁ ▇ █ █ █ █ █ █ █ █
█ ▃ █ █ █ █ █ █ █ █ █ █ █ █ █ █
▃ ▇ █ ▅ █ █ █ █ █ █ █ █ █ █ █▄█ █ █▆█
▅ █ █ █ ▁ ▇ █ █ █ █ █ █ █ █ █ █▆ ▂█████ █▄███
█ █ █ █ ▃ █ █ █ █ ▅ █ █ █ █▂█ █▄████ ██████▆█████
█ ▁ █ █ ▅ █ █ █ █ █ █▄█ █ █▂█████ ▆██████ ████████████
█ █ █ ▇ █ ▁ █ █ █ █ █▆ ▂█████ █▄███████ ███████ ████████████
▃ █ ▇ █ █ ▁ █ █ █ █▂ █▄█ ████ ██████▆█████████ ███████ ████████████
▅ █ █ █ █ █ █ █ ▂█▄█ ██ ████▆████ ████████████████ ███████▁▃████████████
▃ ▇ ▁ █ █ █ █▂█▄█ █▆█ ███████ █████████ ████████████████▅█████████████████████
▁ █ ▃ █ █▂█▄██ █ █████▆███ ███████ █████████▅▇██████████████████████████████████████
█ █ ▅▄█▆█ ██████▁▂███████████▃███████▇█████████████████████████████████████████████████
▁ ▇▆ █▄██ ▂█████▃▅█████████████████████████████████████████████████████████████████████████████
▁▂▃▄████▅▆████▇█████████████████████████████████████████████████████████████████████████████████████
import std.numeric;
import std.range;
import std.stdio;
class Yellowstone {
private bool[int] sequence_;
private int min_ = 1;
private int n_ = 0;
private int n1_ = 0;
private int n2_ = 0;
public this() {
popFront();
}
public bool empty() {
return false;
}
public int front() {
return n_;
}
public void popFront() {
n2_ = n1_;
n1_ = n_;
if (n_ < 3) {
++n_;
} else {
for (n_ = min_;
!(n_ !in sequence_ && gcd(n1_, n_) == 1 && gcd(n2_, n_) > 1);
++n_) {
// empty
}
}
sequence_[n_] = true;
while (true) {
if (min_ !in sequence_) {
break;
}
sequence_.remove(min_);
++min_;
}
}
}
void main() {
new Yellowstone().take(30).writeln();
}
- 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]
Boost.Generics.Collection and Boost.Process are part of DelphiBoostLib.
program Yellowstone_sequence;
{$APPTYPE CONSOLE}
uses
System.SysUtils,
Boost.Generics.Collection,
Boost.Process;
function gdc(x, y: Integer): Integer;
begin
while y <> 0 do
begin
var tmp := x;
x := y;
y := tmp mod y;
end;
Result := x;
end;
function Yellowstone(n: Integer): TArray<Integer>;
var
m: TDictionary<Integer, Boolean>;
a: TArray<Integer>;
begin
m.Init;
SetLength(a, n + 1);
for var i := 1 to 3 do
begin
a[i] := i;
m[i] := True;
end;
var min := 4;
for var c := 4 to n do
begin
var i := min;
repeat
if not m[i, false] and (gdc(a[c - 1], i) = 1) and (gdc(a[c - 2], i) > 1) then
begin
a[c] := i;
m[i] := true;
if i = min then
inc(min);
Break;
end;
inc(i);
until false;
end;
Result := copy(a, 1, length(a));
end;
begin
var x: TArray<Integer>;
SetLength(x, 100);
for var i in Range(100) do
x[i] := i + 1;
var y := yellowstone(High(x));
writeln('The first 30 Yellowstone numbers are:');
for var i := 0 to 29 do
Write(y[i], ' ');
Writeln;
//Plotting
var plot := TPipe.Create('gnuplot -p', True);
plot.WritelnA('unset key; plot ''-''');
for var i := 0 to High(x) do
plot.WriteA('%d %d'#10, [x[i], y[i]]);
plot.WritelnA('e');
writeln('Press enter to close');
Readln;
plot.Kill;
plot.Free;
end.
- 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 Press enter to close
func gcd a b .
if b = 0 : return a
return gcd b (a mod b)
.
proc remove_at i &a[] .
for j = i + 1 to len a[] : a[j - 1] = a[j]
len a[] -1
.
proc yellowstone count &yellow[] .
yellow[] = [ 1 2 3 ]
num = 4
while len yellow[] < count
yell1 = yellow[len yellow[] - 1]
yell2 = yellow[len yellow[]]
for i to len notyellow[]
test = notyellow[i]
if gcd yell1 test > 1 and gcd yell2 test = 1 : break 1
.
if i <= len notyellow[]
yellow[] &= notyellow[i]
remove_at i notyellow[]
else
while gcd yell1 num <= 1 or gcd yell2 num <> 1
notyellow[] &= num
num += 1
.
yellow[] &= num
num += 1
.
.
.
print "First 30 values in the yellowstone sequence:"
yellowstone 30 yellow[]
print yellow[]- 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 ]
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
- 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
(do ;;; Yellowstone sequence - translation of the Pluto sample
(fn gcd [m n]
(var (a b) (values (math.abs m) (math.abs n)))
(while (not= b 0)
(set (b a) (values (% a b) b))
)
a
)
(fn yellowstone [n]
(var (a m) (values [1 2 3] [true true true]))
(for [x 4 n]
(tset a x 0)
(tset m x false)
)
(var minV 4)
(for [c 4 n]
(var (more i) (values true minV))
(while more
(when (and (not (. m i))
(= (gcd (. a (- c 1)) i) 1)
(> (gcd (. a (- c 2)) i) 1)
)
(tset a c i)
(tset m i true)
(when (= i minV)
(set minV (+ 1 minV))
)
(set more false)
)
(set i (+ i 1))
)
)
a
)
(local (ySize perLine) (values 30 10))
(local y (yellowstone ySize))
(io.write (string.format "The first %d Yellowstone numbers are:\n" ySize))
(for [yPos 1 ySize]
(io.write (string.format " %2d" (. y yPos)))
(when (= 0 (% yPos perLine)) (print))
)
)
- 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
: array create cells allot ;
: th cells + ; \ some helper words
30 constant #yellow \ number of yellowstones
#yellow array y \ create array
( n1 n2 -- n3)
: gcd dup if tuck mod recurse exit then drop ;
: init 3 0 do i 1+ y i th ! loop ; ( --)
: show cr #yellow 0 do y i th ? loop ; ( --)
: gcd-y[] - cells y + @ over gcd ; ( k i n -- k gcd )
: loop1 begin 1+ over 2 gcd-y[] 1 = >r over 1 gcd-y[] 1 > r> or 0= until ;
: loop2 over true swap 0 ?do over y i th @ = if 0= leave then loop ;
: yellow #yellow 3 do i 3 begin loop1 loop2 until y rot th ! loop ;
: main init yellow show ;
main
- Output:
main 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 ok
! Yellowstone Sequence - Code Golf Version
! Demonstrates concise Fortran (but keep it readable!)
program yellowstone_golf
implicit none
integer :: y(100), i, c
logical :: u(1000) = .false.
! First 3 terms
y(1:3) = [(i, i=1,3)]
u(1:3) = .true.
! Generate rest
do i = 4, 100
c = 1
do while (u(c) .or. gcd(c,y(i-1))/=1 .or. gcd(c,y(i-2))==1)
c = c + 1
end do
y(i) = c
u(c) = .true.
end do
! Output
print '(A)', "First 30 Yellowstone numbers:"
print '(10I5)', y(1:30)
contains
! pure integer function gcd(a,b) !Euclid's method
! integer, intent(in) :: a, b
! integer :: x, y, t
! x = a; y = b
! do while (y /= 0)
! t = y; y = mod(x,y); x = t
! end do
! gcd = x
! end function
pure integer function gcd(a, b) !Stein's method
! Added and used Steins method of finding the GCD
! simply because everyone else was using Euclid's method
integer, intent(in) :: a, b
integer :: x, y, shift
x = abs(a)
y = abs(b)
! Handle trivial cases
if (a == 0) then
gcd = y
return
else if (b == 0) then
gcd = x
return
end if
! Count common factors of 2
shift = 0
do while ((iand(x,1) == 0) .and. (iand(y,1) == 0))
x = ishft(x,-1)
y = ishft(y,-1)
shift = shift + 1
end do
! Divide x by 2 until odd
do while (iand(x,1) == 0)
x = ishft(x,-1)
end do
do ! Divide y by 2 until odd
do while (iand(y,1) == 0)
y = ishft(y,-1)
end do
! Ensure x <= y
if (x > y) then
! swap
x = x + y
y = x - y
x = x - y
end if
y = y - x
if (y == 0) exit
end do
gcd = ishft(x, shift)
end function gcd
end program yellowstone_golf
- 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
function gcd(a as uinteger, b as uinteger) as uinteger
if b = 0 then return a
return gcd( b, a mod b )
end function
dim as uinteger i, j, k, Y(1 to 100)
Y(1) = 1 : Y(2) = 2: Y(3) = 3
for i = 4 to 100
k = 3
print i
do
k += 1
if gcd( k, Y(i-2) ) = 1 orelse gcd( k, Y(i-1) ) > 1 then continue do
for j = 1 to i-1
if Y(j)=k then continue do
next j
Y(i) = k
exit do
loop
next i
for i = 1 to 30
print str(Y(i))+" ";
next i
print
screen 13
for i = 1 to 100
pset (i, 200-Y(i)), 31
next i
while inkey=""
wend
end
- 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
FB 7.0.34 macOS 14.7.2 Sonoma
include "NSLog.incl"
NSInteger local fn GCD( a as NSInteger, b as NSInteger )
if ( a == 0 ) then return b
return fn GCD( b % a, a )
end fn = 0
void local fn CalcYellowstoneSequence( terms as int )
CFMutableArrayRef sequence = fn MutableArrayWithObjects( @1, @2, @3, NULL )
for int n = 3 to terms - 1
NSInteger prev1 = fn NumberIntegerValue( fn ArrayObjectAtIndex( sequence, n - 1 ) )
NSInteger prev2 = fn NumberIntegerValue( fn ArrayObjectAtIndex( sequence, n - 2 ) )
NSInteger candidate = 1
while (YES)
if ( !fn ArrayContainsObject( sequence, @(candidate) ) && fn GCD( candidate, prev1 ) == 1 && fn GCD( candidate, prev2 ) > 1 )
MutableArrayAddObject( sequence, @(candidate) )
break
end if
candidate++
wend
next
NSLog( @"Yellowstone sequence up to %d terms:", terms )
for CFNumberRef num in sequence
NSLog( @"%@ \b", num )
next
end fn
fn CalcYellowstoneSequence( 30 )
HandleEvents- Output:
Yellowstone sequence up to 30 terms: 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
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]
Or, defining the Yellowstone permutation in terms of iterate, rather than unfoldr,
and displaying a chart of the first 100 terms:
import Codec.Picture
import Data.Bifunctor (second)
import Diagrams.Backend.Rasterific
import Diagrams.Prelude
import Graphics.Rendering.Chart.Backend.Diagrams
import Graphics.Rendering.Chart.Easy
import qualified Graphics.SVGFonts.ReadFont as F
----------------- YELLOWSTONE PERMUTATION ----------------
yellowstone :: [Integer]
yellowstone =
1 :
2 :
(active <$> iterate nextWindow (2, 3, [4 ..]))
where
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 =
renderDia
Rasterific
( RasterificOptions
(mkWidth (fst (envOutputSize env)))
)
$ 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
- 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]
// @test: yellowstone(30L)
gcd x = match x with
| (a, 0) -> a
| (a, b) -> gcd((b, a % b))
coprime x = match x with
| (a, b) -> gcd((a, b)) == 1
usedContains x = match x with
| (val, seq, j) -> if (j >= length(seq)) then false else if (seq[j] == val) then true else usedContains((val, seq, j + 1))
findNext x = match x with
| (seq, prev1, prev2, candidate) -> if (usedContains((candidate, seq, 0L))) then findNext((seq, prev1, prev2, candidate + 1)) else if (coprime((candidate, prev1)) and not(coprime((candidate, prev2)))) then candidate else findNext((seq, prev1, prev2, candidate + 1))
yellowstoneHelper x = match x with
| (seq, i, n) -> if (i >= n) then seq else let next = findNext((seq, seq[i-1], seq[i-2], 1)) in yellowstoneHelper((seq ++ [next], i + 1, n))
yellowstone n = yellowstoneHelper(([1, 2, 3], 3L, n))
- 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]
tacit
Until=: 2 :'u^:(0-:v)^:_'
assert 44 -: >:Until(>&43) 32 NB. increment until exceeding 43
gcd=: +.
coprime=: 1 = gcd
prepare=:1 2 3"_ NB. start with the vector 1 2 3
condition=: 0 1 -: (coprime _2&{.) NB. trial coprime most recent 2, nay and yay
append=: , NB. concatenate
novel=: -.@e. NB. x is not a member of y
term=: >:@:]Until((condition *. novel)~) 4:
ys=: (append term)@]^:(0 >. _3+[) prepare
assert (ys 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
explicit
GCD=: +.
relatively_prime=: 1 = GCD
yellowstone=: {{
s=. 1 2 3 NB. initial sequence
while. y > # s do.
z=. <./(1+s)-.s NB. lowest positive inteeger not in sequence
while. if. 0 1 -: z relatively_prime _2{.s do. z e. s end. do.
z=. z+1
end. NB. find next value for sequence
s=. s, z
end.
}}
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
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);
}
}
- 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]
(() => {
'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();
})();
- 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
# jq optimizes the recursive call of _gcd in the following:
def gcd(a;b):
def _gcd:
if .[1] != 0 then [.[1], .[0] % .[1]] | _gcd else .[0] end;
[a,b] | _gcd ;
# emit the yellowstone sequence as a stream
def yellowstone:
1,2,3,
({ a: [2, 3], # the last two items only
b: {"1": true, "2": true, "3" : true}, # a record, to avoid having to save the entire history
start: 4 }
| foreach range(1; infinite) as $n (.;
first(
.b as $b
| .start = first( range(.start;infinite) | select($b[tostring]|not) )
| foreach range(.start; infinite) as $i (.;
.emit = null
| ($i|tostring) as $is
| if .b[$is] then .
# "a(n) is relatively prime to a(n-1) and is not relatively prime to a(n-2)"
elif (gcd($i; .a[1]) == 1) and (gcd($i; .a[0]) > 1)
then .emit = $i
| .a = [.a[1], $i]
| .b[$is] = true
else .
end;
select(.emit)) );
.emit ));The task
"The first 30 entries of the Yellowstone permutation:",
[limit(30;yellowstone)]- 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]
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)
- 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]
GcdStep = b~, a mod b, a mod b != 0
Gcd = GcdStep.while(a, b):1
FindNext(history, pre1, pre2) = {
IsYSCandidate = not history.contains(candidate) and
Gcd(candidate, pre1) == 1 and Gcd(candidate, pre2) != 1
FindStep = candidate + 1, not IsYSCandidate
FindStep.while(1):0
}
YSStep(history, pre2, pre1) = {
Next = FindNext(history, pre1, pre2)
(history..., Next), pre1, Next
}
YSStep.repeat(27, (1, 2, 3), 2, 3):0fun main() {
println("First 30 values in the yellowstone sequence:")
println(yellowstoneSequence(30))
}
private fun yellowstoneSequence(sequenceCount: Int): List<Int> {
val yellowstoneList = mutableListOf(1, 2, 3)
var num = 4
val notYellowstoneList = mutableListOf<Int>()
var yellowSize = 3
while (yellowSize < sequenceCount) {
var found = -1
for (index in notYellowstoneList.indices) {
val test = notYellowstoneList[index]
if (gcd(yellowstoneList[yellowSize - 2], test) > 1 && gcd(
yellowstoneList[yellowSize - 1], test
) == 1
) {
found = index
break
}
}
if (found >= 0) {
yellowstoneList.add(notYellowstoneList.removeAt(found))
yellowSize++
} else {
while (true) {
if (gcd(yellowstoneList[yellowSize - 2], num) > 1 && gcd(
yellowstoneList[yellowSize - 1], num
) == 1
) {
yellowstoneList.add(num)
yellowSize++
num++
break
}
notYellowstoneList.add(num)
num++
}
}
}
return yellowstoneList
}
private fun gcd(a: Int, b: Int): Int {
return if (b == 0) {
a
} else gcd(b, a % b)
}
- 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]
function gcd(a, b)
if b == 0 then
return a
end
return gcd(b, a % b)
end
function printArray(a)
io.write('[')
for i,v in pairs(a) do
if i > 1 then
io.write(', ')
end
io.write(v)
end
io.write(']')
return nil
end
function removeAt(a, i)
local na = {}
for j,v in pairs(a) do
if j ~= i then
table.insert(na, v)
end
end
return na
end
function yellowstone(sequenceCount)
local yellow = {1, 2, 3}
local num = 4
local notYellow = {}
local yellowSize = 3
while yellowSize < sequenceCount do
local found = -1
for i,test in pairs(notYellow) do
if gcd(yellow[yellowSize - 1], test) > 1 and gcd(yellow[yellowSize - 0], test) == 1 then
found = i
break
end
end
if found >= 0 then
table.insert(yellow, notYellow[found])
notYellow = removeAt(notYellow, found)
yellowSize = yellowSize + 1
else
while true do
if gcd(yellow[yellowSize - 1], num) > 1 and gcd(yellow[yellowSize - 0], num) == 1 then
table.insert(yellow, num)
yellowSize = yellowSize + 1
num = num + 1
break
end
table.insert(notYellow, num)
num = num + 1
end
end
end
return yellow
end
function main()
print("First 30 values in the yellowstone sequence:")
printArray(yellowstone(30))
print()
end
main()
- 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]
import "int"
yellowstone = function(n)
m = {}
a = [0] * (n + 1)
for i in range(1, 3)
a[i] = i
m[i] = true
end for
min = 4
for c in range(4, n, 1)
i = min
while true
if not m.hasIndex(i) and int.gcd(a[c-1], i) == 1 and int.gcd(a[c-2], i) > 1 then
a[c] = i
m[i] = true
if i == min then min += 1
break
end if
i += 1
end while
end for
return a[1:]
end function
print "The first 30 Yellowstone numbers are:"
print yellowstone(30)
- 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]
state = {1, 2, 3};
MakeNext[state_List] := Module[{i = First[state], done = False, out},
While[! done,
If[FreeQ[state, i],
If[GCD[Last[state], i] == 1,
If[GCD[state[[-2]], i] > 1,
out = Append[state, i];
done = True;
]
]
];
i++;
];
out
]
Nest[MakeNext, state, 30 - 3]
ListPlot[%]
- 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}
(* Graphical visualisation of the data *)
Procedure version
This version uses a set and, so, is limited to 65536 elements. It is easy to change this limit by using a HashSet (standard module “sets”) instead of a set. See the iterator version which uses such a HashSet.
import math
proc yellowstone(n: int): seq[int] =
assert n >= 3
result = @[1, 2, 3]
var present = {1, 2, 3}
var start = 4
while result.len < n:
var candidate = start
while true:
if candidate notin present and gcd(candidate, result[^1]) == 1 and gcd(candidate, result[^2]) != 1:
result.add candidate
present.incl candidate
while start in present: inc start
break
inc candidate
echo yellowstone(30)
- 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]
Iterator version
This version uses a HashSet, but using a set as in the previous version is possible if we accept the limit of 65536 elements.
import math, sets
iterator yellowstone(n: int): int =
assert n >= 3
for i in 1..3: yield i
var present = [1, 2, 3].toHashSet
var prevLast = 2
var last = 3
var start = 4
for _ in 4..n:
var candidate = start
while true:
if candidate notin present and gcd(candidate, last) == 1 and gcd(candidate, prevLast) != 1:
yield candidate
present.incl candidate
prevLast = last
last = candidate
while start in present: inc start
break
inc candidate
for n in yellowstone(30):
stdout.write " ", n
echo()
- 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
As Oberon-07 doesn't have dynamic arrays, either fixed size "large enough" arrays can be used or alternatively, dynamically allocated structures. This sample uses "large enough" arrays. see the Ulam numbers task for an alternative method.
MODULE YellowstoneSequence; (* Show some of the Yellowstone sequence *)
(* - translation of the Pluto sample via Agena *)
IMPORT Out;
(* returns the gcd of a and b *)
PROCEDURE gcd( a, b : INTEGER ) : INTEGER;
VAR t, m, n : INTEGER;
BEGIN
m := a;
n := b;
WHILE m # 0 DO t := m; m := n MOD m; n := t END
RETURN n
END gcd ;
(* sets a to the first elements of the Yellowstone sequence; *)
(* a must have at least 4 elements *)
(* element 0 of a is not used; i.e., a is indexed from 1 *)
(* the required size of the work array m will depend on the size of a *)
(* experimentation suggests that for sizes of a up to 100 000, *)
(* 10 * the size of a should be OK - it would probably be OK for *)
(* 1 000 000 but it will take a long time to construct a 1 000 000 *)
(* element sequence so I haven't tried it... *)
PROCEDURE yellowstone( VAR a : ARRAY OF INTEGER; VAR m : ARRAY OF BOOLEAN );
VAR n, minV, c, i : INTEGER;
more : BOOLEAN;
BEGIN
n := LEN( a ) - 1;
FOR i := 1 TO 3 DO a[ i ] := i; m[ i ] := TRUE END;
FOR i := 4 TO n DO a[ i ] := 0 END;
FOR i := 4 TO LEN( m ) - 1 DO m[ i ] := FALSE END;
minV := 4;
FOR c := 4 TO n DO
more := TRUE;
i := minV - 1;
WHILE more DO
INC( i );
IF ( ~ m[ i ] ) & ( gcd( a[ c - 1 ], i ) = 1 ) & ( gcd( a[ c - 2 ], i ) > 1 ) THEN
a[ c ] := i; m[ i ] := TRUE;
IF i = minV THEN INC( minV ) END;
more := FALSE
END
END
END
END yellowstone;
(* show the first 30 elements of the sequence *)
PROCEDURE task;
CONST ySize = 30;
VAR y : ARRAY ySize + 2 OF INTEGER;
w : ARRAY ySize * 10 OF BOOLEAN;
yPos : INTEGER;
BEGIN
yellowstone( y, w );
Out.String( "The first ");Out.Int( ySize, 0 );Out.String( " Yellowstone numbers are:" );Out.Ln;
FOR yPos := 1 TO ySize DO
Out.String( " " );Out.Int( y[ yPos ], 0 )
END;
Out.Ln
END task;
BEGIN
task
END YellowstoneSequence.
- 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
yellowstone(n) = {
my(a=3, o=2, u=[]);
if(n<3, return(n)); \\ Base case: return n if it is less than 3
print1("1, 2"); \\ Print initial values
for(i = 4, n, \\ Iterate from 4 to n
print1(", "a); \\ Print current value of a
u = setunion(u, Set(a)); \\ Add a to the set u
\\ Remove consecutive elements from u
while(#u > 1 && u[2] == u[1] + 1,
u = vecextract(u, "^1")
);
\\ Find next value of a
for(k = u[1] + 1, 1e10,
if(gcd(k, o) <= 1, next); \\ Skip if gcd(k, o) is greater than 1
if(setsearch(u, k), next); \\ Skip if k is in set u
if(gcd(k, a) != 1, next); \\ Skip if gcd(k, a) is not 1
o = a; \\ Update o to current a
a = k; \\ Update a to k
break
)
);
a \\ Return the final value of a
}
yellowstone(20); \\ Call the function with n = 20- Output:
1, 2, 3, 4, 9, 8, 15, 14, 5, 6, 25, 12, 35, 16, 7, 10, 21, 20, 27
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
You can run this online here.
-- -- demo\rosetta\Yellowstone_sequence.exw -- with javascript_semantics requires("1.0.2") 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)}) -- a simple plot: include pGUI.e include IupGraph.e function get_data(Ihandle graph) sequence y500 = yellowstone(500) integer {w,h} = IupGetIntInt(graph,"DRAWSIZE") IupSetInt(graph,"XTICK",iff(w<640?iff(h<300?100:50):20)) IupSetInt(graph,"YTICK",iff(h<250?iff(h<140?iff(h<120?700:350):200):100)) return {{tagset(500),y500,CD_RED}} end function IupOpen() Ihandle graph = IupGraph(get_data,"RASTERSIZE=960x600") IupSetAttributes(graph,`GTITLE="Yellowstone Numbers"`) IupSetInt(graph,"TITLESTYLE",CD_ITALIC) IupSetAttributes(graph,`XNAME="n", YNAME="a(n)"`) IupSetAttributes(graph,"XTICK=20,XMIN=0,XMAX=500") IupSetAttributes(graph,"YTICK=100,YMIN=0,YMAX=1400") Ihandle dlg = IupDialog(graph,`TITLE="Yellowstone Names"`) IupSetAttributes(dlg,"MINSIZE=290x140") IupShow(dlg) if platform()!=JS then IupMainLoop() IupClose() end if
- 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}
Require Utilitys library version 1.3
include ..\Utilitys.pmt
def gcd /# u v -- n #/
abs int swap abs int swap
dup
while
over over mod rot drop dup
endwhile
drop
enddef
def test enddef
def yellow var n
( 1 2 3 ) var a
newd ( 1 true ) setd ( 2 true ) setd ( 3 true ) setd var b
4 var i
test
while
b i getd "Unfound" == >ps
a -1 get >ps -2 get
i gcd 1 > ps> i gcd 1 == ps>
and and if
i 0 put var a
( i true ) setd var b
4 var i
else
drop drop
endif
i 1 + var i
test
endwhile
a
enddef
def test n a len nip > enddef
"The first 30 entries of the Yellowstone permutation:" ? 30 yellow ?- 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] === Press any key to exit ===
(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))- 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)
require "table2"
local int = require "int"
local fmt = require "fmt"
local function yellowstone(n)
local m = {}
local a = table.rep(n, 0)
for i = 1, 3 do
a[i] = i
m[i] = true
end
local min = 4
for c = 4, n do
local i = min
while true do
if !m[i] and int.gcd(a[c - 1], i) == 1 and int.gcd(a[c - 2], i) > 1 then
a[c] = i
m[i] = true
if i == min then ++min end
break
end
++i
end
end
return a
end
local x = {}
for i = 1, 30 do x[i] = i end
local y = yellowstone(30)
print("The first 30 Yellowstone numbers are:")
fmt.tprint("%2d ", y, 10)
- 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
Procedure.i gcd(x.i,y.i)
While y<>0 : t=x : x=y : y=t%y : Wend : ProcedureReturn x
EndProcedure
If OpenConsole()
Dim Y.i(100)
For i=1 To 100
If i<=3 : Y(i)=i : Continue : EndIf : k=3
Repeat
RepLoop:
k+1
For j=1 To i-1 : If Y(j)=k : Goto RepLoop : EndIf : Next
If gcd(k,Y(i-2))=1 Or gcd(k,Y(i-1))>1 : Continue : EndIf
Y(i)=k : Break
ForEver
Next
For i=1 To 30 : Print(Str(Y(i))+" ") : Next : Input()
EndIf
- 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
'''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()
- 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]
gcd is defined at Greatest common divisor#Quackery.
[ stack ] is seqbits ( --> s )
[ bit
seqbits take |
seqbits put ] is seqadd ( n --> )
[ bit
seqbits share & not ] is notinseq ( n --> b )
[ temp put
' [ 1 2 3 ]
7 seqbits put
4
[ dip
[ dup -1 peek
over -2 peek ]
dup dip
[ tuck gcd 1 !=
unrot gcd 1 =
and ]
swap if
[ dup dip join
seqadd
3 ]
[ 1+
dup notinseq until ]
over size temp share
< not until ]
drop
seqbits release
temp take split drop ] is yellowstones ( n --> [ )
30 yellowstones echo- 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 ]
#lang racket
(require plot)
(define a098550
(let ((hsh# (make-hash '((1 . 1) (2 . 2) (3 . 3))))
(rev# (make-hash '((1 . 1) (2 . 2) (3 . 3)))))
(λ (n)
(hash-ref hsh# n
(λ ()
(let ((a_n (for/first ((i (in-naturals 4))
#:unless (hash-has-key? rev# i)
#:when (and (= (gcd i (a098550 (- n 1))) 1)
(> (gcd i (a098550 (- n 2))) 1)))
i)))
(hash-set! hsh# n a_n)
(hash-set! rev# a_n n)
a_n))))))
(map a098550 (range 1 (add1 30)))
(plot (points
(map (λ (i) (vector i (a098550 i))) (range 1 (add1 100)))))
- Output:
Just the output text... you'll have to run this yourself in racket to see the plot!
'(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)
(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.
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;
- 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
Include: How to use
Include: Source code
-- 13 Jun 2026
include Setting
numeric digits 100
arg xx
if xx = '' then
xx = 30
say 'YELLOWSTONE SEQUENCE'
say version
say
call GetYellow xx
call DisplayList xx
call Timer
exit
GetYellow:
procedure expose yell. work.
arg xx
say 'Get yellowstones...'
yell. = 0; work. = 0; n = 0
do i = 1 until n = xx
p = n-1
if i < 5 then do
n = n+1; yell.n = i; work.n = i; work.i = 1
iterate i
end
do j = 1
if work.j then
iterate j
if Gcd(j,yell.p) = 1 then
iterate j
if Gcd(j,yell.n) <> 1 then
iterate j
n = n+1; yell.n = j; work.j = 1
leave j
end
end
say
return xx
DisplayList:
procedure expose yell.
arg xx
say 'Yellowstone sequence...'
do i = 1 to xx
call Charout ,Right(yell.i,5)
if i//10 = 0 then
say
end
say
return
include Math
- Output:
YELLOWSTONE SEQUENCE
REXX-Regina_3.9.7(MT) 5.00 18 Mar 2025
Get yellowstones...
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
0.001 seconds
see "working..." + nl
row = 3
num = 2
numbers = 1:51
first = 2
second = 3
see "Yellowstone numbers are:" + nl
see "1 " + first + " " + second + " "
for n = 4 to len(numbers)
flag1 = 1
flag2 = 1
if first < numbers[n]
min = first
else
min = numbers[n]
ok
for m = 2 to min
if first%m = 0 and numbers[n]%m = 0
flag1 = 0
exit
ok
next
if second < numbers[n]
min = second
else
min = numbers[n]
ok
for m = 2 to min
if second%m = 0 and numbers[n]%m = 0
flag2 = 0
exit
ok
next
if flag1 = 0 and flag2 = 1
see "" + numbers[n] + " "
first = second
second = numbers[n]
del(numbers,n)
row = row+1
if row%10 = 0
see nl
ok
num = num + 1
if num = 29
exit
ok
n = 3
ok
next
see "Found " + row + " Yellowstone numbers" + nl
see "done..." + nl- Output:
working... 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 Found 30 Yellowstone numbers done...
| Code | Comments |
|---|---|
≪ IF DUP2 < THEN SWAP END WHILE DUP REPEAT SWAP OVER MOD END DROP ≫ 'GCD' STO ≪ DUP SIZE 1 - GETI ROT ROT GET → am2 am1 ≪ 3 DO DO 1 + UNTIL DUP2 POS NOT END UNTIL DUP am1 GCD 1 == OVER am2 GCD 1 ≠ AND END + ≫ ≫ 'YELLO' STO |
( a b -- gcd(a,b) )
Ensure a > b
Euclidean algorithm
( { a(1)..a(n-1) } -- { a(1)..a(n) } )
Store locally a(n-1) and a(n-2)
Find smallest number not already in sequence
Check primality requirements
Add a(n) to the sequence
|
The following words in the command line deliver what is required:
{ 1 2 3 }
≪ 1 27 START YELLO NEXT ≫ EVAL
- Output:
1: { 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 ]
Recent RPL version
« DUP REVLIST 1 2 SUB -> a12
« 3
DO DO 1 + UNTIL DUP2 POS NOT END
UNTIL a12 OVER GCD EVAL 1 ≠ * 1 == END
+
» » 'YELLO' STO
« DO YELLO UNTIL DUP SIZE 30 ≥ END
» 'TASK' STO
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)
- 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]
// [dependencies]
// num = "0.3"
// plotters = "^0.2.15"
use num::integer::gcd;
use plotters::prelude::*;
use std::collections::HashSet;
fn yellowstone_sequence() -> impl std::iter::Iterator<Item = u32> {
let mut sequence: HashSet<u32> = HashSet::new();
let mut min = 1;
let mut n = 0;
let mut n1 = 0;
let mut n2 = 0;
std::iter::from_fn(move || {
n2 = n1;
n1 = n;
if n < 3 {
n += 1;
} else {
n = min;
while !(!sequence.contains(&n) && gcd(n1, n) == 1 && gcd(n2, n) > 1) {
n += 1;
}
}
sequence.insert(n);
while sequence.contains(&min) {
sequence.remove(&min);
min += 1;
}
Some(n)
})
}
// Based on the example in the "Quick Start" section of the README file for
// the plotters library.
fn plot_yellowstone(filename: &str) -> Result<(), Box<dyn std::error::Error>> {
let root = BitMapBackend::new(filename, (800, 600)).into_drawing_area();
root.fill(&WHITE)?;
let mut chart = ChartBuilder::on(&root)
.caption("Yellowstone Sequence", ("sans-serif", 24).into_font())
.margin(10)
.x_label_area_size(20)
.y_label_area_size(20)
.build_ranged(0usize..100usize, 0u32..180u32)?;
chart.configure_mesh().draw()?;
chart.draw_series(LineSeries::new(
yellowstone_sequence().take(100).enumerate(),
&BLUE,
))?;
Ok(())
}
fn main() {
println!("First 30 Yellowstone numbers:");
for y in yellowstone_sequence().take(30) {
print!("{} ", y);
}
println!();
match plot_yellowstone("yellowstone.png") {
Ok(()) => {}
Err(error) => eprintln!("Error: {}", error),
}
}
- Output:
A plot of the first 100 Yellowstone numbers is saved to the file "yellowstone.png".
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
Media:Yellowstone sequence rust.png
import scala.util.control.Breaks._
object YellowstoneSequence extends App {
println(s"First 30 values in the yellowstone sequence:\n${yellowstoneSequence(30)}")
def yellowstoneSequence(sequenceCount: Int): List[Int] = {
var yellowstoneList = List(1, 2, 3)
var num = 4
var notYellowstoneList = List[Int]()
while (yellowstoneList.size < sequenceCount) {
val foundIndex = notYellowstoneList.indexWhere(test =>
gcd(yellowstoneList(yellowstoneList.size - 2), test) > 1 &&
gcd(yellowstoneList.last, test) == 1
)
if (foundIndex >= 0) {
yellowstoneList = yellowstoneList :+ notYellowstoneList(foundIndex)
notYellowstoneList = notYellowstoneList.patch(foundIndex, Nil, 1)
} else {
breakable({
while (true) {
if (gcd(yellowstoneList(yellowstoneList.size - 2), num) > 1 &&
gcd(yellowstoneList.last, num) == 1) {
yellowstoneList = yellowstoneList :+ num
num += 1
// break the inner while loop
break
}
notYellowstoneList = notYellowstoneList :+ num
num += 1
}
});
}
}
yellowstoneList
}
def gcd(a: Int, b: Int): Int = if (b == 0) a else gcd(b, a % b)
}
- Output:
First 30 values in the yellowstone sequence: List(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)
proc gcd {a b} {
while {$b} {
lassign [list $b [expr {$a % $b}]] a b
}
return $a
}
proc gen_yellowstones {{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]
- 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
Dim @y(30)
@y(0) = 1
@y(1) = 2
@y(2) = 3
For i = 3 To 29
k = 3
Do
k = k + 1
If (FUNC(_gcd(k, @y(i-2))) = 1) + (FUNC(_gcd(k, @y(i-1))) > 1) Then
Continue
EndIf
For j = 0 To i - 1
If @y(j) = k Then Unloop : Continue
Next
@y(i) = k : Break
Loop
Next
For i = 0 To 29
Print @y(i); " ";
Next
Print : End
_gcd Param (2)
If b@ = 0 Then Return (a@)
Return (FUNC(_gcd(b@, a@ % b@)))
- 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 0 OK, 0:670
Very simple approach, probably lots of room for improvement. Run it in Uiua Pad to see the plot.
Gcd ← ⊙◌⍢(˜⟜◿|⊃⋅±⊙∘)
# Experimental!
Gcd ← ∨
NextY ← ⍢(+1|⨬(¬≍0_1=1≡Gcd⊙↙₋₂|1)˜∊⊃⊙∘⊙∘)1 # [a0 a1...an] -> an+1
N ← 200
⍢(˜⊂NextY|<N⧻)[1 2 3]
⟜↙₃₀ # First 30
▽⟜≡▽2⊞= ⇌⇡N # Plot 200- 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]
Function gcd(a As Long, b As Long) As Long
If b = 0 Then
gcd = a
Exit Function
End If
gcd = gcd(b, a Mod b)
End Function
Sub Yellowstone()
Dim i As Long, j As Long, k As Long, Y(1 To 30) As Long
Y(1) = 1
Y(2) = 2
Y(3) = 3
For i = 4 To 30
k = 3
Do
k = k + 1
If gcd(k, Y(i - 2)) = 1 Or gcd(k, Y(i - 1)) > 1 Then GoTo EndLoop:
For j = 1 To i - 1
If Y(j) = k Then GoTo EndLoop:
Next j
Y(i) = k
Exit Do
EndLoop:
Loop
Next i
For i = 1 To 30
Debug.Print Y(i) & " ";
Next i
End Sub
- 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
fn gcd(xx int, yy int) int {
mut x := xx
mut y := yy
for y != 0 {
x, y = y, x%y
}
return x
}
fn yellowstone(n int) []int {
mut m := map[int]bool{}
mut a := []int{len: n+1}
for i in 1..4 {
a[i] = i
m[i] = true
}
mut 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..]
}
fn main() {
mut x := []int{len: 100}
for i in 0..100 {
x[i] = i + 1
}
y := yellowstone(100)
println("The first 30 Yellowstone numbers are:")
println(y[..30])
}- 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]
Without the extra credit part.
import "./math" for Int
var yellowstone = Fn.new { |n|
var m = {}
var a = List.filled(n + 1, 0)
for (i in 1..3) {
a[i] = i
m[i] = true
}
var min = 4
for (c in 4..n) {
var i = min
while (true) {
if (!m[i] && Int.gcd(a[c-1], i) == 1 && Int.gcd(a[c-2], i) > 1) {
a[c] = i
m[i] = true
if (i == min) min = min + 1
break
}
i = i + 1
}
}
return a[1..-1]
}
var y = yellowstone.call(30)
System.print("The first 30 Yellowstone numbers are:")
System.print(y)
- 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]
func GCD(N, D); \Return the greatest common divisor of N and D
int N, D, R; \numerator, denominator, remainder
[if D > N then
[R:=D; D:=N; N:=R]; \swap D and N
while D > 0 do
[R:= rem(N/D);
N:= D;
D:= R;
];
return N;
];
int I, A(30+1), N, T;
[for I:= 1 to 3 do A(I):= I; \givens
N:= 4;
repeat T:= 4;
loop [if GCD(T, A(N-1)) = 1 and \relatively prime
GCD(T, A(N-2)) # 1 then \not relatively prime
[loop [for I:= 1 to N-1 do \test if in sequence
if T = A(I) then quit;
quit;
];
if I = N then \T is not in sequence so
[A(N):= T; \ add it in
N:= N+1;
quit;
];
];
T:= T+1; \next trial
];
until N > 30;
for N:= 1 to 30 do
[IntOut(0, A(N)); ChOut(0, ^ )];
\\for N:= 1 to 100 do Point(N, A(N)); \plot demonstration
]- 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
!ys-0
defn main(n=30):
say: n:yellow
defn yellow(n):
a =: +[1 2 3]
b =: (1 .. 3).zipmap([true]:cycle)
i =: 4
loop a a, b b, i i:
if n > a.#:
? if b.$i.! &&
(i.gcd(a.$) == 1) &&
(i.gcd(a:butlast.$) > 1)
: recur:
a.conj(i), b.assoc(i true), 5
recur: a, b, i.++
=>: a
defn gcd(a b):
loop a a, b b:
if b.?:
recur b: a % b
else: a
- Output:
$ ys yellowstone-sequence.ys [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 "std/set.zc"
fn gcd(x: int, y: int) -> int {
while y {
let t = y;
y = x % y;
x = t;
}
return x >= 0 ? x :-x;
}
fn yellowstone(n: int, a: int*) {
let s = Set<int>::new();
for i in 1..4 {
a[i] = i;
s.add(i);
}
let min = 4;
for c in 4..=n {
for let i = min; ; ++i {
if !s.contains(i) && gcd(a[c - 1], i) == 1 && gcd(a[c - 2], i) > 1 {
a[c] = i;
s.add(i);
if i == min { min++; }
break;
}
}
}
}
fn main() {
let x: int[30];
for i in 0..30 { x[i] = i + 1; }
let a: [int; 31];
yellowstone(30, (int*)a);
println "The first 30 Yellowstone numbers are:";
for i in 1..=30 {
print "{a[i]:2d} ";
if !(i % 10) { println ""; }
}
}
- 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
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
