/* Predicate functions that checks wether an integer is a Duffinian number or not */
duffinianp(n):=if n#1 and not primep(n) and gcd(n,divsum(n))=1 then true$

/* Function that returns a list of the first len Duffinian numbers */
duffinian_count(len):=block(
    [i:1,count:0,result:[]],
    while count<len do (if duffinianp(i) then (result:endcons(i,result),count:count+1),i:i+1),
    result)$

/* Function that returns a list of the first len Duffinian triples */
duffinian_triples_count(len):=block(
    [i:1,count:0,result:[]],
    while count<len do (if map(duffinianp,[i,i+1,i+2])=[true,true,true] then (result:endcons([i,i+1,i+2],result),count:count+1),i:i+1),
    result)$

/* Test cases */
/* First 50 Duffinian numbers */
duffinian_count(50);

/* First 15 Duffinian triples */
duffinian_triples_count(15);