sternBrocot: procedure options(main);
    %replace MAX by 1200;
    declare S(1:MAX) fixed;

    /* find the first occurrence of N in S */
    findFirst: procedure(n) returns(fixed);
        declare (n, i) fixed;
        do i=1 to MAX;
            if S(i)=n then return(i);
        end;
    end findFirst;

    /* find the greatest common divisor of A and B */
    gcd: procedure(a, b) returns(fixed) recursive;
        declare (a, b) fixed;
        if b = 0 then return(a);
        return(gcd(b, mod(a, b)));
    end gcd;

    /* calculate S(i) up to MAX */
    declare i fixed;
    S(1) = 1; S(2) = 1;
    do i=2 to MAX/2;
        S(i*2-1) = S(i) + S(i-1);
        S(i*2) = S(i);
    end;

    /* print first 15 elements */
    put skip list('First 15 elements: ');
    do i=1 to 15;
        put edit(S(i)) (F(2));
    end;

    /* find first occurrences of 1..10 and 100 */
    do i=1 to 10;
        put skip list('First',i,'at',findFirst(i));
    end;
    put skip list('First   ',100,'at',findFirst(100));

    /* check GCDs of adjacent pairs up to 1000th element */
    do i=2 to 1000;
        if gcd(S(i-1),S(i)) ^= 1 then do;
            put skip list('GCD of adjacent pair not 1 at i=',i);
            stop;
        end;
    end;
    put skip list('All GCDs of adjacent pairs are 1.');
end sternBrocot;