program duffinian_numbers;
    init sigma := divisor_sum_table(20000);

    print("First 50 Duffinian numbers:");
    loop for n in first(50, routine is_duffinian) do
        nprint(lpad(str n, 6));
        if (col +:= 1) mod 10 = 0 then print; end if;
    end loop;

    print;
    print("First 15 Duffinian triplets:");
    loop for n in first(15, routine is_duffinian_triplet) do
        print(+/[lpad(str i, 6) : i in [n,n+1,n+2]]);
    end loop;

    proc first(num, pred);
        ls := [];
        loop while #ls < num do
            loop until call(pred, n) do
                n +:= 1;
            end loop;
            ls with:= n;
        end loop;
        return ls;
    end proc;

    proc is_duffinian_triplet(n);
        return and/[is_duffinian(i) : i in [n,n+1,n+2]];
    end proc;

    proc is_duffinian(n);
        return sigma(n) > n+1 and gcd(n, sigma(n)) = 1;
    end proc;

    proc gcd(a,b);
        return if b=0 then a else gcd(b, a mod b) end;
    end proc;

    proc divisor_sum_table(sz);
        ds := [0] * sz;
        loop for i in [1..sz] do
            loop for j in [i,i+i..sz] do
                ds(j) +:= i;
            end loop;
        end loop;
        return ds;
    end proc;
end program;