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.