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