duffinianNumbers: procedure options(main);
    %replace MAXSIGMA by 10000;
    declare sigma (1:MAXSIGMA) fixed;

    calculateSigmaTable: procedure;
        declare (i, j) fixed;
        do i=1 to MAXSIGMA;
            sigma(i) = 0;
        end;
        do i=1 to MAXSIGMA;
            do j=i to MAXSIGMA by i;
                sigma(j) = sigma(j) + i;
            end;
        end;
    end calculateSigmaTable;

    gcd: procedure(aa,bb) returns(fixed);
        declare (a, aa, b, bb, c) fixed;
        a = aa;
        b = bb;
        do while(b > 0);
            c = mod(a,b);
            a = b;
            b = c;
        end;
        return(a);
    end gcd;

    duffinian: procedure(n) returns(bit);
        declare n fixed;
        return(sigma(n) > n+1 & gcd(n, sigma(n)) = 1);
    end duffinian;

    triplet: procedure(n) returns(bit);
        declare n fixed;
        return(duffinian(n) & duffinian(n+1) & duffinian(n+2));
    end triplet;

    declare (i, n) fixed;

    call calculateSigmaTable;
    put skip list('First 50 Duffinian numbers:');
    put skip;
    n=0;
    do i=1 to 50;
        do n=n+1 repeat(n+1) while(^duffinian(n)); end;
        put edit(n) (F(5));
        if mod(i,10) = 0 then put skip;
    end;

    put skip;
    put skip list('First 15 Duffinian triplets:');
    n=0;
    do i=1 to 15;
        do n=n+1 repeat(n+1) while(^triplet(n)); end;
        put skip edit(n, n+1, n+2) (F(7),F(7),F(7));
    end;
end duffinianNumbers;