procedure main(A)
   write("The first 50 Duffian numbers:")
   every writes(" ",isduffinian(seq())\50)
   write("\n")
   limit := 15
   write("The first ",limit," Duffian triples:")
   every (isduffinian(n := seq()),isduffinian(n+1),isduffinian(n+2))\limit do
      write("\t(",n,",",n+1,",",n+2,")")
end

procedure isduffinian(n)
   x := \cfact(n)**\afact(sumfactors(n))
   return (*\x = 1, !x = 1, n)
end

procedure cfact(n)  # all factors of n if n is a composite number
   return (*(f := afact(n)) > 2, f)
end

procedure afact(n)  # all factors of n
   every (f := set(), i := 1 to n, n%i = 0,insert(f,i))
   return f
end

procedure sumfactors(n)
   every (s := 0, i := 1 to n, n%i = 0) do s +:= i
   return s
end