Duffinian numbers

From Rosetta Code
Task
Duffinian numbers
You are encouraged to solve this task according to the task description, using any language you may know.

A Duffinian number is a composite number k that is relatively prime to its sigma sum σ.

The sigma sum of k is the sum of the divisors of k.


E.G.

161 is a Duffinian number.

  • It is composite. (7 × 23)
  • The sigma sum 192 (1 + 7 + 23 + 161) is relatively prime to 161.


Duffinian numbers are very common.

It is not uncommon for two consecutive integers to be Duffinian (a Duffinian twin) (8, 9), (35, 36), (49, 50), etc.

Less common are Duffinian triplets; three consecutive Duffinian numbers. (63, 64, 65), (323, 324, 325), etc.

Much, much less common are Duffinian quadruplets and quintuplets. The first Duffinian quintuplet is (202605639573839041, 202605639573839042, 202605639573839043, 202605639573839044, 202605639573839045).

It is not possible to have six consecutive Duffinian numbers


Task
  • Find and show the first 50 Duffinian numbers.
  • Find and show at least the first 15 Duffinian triplets.


See also



ALGOL 68

Constructs a table of divisor counts without doing any divisions.

BEGIN # find Duffinian numbers: non-primes relatively prime to their divisor count #
    INT max number            := 500 000;      # largest number we will consider   #
    # iterative Greatest Common Divisor routine, returns the gcd of m and n        #
    PROC gcd = ( INT m, n )INT:
         BEGIN
            INT a := ABS m, b := ABS n;
            WHILE b /= 0 DO
                INT new a = b;
                b        := a MOD b;
                a        := new a
            OD;
            a
         END # gcd # ;
    # construct a table of the divisor counts                                      #
    [ 1 : max number ]INT ds; FOR i TO UPB ds DO ds[ i ] := 1 OD;
    FOR i FROM 2 TO UPB ds
        DO FOR j FROM i BY i TO UPB ds DO ds[ j ] +:= i OD
    OD;
    # set the divisor counts of non-Duffinian numbers to 0                         #
    ds[ 1 ] := 0;  # 1 is not Duffinian                                            #
    FOR n FROM 2 TO UPB ds DO
        IF INT nds = ds[ n ];
           IF nds = n + 1 THEN TRUE ELSE gcd( n, nds ) /= 1 FI
        THEN
            # n is prime or is not relatively prime to its divisor sum             #
            ds[ n ] := 0
        FI
    OD;
    # show the first 50 Duffinian numbers                                          #
    print( ( "The first 50 Duffinian numbers:", newline ) );
    INT dcount := 0;
    FOR n WHILE dcount < 50 DO
        IF ds[ n ] /= 0
        THEN # found a Duffinian number                                            #
            print( ( " ", whole( n, -3) ) );
            IF ( dcount +:= 1 ) MOD 25 = 0 THEN print( ( newline ) ) FI
        FI
    OD;
    print( ( newline ) );
    # show the duffinian triplets below UPB ds                                     #
    print( ( "The Duffinian triplets up to ", whole( UPB ds, 0 ), ":", newline ) );
    dcount := 0;
    FOR n FROM 3 TO UPB ds DO
        IF ds[ n - 2 ] /= 0 AND ds[ n - 1 ] /= 0 AND ds[ n ] /= 0
        THEN # found a Duffinian triplet                                           #
            print( ( " (", whole( n - 2, -7 ), " ", whole( n - 1, -7 ), " ", whole( n, -7 ), ")" ) );
            IF ( dcount +:= 1 ) MOD 4 = 0 THEN print( ( newline ) ) FI
        FI
    OD
END
Output:
The first 50 Duffinian numbers:
   4   8   9  16  21  25  27  32  35  36  39  49  50  55  57  63  64  65  75  77  81  85  93  98 100
 111 115 119 121 125 128 129 133 143 144 155 161 169 171 175 183 185 187 189 201 203 205 209 215 217

The Duffinian triplets up to 500000:
 (     63      64      65) (    323     324     325) (    511     512     513) (    721     722     723)
 (    899     900     901) (   1443    1444    1445) (   2303    2304    2305) (   2449    2450    2451)
 (   3599    3600    3601) (   3871    3872    3873) (   5183    5184    5185) (   5617    5618    5619)
 (   6049    6050    6051) (   6399    6400    6401) (   8449    8450    8451) (  10081   10082   10083)
 (  10403   10404   10405) (  11663   11664   11665) (  12481   12482   12483) (  13447   13448   13449)
 (  13777   13778   13779) (  15841   15842   15843) (  17423   17424   17425) (  19043   19044   19045)
 (  26911   26912   26913) (  30275   30276   30277) (  36863   36864   36865) (  42631   42632   42633)
 (  46655   46656   46657) (  47523   47524   47525) (  53137   53138   53139) (  58563   58564   58565)
 (  72961   72962   72963) (  76175   76176   76177) (  79523   79524   79525) (  84099   84100   84101)
 (  86527   86528   86529) (  94177   94178   94179) ( 108899  108900  108901) ( 121103  121104  121105)
 ( 125315  125316  125317) ( 128017  128018  128019) ( 129599  129600  129601) ( 137287  137288  137289)
 ( 144399  144400  144401) ( 144721  144722  144723) ( 154567  154568  154569) ( 158403  158404  158405)
 ( 166463  166464  166465) ( 167041  167042  167043) ( 175231  175232  175233) ( 177607  177608  177609)
 ( 181475  181476  181477) ( 186623  186624  186625) ( 188497  188498  188499) ( 197191  197192  197193)
 ( 199711  199712  199713) ( 202499  202500  202501) ( 211249  211250  211251) ( 230399  230400  230401)
 ( 231199  231200  231201) ( 232561  232562  232563) ( 236195  236196  236197) ( 242063  242064  242065)
 ( 243601  243602  243603) ( 248003  248004  248005) ( 260099  260100  260101) ( 260641  260642  260643)
 ( 272483  272484  272485) ( 274575  274576  274577) ( 285155  285156  285157) ( 291599  291600  291601)
 ( 293763  293764  293765) ( 300303  300304  300305) ( 301087  301088  301089) ( 318095  318096  318097)
 ( 344449  344450  344451) ( 354481  354482  354483) ( 359551  359552  359553) ( 359999  360000  360001)
 ( 367235  367236  367237) ( 374543  374544  374545) ( 403201  403202  403203) ( 406801  406802  406803)
 ( 417697  417698  417699) ( 419903  419904  419905) ( 423199  423200  423201) ( 435599  435600  435601)
 ( 468511  468512  468513) ( 470449  470450  470451) ( 488071  488072  488073)

APL

duffinian_numbers{
    sigma  +/(0=⍳|⊢)
    duff  sigma((1=∨)∧⊣>1+⊢)
    'First 50 Duffinian numbers:'
    5 10((/⍨)duff¨)220
    'First 15 Duffinian triplets:'
    (0 1 2∘.+(/⍨)0 1 2(⊃∧.)(duff¨))8500
}
Output:
First 50 Duffinian numbers:
  4   8   9  16  21  25  27  32  35  36
 39  49  50  55  57  63  64  65  75  77
 81  85  93  98 100 111 115 119 121 125
128 129 133 143 144 155 161 169 171 175
183 185 187 189 201 203 205 209 215 217
First 15 Duffinian triplets:
  63   64   65
 323  324  325
 511  512  513
 721  722  723
 899  900  901
1443 1444 1445
2303 2304 2305
2449 2450 2451
3599 3600 3601
3871 3872 3873
5183 5184 5185
5617 5618 5619
6049 6050 6051
6399 6400 6401
8449 8450 8451

AppleScript

As is often the case with these tasks, it takes as much code to format the output as it does to get the numbers.  :)

on aliquotSum(n)
    if (n < 2) then return 0
    set sum to 1
    set sqrt to n ^ 0.5
    set limit to sqrt div 1
    if (limit = sqrt) then
        set sum to sum + limit
        set limit to limit - 1
    end if
    repeat with i from 2 to limit
        if (n mod i is 0) then set sum to sum + i + n div i
    end repeat
    
    return sum
end aliquotSum

on hcf(a, b)
    repeat until (b = 0)
        set x to a
        set a to b
        set b to x mod b
    end repeat
    
    if (a < 0) then return -a
    return a
end hcf

on isDuffinian(n)
    set aliquot to aliquotSum(n) -- = sigma sum - n. = 1 if n's prime.
    return ((aliquot > 1) and (hcf(n, aliquot + n) = 1))
end isDuffinian

-- Task code:
on matrixToText(matrix, w)
    script o
        property matrix : missing value
        property row : missing value
    end script
    
    set o's matrix to matrix
    set padding to "          "
    repeat with r from 1 to (count o's matrix)
        set o's row to o's matrix's item r
        repeat with i from 1 to (count o's row)
            set o's row's item i to text -w thru end of (padding & o's row's item i)
        end repeat
        set o's matrix's item r to join(o's row, "")
    end repeat
    
    return join(o's matrix, linefeed)
end matrixToText

on join(lst, delim)
    set astid to AppleScript's text item delimiters
    set AppleScript's text item delimiters to delim
    set txt to lst as text
    set AppleScript's text item delimiters to astid
    return txt
end join

on task(duffTarget, tupTarget, tupSize)
    if ((duffTarget < 1) or (tupTarget < 1) or (tupSize < 2)) then error "Duff parameter(s)."
    script o
        property duffinians : {}
        property tuplets : {}
    end script
    
    -- Populate o's duffinians and tuplets lists.
    set n to 1
    set tuplet to {}
    repeat while (((count o's tuplets) < tupTarget) or ((count o's duffinians) < duffTarget))
        if (isDuffinian(n)) then
            if ((count o's duffinians) < duffTarget) then set end of o's duffinians to n
            if (tuplet ends with n - 1) then
                set end of tuplet to n
            else
                if ((count tuplet) = tupSize) then set end of o's tuplets to tuplet
                set tuplet to {n}
            end if
        end if
        set n to n + 1
    end repeat
    
    -- Format for output.
    set duffinians to {}
    repeat with i from 1 to duffTarget by 20
        set j to i + 19
        if (j > duffTarget) then set j to duffTarget
        set end of duffinians to items i thru j of o's duffinians
    end repeat
    set part1 to "First " & duffTarget & " Duffinian numbers:" & linefeed & ¬
        matrixToText(duffinians, (count (end of o's duffinians as text)) + 2)
    set tupletTypes to {missing value, "twins", "triplets:", "quadruplets:", "quintuplets:"}
    set part2 to "First " & tupTarget & " Duffinian " & item tupSize of tupletTypes & linefeed & ¬
        matrixToText(o's tuplets, (count (end of end of o's tuplets as text)) + 2)
    
    return part1 & (linefeed & linefeed & part2)
end task

return task(50, 20, 3) -- First 50 Duffinians, first 20 3-item tuplets.
Output:
"First 50 Duffinian numbers:
    4    8    9   16   21   25   27   32   35   36   39   49   50   55   57   63   64   65   75   77
   81   85   93   98  100  111  115  119  121  125  128  129  133  143  144  155  161  169  171  175
  183  185  187  189  201  203  205  209  215  217

First 20 Duffinian triplets:
     63     64     65
    323    324    325
    511    512    513
    721    722    723
    899    900    901
   1443   1444   1445
   2303   2304   2305
   2449   2450   2451
   3599   3600   3601
   3871   3872   3873
   5183   5184   5185
   5617   5618   5619
   6049   6050   6051
   6399   6400   6401
   8449   8450   8451
  10081  10082  10083
  10403  10404  10405
  11663  11664  11665
  12481  12482  12483
  13447  13448  13449"

Arturo

duffinian?: function [n]->
    and? [not? prime? n]
         [
             fn: factors n
             [1] = intersection factors sum fn fn
         ]

first50: new []
i: 0
while [50 > size first50][
    if duffinian? i -> 'first50 ++ i
    i: i + 1
]

print "The first 50 Duffinian numbers:"
loop split.every: 10 first50 'row [
    print map to [:string] row 'item -> pad item 3
]

first15: new []
i: 0
while [15 > size first15][
    if every? i..i+2 => duffinian? [
        'first15 ++ @[@[i, i+1, i+2]]
        i: i+2
    ]
    i: i + 1
]

print ""
print "The first 15 Duffinian triplets:"
loop split.every: 5 first15 'row [
    print map row 'item -> pad.right as.code item 17
]
Output:
The first 50 Duffinian numbers:
  1   4   8   9  16  21  25  27  32  35 
 36  39  49  50  55  57  63  64  65  75 
 77  81  85  93  98 100 111 115 119 121 
125 128 129 133 143 144 155 161 169 171 
175 183 185 187 189 201 203 205 209 215 

The first 15 Duffinian triplets:
[63 64 65]        [323 324 325]     [511 512 513]     [721 722 723]     [899 900 901]     
[1443 1444 1445]  [2303 2304 2305]  [2449 2450 2451]  [3599 3600 3601]  [3871 3872 3873]  
[5183 5184 5185]  [5617 5618 5619]  [6049 6050 6051]  [6399 6400 6401]  [8449 8450 8451]

BASIC

Applesoft BASIC

Translation of: FreeBASIC
 100  DEF  FN MOD(NUM) = NUM -  INT (NUM / DIV) * DIV: REM NUM MOD DIV
 110  M = 50:N = 4
 120  PRINT "FIRST "M" DUFFINIAN NUMBERS:"
 130  FOR C = 0 TO M STEP 0
 140      GOSUB 600"DUFF
 150      IF DUFF THEN  PRINT  RIGHT$ ("   " +  STR$ (N),4);:C = C + 1
 160      N = N + 1
 170  NEXT C
 180  M = 15:S = 4:M$ =  CHR$ (13)
 190  PRINT M$M$"FIRST "M" DUFFINIAN TRIPLETS:"
 200  FOR C = 0 TO M STEP 0
 210      FOR D = 2 TO 0 STEP  - 1:N = S + D: GOSUB 600: IF DUFF THEN  NEXT D
 220      IF D < 0 THEN C = C + 1: PRINT  RIGHT$ ("   " +  STR$ (S) + "-",5) LEFT$ ( STR$ (S + 2) + "    ",5);:D = 0
 230      S = S + D + 1
 240  NEXT C
 250  END

 REM ISPRIME V RETURNS ISPRIME
 260  ISPRIME = FALSE: IF V < 2 THEN  RETURN
 270  DIV = 2:ISPRIME =  FN MOD(V): IF  NOT ISPRIME THEN ISPRIME = V = 2: RETURN
 280  LIMIT =  SQR (V): IF LIMIT >  = 3 THEN  FOR DIV = 3 TO LIMIT STEP 2:ISPRIME =  FN MOD(V): IF ISPRIME THEN  NEXT DIV
 290  RETURN

 REM GREATEST COMMON DIVISOR (GCD) A B RETURNS GCD
 300  A =  ABS ( INT (A))
 310  B =  ABS ( INT (B))
 320  GCD = A *  NOT  NOT B
 330  FOR B = B + A *  NOT B TO 0 STEP 0
 340      A = GCD
 350      GCD = B
 360      B = A -  INT (A / GCD) * GCD
 370  NEXT B
 380  RETURN

 REM SUMDIV NUM RETURNS SUM
 400  DIV = 2
 410  SUM = 0
 420  QUOT = NUM / DIV
 430  IF DIV > QUOT THEN SUM = 1: RETURN
 440  FOR LOOP = 0 TO 1 STEP 0
 450      IF  FN MOD(NUM) = 0 THEN SUM = SUM + DIV: IF DIV <  > QUOT THEN SUM = SUM + QUOT
 460      DIV = DIV + 1
 470      QUOT = NUM / DIV
 480      LOOP = DIV > QUOT
 500  NEXT LOOP
 510  SUM = SUM + 1
 520  RETURN

 REM DUFF N RETURNS DUFF
 600  DUFF = FALSE
 610  V = N: GOSUB 260"ISPRIME
 620  IF ISPRIME THEN  RETURN
 630  NUM = N: GOSUB 400"SUMDIV
 640  A = SUM:B = N: GOSUB 300"GCD
 650  DUFF = GCD = 1
 660  RETURN
Output:
FIRST 50 DUFFINIAN NUMBERS:
   4   8   9  16  21  25  27  32  35  36  39  49  50  55  57  63  64  65  75  77  81  85  93  98 100 111 115 119 121 125 128 129 133 143 144 155 161 169 171 175 183 185 187 189 201 203 205 209 215 217

FIRST 15 DUFFINIAN TRIPLETS:
  63-65    323-325   511-513   721-723   899-901  1443-1445 2303-2305 2449-2451 3599-3601 3871-3873 5183-5185 5617-5619 6049-6051 6399-6401 8449-8451 

FreeBASIC

Translation of: XPL0
#include "isprime.bas"

Function GCD(p As Integer, q As Integer) As Integer
    Return Iif(q = 0, p, GCD(q, p Mod q))
End Function

Function SumDiv(Num As Uinteger) As Uinteger
    Dim As Uinteger Div = 2, Sum = 0, Quot
    Do
        Quot = Num / Div
        If Div > Quot Then Exit Do
        If Num Mod Div = 0 Then
            Sum += Div
            If Div <> Quot Then Sum += Quot
        End If
        Div += 1
    Loop
    Return Sum+1
End Function

Function Duff(N As Uinteger) As Boolean
    Return Iif(isPrime(N), False, GCD(SumDiv(N), N) = 1)
End Function

Dim As Integer C = 0, N = 4
Print "First 50 Duffinian numbers:"
Do
    If Duff(N) Then
        Print Using "####"; N;
        C += 1
        If C Mod 20 = 0 Then Print
    End If
    N += 1
Loop Until C >= 50

C = 0 : N = 4
Print !"\n\nFirst 50 Duffinian triplets:"
Do
    If Duff(N) And Duff(N+1) And Duff(N+2) Then
        Print Using !" [###### ###### ######]\t"; N; N+1; N+2;
        C += 1
        If C Mod 4 = 0 Then Print
    End If
    N += 1
Loop Until C >= 50

Sleep
Output:
First 50 Duffinian numbers:
   4   8   9  16  21  25  27  32  35  36  39  49  50  55  57  63  64  65  75  77
  81  85  93  98 100 111 115 119 121 125 128 129 133 143 144 155 161 169 171 175
 183 185 187 189 201 203 205 209 215 217

First 50 Duffinian triplets:
 [    63     64     65]  [   323    324    325]  [   511    512    513]  [   721    722    723]
 [   899    900    901]  [  1443   1444   1445]  [  2303   2304   2305]  [  2449   2450   2451]
 [  3599   3600   3601]  [  3871   3872   3873]  [  5183   5184   5185]  [  5617   5618   5619]
 [  6049   6050   6051]  [  6399   6400   6401]  [  8449   8450   8451]  [ 10081  10082  10083]
 [ 10403  10404  10405]  [ 11663  11664  11665]  [ 12481  12482  12483]  [ 13447  13448  13449]
 [ 13777  13778  13779]  [ 15841  15842  15843]  [ 17423  17424  17425]  [ 19043  19044  19045]
 [ 26911  26912  26913]  [ 30275  30276  30277]  [ 36863  36864  36865]  [ 42631  42632  42633]
 [ 46655  46656  46657]  [ 47523  47524  47525]  [ 53137  53138  53139]  [ 58563  58564  58565]
 [ 72961  72962  72963]  [ 76175  76176  76177]  [ 79523  79524  79525]  [ 84099  84100  84101]
 [ 86527  86528  86529]  [ 94177  94178  94179]  [108899 108900 108901]  [121103 121104 121105]
 [125315 125316 125317]  [128017 128018 128019]  [129599 129600 129601]  [137287 137288 137289]
 [144399 144400 144401]  [144721 144722 144723]  [154567 154568 154569]  [158403 158404 158405]
 [166463 166464 166465]  [167041 167042 167043]

BCPL

get "libhdr"

let calcsigmas(sig, n) be
$(  sig!0 := 0
    for i = 0 to n do sig!i := 0
    for i = 1 to n/2 do
    $(  let j = i
        while 0 < j <= n do
        $(  sig!j := sig!j + i
            j := j + i
        $)
    $)
$)

let gcd(m, n) = n=0 -> m, gcd(n, m rem n)

let duff(sig, n) = sig!n > n+1 & gcd(n, sig!n) = 1
let triple(sig, n) = duff(sig, n) & duff(sig, n+1) & duff(sig, n+2)

let first(sig, f, max, cb) be
$(  let n = 0
    for i = 1 to max
    $(  n := n+1 repeatuntil f(sig, n)
        cb(i, n)
    $)
$)

let start() be
$(  let showsingle(i, n) be
    $(  writef("%I4", n)
        if i rem 10=0 then wrch('*N')
    $)

    let showtriple(i, n) be writef("%I2: %I6 %I6 %I6*N", i, n, n+1, n+2)

    let sig = getvec(20000)
    calcsigmas(sig, 20000)

    writes("First 50 Duffinian numbers:*N")
    first(sig, duff, 50, showsingle)

    writes("*NFirst 15 Duffinian triples:*N")
    first(sig, triple, 15, showtriple)
    freevec(sig)
$)
Output:
First 50 Duffinian numbers:
   4   8   9  16  21  25  27  32  35  36
  39  49  50  55  57  63  64  65  75  77
  81  85  93  98 100 111 115 119 121 125
 128 129 133 143 144 155 161 169 171 175
 183 185 187 189 201 203 205 209 215 217

First 15 Duffinian triples:
 1:     63     64     65
 2:    323    324    325
 3:    511    512    513
 4:    721    722    723
 5:    899    900    901
 6:   1443   1444   1445
 7:   2303   2304   2305
 8:   2449   2450   2451
 9:   3599   3600   3601
10:   3871   3872   3873
11:   5183   5184   5185
12:   5617   5618   5619
13:   6049   6050   6051
14:   6399   6400   6401
15:   8449   8450   8451

C++

#include <iomanip>
#include <iostream>
#include <numeric>
#include <sstream>

bool duffinian(int n) {
    if (n == 2)
        return false;
    int total = 1, power = 2, m = n;
    for (; (n & 1) == 0; power <<= 1, n >>= 1)
        total += power;
    for (int p = 3; p * p <= n; p += 2) {
        int sum = 1;
        for (power = p; n % p == 0; power *= p, n /= p)
            sum += power;
        total *= sum;
    }
    if (m == n)
        return false;
    if (n > 1)
        total *= n + 1;
    return std::gcd(total, m) == 1;
}

int main() {
    std::cout << "First 50 Duffinian numbers:\n";
    for (int n = 1, count = 0; count < 50; ++n) {
        if (duffinian(n))
            std::cout << std::setw(3) << n << (++count % 10 == 0 ? '\n' : ' ');
    }
    std::cout << "\nFirst 50 Duffinian triplets:\n";
    for (int n = 1, m = 0, count = 0; count < 50; ++n) {
        if (duffinian(n))
            ++m;
        else
            m = 0;
        if (m == 3) {
            std::ostringstream os;
            os << '(' << n - 2 << ", " << n - 1 << ", " << n << ')';
            std::cout << std::left << std::setw(24) << os.str()
                      << (++count % 3 == 0 ? '\n' : ' ');
        }
    }
    std::cout << '\n';
}
Output:
First 50 Duffinian numbers:
  4   8   9  16  21  25  27  32  35  36
 39  49  50  55  57  63  64  65  75  77
 81  85  93  98 100 111 115 119 121 125
128 129 133 143 144 155 161 169 171 175
183 185 187 189 201 203 205 209 215 217

First 50 Duffinian triplets:
(63, 64, 65)             (323, 324, 325)          (511, 512, 513)         
(721, 722, 723)          (899, 900, 901)          (1443, 1444, 1445)      
(2303, 2304, 2305)       (2449, 2450, 2451)       (3599, 3600, 3601)      
(3871, 3872, 3873)       (5183, 5184, 5185)       (5617, 5618, 5619)      
(6049, 6050, 6051)       (6399, 6400, 6401)       (8449, 8450, 8451)      
(10081, 10082, 10083)    (10403, 10404, 10405)    (11663, 11664, 11665)   
(12481, 12482, 12483)    (13447, 13448, 13449)    (13777, 13778, 13779)   
(15841, 15842, 15843)    (17423, 17424, 17425)    (19043, 19044, 19045)   
(26911, 26912, 26913)    (30275, 30276, 30277)    (36863, 36864, 36865)   
(42631, 42632, 42633)    (46655, 46656, 46657)    (47523, 47524, 47525)   
(53137, 53138, 53139)    (58563, 58564, 58565)    (72961, 72962, 72963)   
(76175, 76176, 76177)    (79523, 79524, 79525)    (84099, 84100, 84101)   
(86527, 86528, 86529)    (94177, 94178, 94179)    (108899, 108900, 108901)
(121103, 121104, 121105) (125315, 125316, 125317) (128017, 128018, 128019)
(129599, 129600, 129601) (137287, 137288, 137289) (144399, 144400, 144401)
(144721, 144722, 144723) (154567, 154568, 154569) (158403, 158404, 158405)
(166463, 166464, 166465) (167041, 167042, 167043) 

Delphi

Works with: Delphi version 6.0


{These subroutines would normally be in a library, but is included here for clarity}

function GetAllProperDivisors(N: Integer;var IA: TIntegerDynArray): integer;
{Make a list of all the "proper dividers" for N}
{Proper dividers are the of numbers the divide evenly into N}
var I: integer;
begin
SetLength(IA,0);
for I:=1 to N-1 do
 if (N mod I)=0 then
	begin
	SetLength(IA,Length(IA)+1);
	IA[High(IA)]:=I;
	end;
Result:=Length(IA);
end;


function GetAllDivisors(N: Integer;var IA: TIntegerDynArray): integer;
{Make a list of all the "proper dividers" for N, Plus N itself}
begin
Result:=GetAllProperDivisors(N,IA)+1;
SetLength(IA,Length(IA)+1);
IA[High(IA)]:=N;
end;



function IsDuffinianNumber(N: integer): boolean;
{Test number to see if it a Duffinian number}
var Facts1,Facts2: TIntegerDynArray;
var Sum,I,J: integer;
begin
Result:=False;
{Must be a composite number}
if IsPrime(N) then exit;
{Get all divisors}
GetAllDivisors(N,Facts1);
{Get sum of factors}
Sum:=0;
for I:=0 to High(Facts1) do
 Sum:=Sum+Facts1[I];
{Get all factor of Sum}
GetAllDivisors(Sum,Facts2);
{Test if the two number share any factors}
for I:=1 to High(Facts1) do
 for J:=1 to High(Facts2) do
	if Facts1[I]=Facts2[J] then exit;
{If not, they are relatively prime}
Result:=True;
end;

procedure ShowDuffinianNumbers(Memo: TMemo);
var N,Cnt,D1,D2,D3: integer;
var S: string;
begin
Cnt:=0;
S:='';
Memo.Lines.Add('First 50 Duffinian Numbers');
for N:=2 to high(integer) do
 if IsDuffinianNumber(N) then
	begin
	Inc(Cnt);
	S:=S+Format('%5d',[N]);
	if (Cnt mod 10)=0 then S:=S+CRLF;
	if Cnt>=50 then break;
	end;
Memo.Lines.Add(S);
D1:=0; D2:=-10; D3:=0;
S:=''; Cnt:=0;
Memo.Lines.Add('First 15 Duffinian Triples');
for N:=2 to high(integer) do
 if IsDuffinianNumber(N) then
	begin
	D1:=D2; D2:=D3;
	D3:=N;
	if ((D2-D1)=1) and ((D3-D2)=1) then
		begin
		Inc(Cnt);
		S:=S+Format('(%5d%5d%5d) ',[D1,D2,D3]);
		if (Cnt mod 3)=0 then S:=S+CRLF;
		if Cnt>=15 then break;
		end;
	end;
Memo.Lines.Add(S);
end;
Output:
First 50 Duffinian Numbers
    4    8    9   16   21   25   27   32   35   36
   39   49   50   55   57   63   64   65   75   77
   81   85   93   98  100  111  115  119  121  125
  128  129  133  143  144  155  161  169  171  175
  183  185  187  189  201  203  205  209  215  217

First 20 Duffinian Triples
(    63    64    65) (   323   324   325) (   511   512   513) 
(   721   722   723) (   899   900   901) (  1443  1444  1445) 
(  2303  2304  2305) (  2449  2450  2451) (  3599  3600  3601) 
(  3871  3872  3873) (  5183  5184  5185) (  5617  5618  5619) 
(  6049  6050  6051) (  6399  6400  6401) (  8449  8450  8451) 
( 10081 10082 10083) ( 10403 10404 10405) ( 11663 11664 11665) 
( 12481 12482 12483) ( 13447 13448 13449) 

Elapsed Time: 825.340 ms.


Draco

word MAXSIGMA = 10000;
[MAXSIGMA+1]word sigma;

proc calcsigma() void:
    word i, j;
    for i from 0 upto MAXSIGMA do sigma[i] := 0 od;
    for i from 1 upto MAXSIGMA do
        for j from i by i upto MAXSIGMA do
            sigma[j] := sigma[j] + i
        od
    od
corp

proc gcd(word a, b) word:
    word c;
    while b > 0 do
        c := a % b;
        a := b;
        b := c;
    od;
    a
corp

proc duff(word n) bool:
    sigma[n] > n+1 and gcd(n, sigma[n]) = 1
corp

proc triplet(word n) bool:
    duff(n) and duff(n+1) and duff(n+2)
corp

proc first(word n; proc(word n)bool pred; proc(word i,n)void cb) void:
    word i, cur;
    cur := 0;
    for i from 1 upto n do
        while cur := cur + 1; not pred(cur) do od;
        cb(i, cur)
    od
corp

proc tablenum(word i, n) void:
    write(n:5);
    if i%10 = 0 then writeln() fi
corp

proc tripletline(word i, n) void:
    writeln(i:2, ' ', n:6, n+1:6, n+2:6)
corp

proc main() void:
    calcsigma();
    writeln("First 50 Duffinian numbers:");
    first(50, duff, tablenum);
    writeln();

    writeln("First 15 Duffinian triplets:");
    first(15, triplet, tripletline)
corp
Output:
First 50 Duffinian numbers:
    4    8    9   16   21   25   27   32   35   36
   39   49   50   55   57   63   64   65   75   77
   81   85   93   98  100  111  115  119  121  125
  128  129  133  143  144  155  161  169  171  175
  183  185  187  189  201  203  205  209  215  217

First 15 Duffinian triplets:
 1     63    64    65
 2    323   324   325
 3    511   512   513
 4    721   722   723
 5    899   900   901
 6   1443  1444  1445
 7   2303  2304  2305
 8   2449  2450  2451
 9   3599  3600  3601
10   3871  3872  3873
11   5183  5184  5185
12   5617  5618  5619
13   6049  6050  6051
14   6399  6400  6401
15   8449  8450  8451

EasyLang

fastfunc isprim num .
   i = 2
   while i <= sqrt num
      if num mod i = 0
         return 0
      .
      i += 1
   .
   return 1
.
func gcd a b .
   while b <> 0
      h = b
      b = a mod b
      a = h
   .
   return a
.
func sumdiv num .
   d = 2
   repeat
      quot = num div d
      until d > quot
      if num mod d = 0
         sum += d
         if d <> quot
            sum += quot
         .
      .
      d += 1
   .
   return sum + 1
.
func isduff n .
   if isprim n = 0 and gcd sumdiv n n = 1
      return 1
   .
   return 0
.
proc duffs . .
   print "First 50 Duffinian numbers:"
   n = 4
   repeat
      if isduff n = 1
         write n & " "
         cnt += 1
      .
      until cnt = 50
      n += 1
   .
   cnt = 0
   n = 4
   print "\n\nFirst 15 Duffinian triplets:"
   repeat
      if isduff n = 1 and isduff (n + 1) = 1 and isduff (n + 2) = 1
         print n & " - " & n + 2
         cnt += 1
      .
      until cnt = 15
      n += 1
   .
.
duffs

Factor

Works with: Factor version 0.99 2021-06-02
USING: combinators.short-circuit.smart grouping io kernel lists
lists.lazy math math.primes math.primes.factors math.statistics
prettyprint sequences sequences.deep ;

: duffinian? ( n -- ? )
    { [ prime? not ] [ dup divisors sum simple-gcd 1 = ] } && ;

: duffinians ( -- list ) 3 lfrom [ duffinian? ] lfilter ;

: triples ( -- list )
    duffinians dup cdr dup cdr lzip lzip [ flatten ] lmap-lazy
    [ differences { 1 1 } = ] lfilter ;

"First 50 Duffinian numbers:" print
50 duffinians ltake list>array 10 group simple-table. nl

"First 15 Duffinian triplets:" print
15 triples ltake list>array simple-table.
Output:
First 50 Duffinian numbers:
4   8   9   16  21  25  27  32  35  36
39  49  50  55  57  63  64  65  75  77
81  85  93  98  100 111 115 119 121 125
128 129 133 143 144 155 161 169 171 175
183 185 187 189 201 203 205 209 215 217

First 15 Duffinian triplets:
63   64   65
323  324  325
511  512  513
721  722  723
899  900  901
1443 1444 1445
2303 2304 2305
2449 2450 2451
3599 3600 3601
3871 3872 3873
5183 5184 5185
5617 5618 5619
6049 6050 6051
6399 6400 6401
8449 8450 8451

FutureBasic

local fn IsPrime( n as NSUInteger ) as BOOL
  BOOL       isPrime = YES
  NSUInteger i
  
  if n < 2        then exit fn = NO
  if n = 2        then exit fn = YES
  if n mod 2 == 0 then exit fn = NO
  for i = 3 to int(n^.5) step 2
    if n mod i == 0 then exit fn = NO
  next
end fn = isPrime

local fn GCD( a as long, b as long ) as long
  long r
  
  if ( a == 0 ) then r = b else r = fn GCD( b mod a, a )
end fn = r

local fn SumDiv( num as NSUInteger ) as NSUInteger
  NSUInteger div = 2, sum = 0, quot, result
  
  while (1)
    quot = num / div
    if ( div > quot ) then result = 0 : exit while
    if ( num mod div == 0 )
      sum += div
      if ( div != quot ) then sum += quot
    end if
    div++
  wend
  result = sum + 1
end fn = result

local fn IsDuffinian( n as NSUInteger) as BOOL
  BOOL result = NO
  
  if ( fn IsPrime(n) == NO && fn GCD( fn SumDiv(n), n ) == 1 ) then exit fn = YES
end fn = result

local fn FindDuffinians
  long c = 0, n = 4
  
  print "First 50 Duffinian numbers:"
  do
    if ( fn IsDuffinian(n) )
      printf @"%4d \b", n
      c++
      if ( c mod 10 == 0 ) then print
    end if
    n++
  until ( c >= 50 )
  
  c = 0 : n = 4
  printf @"\n\nFirst 56 Duffinian triplets:"
  do
    if ( fn IsDuffinian(n) and fn IsDuffinian(n + 1) and fn IsDuffinian(n + 2) )
      printf @" [%6ld %6ld %6ld] \b", n, n+1, n+2
      c++
      if ( c mod 4 == 0 ) then print
    end if
    n++
  until ( c >= 56 )
end fn

CFTimeInterval t
t = fn CACurrentMediaTime
fn FindDuffinians
printf @"\nCompute time: %.3f ms",(fn CACurrentMediaTime-t)*1000

HandleEvents
Output:
First 50 Duffinian numbers:
  4   8   9  16  21  25  27  32  35  36 
 39  49  50  55  57  63  64  65  75  77 
 81  85  93  98 100 111 115 119 121 125 
128 129 133 143 144 155 161 169 171 175 
183 185 187 189 201 203 205 209 215 217 

First 56 Duffinian triplets:
[    63     64     65] [   323    324    325] [   511    512    513] [   721    722    723]
[   899    900    901] [  1443   1444   1445] [  2303   2304   2305] [  2449   2450   2451]
[  3599   3600   3601] [  3871   3872   3873] [  5183   5184   5185] [  5617   5618   5619]
[  6049   6050   6051] [  6399   6400   6401] [  8449   8450   8451] [ 10081  10082  10083]
[ 10403  10404  10405] [ 11663  11664  11665] [ 12481  12482  12483] [ 13447  13448  13449]
[ 13777  13778  13779] [ 15841  15842  15843] [ 17423  17424  17425] [ 19043  19044  19045]
[ 26911  26912  26913] [ 30275  30276  30277] [ 36863  36864  36865] [ 42631  42632  42633]
[ 46655  46656  46657] [ 47523  47524  47525] [ 53137  53138  53139] [ 58563  58564  58565]
[ 72961  72962  72963] [ 76175  76176  76177] [ 79523  79524  79525] [ 84099  84100  84101]
[ 86527  86528  86529] [ 94177  94178  94179] [108899 108900 108901] [121103 121104 121105]
[125315 125316 125317] [128017 128018 128019] [129599 129600 129601] [137287 137288 137289]
[144399 144400 144401] [144721 144722 144723] [154567 154568 154569] [158403 158404 158405]
[166463 166464 166465] [167041 167042 167043] [175231 175232 175233] [177607 177608 177609]
[181475 181476 181477] [186623 186624 186625] [188497 188498 188499] [197191 197192 197193]

Compute time: 2963.753 ms

Go

Translation of: Wren
Library: Go-rcu
package main

import (
    "fmt"
    "math"
    "rcu"
)

func isSquare(n int) bool {
    s := int(math.Sqrt(float64(n)))
    return s*s == n
}

func main() {
    limit := 200000 // say
    d := rcu.PrimeSieve(limit-1, true)
    d[1] = false
    for i := 2; i < limit; i++ {
        if !d[i] {
            continue
        }
        if i%2 == 0 && !isSquare(i) && !isSquare(i/2) {
            d[i] = false
            continue
        }
        sigmaSum := rcu.SumInts(rcu.Divisors(i))
        if rcu.Gcd(sigmaSum, i) != 1 {
            d[i] = false
        }
    }

    var duff []int
    for i := 1; i < len(d); i++ {
        if d[i] {
            duff = append(duff, i)
        }
    }
    fmt.Println("First 50 Duffinian numbers:")
    rcu.PrintTable(duff[0:50], 10, 3, false)

    var triplets [][3]int
    for i := 2; i < limit; i++ {
        if d[i] && d[i-1] && d[i-2] {
            triplets = append(triplets, [3]int{i - 2, i - 1, i})
        }
    }
    fmt.Println("\nFirst 56 Duffinian triplets:")
    for i := 0; i < 14; i++ {
        s := fmt.Sprintf("%6v", triplets[i*4:i*4+4])
        fmt.Println(s[1 : len(s)-1])
    }
}
Output:
First 50 Duffinian numbers:
  4   8   9  16  21  25  27  32  35  36 
 39  49  50  55  57  63  64  65  75  77 
 81  85  93  98 100 111 115 119 121 125 
128 129 133 143 144 155 161 169 171 175 
183 185 187 189 201 203 205 209 215 217 

First 56 Duffinian triplets:
[    63     64     65] [   323    324    325] [   511    512    513] [   721    722    723]
[   899    900    901] [  1443   1444   1445] [  2303   2304   2305] [  2449   2450   2451]
[  3599   3600   3601] [  3871   3872   3873] [  5183   5184   5185] [  5617   5618   5619]
[  6049   6050   6051] [  6399   6400   6401] [  8449   8450   8451] [ 10081  10082  10083]
[ 10403  10404  10405] [ 11663  11664  11665] [ 12481  12482  12483] [ 13447  13448  13449]
[ 13777  13778  13779] [ 15841  15842  15843] [ 17423  17424  17425] [ 19043  19044  19045]
[ 26911  26912  26913] [ 30275  30276  30277] [ 36863  36864  36865] [ 42631  42632  42633]
[ 46655  46656  46657] [ 47523  47524  47525] [ 53137  53138  53139] [ 58563  58564  58565]
[ 72961  72962  72963] [ 76175  76176  76177] [ 79523  79524  79525] [ 84099  84100  84101]
[ 86527  86528  86529] [ 94177  94178  94179] [108899 108900 108901] [121103 121104 121105]
[125315 125316 125317] [128017 128018 128019] [129599 129600 129601] [137287 137288 137289]
[144399 144400 144401] [144721 144722 144723] [154567 154568 154569] [158403 158404 158405]
[166463 166464 166465] [167041 167042 167043] [175231 175232 175233] [177607 177608 177609]
[181475 181476 181477] [186623 186624 186625] [188497 188498 188499] [197191 197192 197193]

J

Implementation:

sigmasum=: >:@#.~/.~&.q:
composite=: 1&< * 0 = 1&p: 
duffinian=: composite * 1 = ] +. sigmasum

Task examples:

   5 10$(#~ duffinian) 1+i.1000
  4   8   9  16  21  25  27  32  35  36
 39  49  50  55  57  63  64  65  75  77
 81  85  93  98 100 111 115 119 121 125
128 129 133 143 144 155 161 169 171 175
183 185 187 189 201 203 205 209 215 217
   (i.3)+/~15 {.(#~ 1 1 1 E. duffinian) 1+i.100000
  63   64   65
 323  324  325
 511  512  513
 721  722  723
 899  900  901
1443 1444 1445
2303 2304 2305
2449 2450 2451
3599 3600 3601
3871 3872 3873
5183 5184 5185
5617 5618 5619
6049 6050 6051
6399 6400 6401
8449 8450 8451

Java

import java.util.Arrays;

public final class DuffianNumbers {

	public static void main(String[] aArgs) {
		int[] duffians = createDuffians(11_000);
		
		System.out.println("The first 50 Duffinian numbers:");
		int count = 0;
		int n = 1;
		while ( count < 50 ) {
			if ( duffians[n] > 0 ) {
				System.out.print(String.format("%4d%s", n, ( ++count % 25 == 0 ? "\n" : "" )));
			}
			n += 1;				
		}
		System.out.println();
		
		System.out.println("The first 16 Duffinian triplets:");
		count = 0;
		n = 3;
		while( count < 16 ) {
		    if ( duffians[n - 2] > 0 && duffians[n - 1] > 0 && duffians[n] > 0 ) { 
		    	System.out.print(String.format("%22s%s", 
		    		"(" + ( n - 2 ) + ", " + ( n - 1 ) + ", " + n + ")", ( ++count % 4 == 0 ? "\n" : "" )));  
		    }
		    n += 1;
		}
		System.out.println();
	}
	
	private static int[] createDuffians(int aLimit) {
		// Create a list where list[i] is the divisor sum of i.
		int[] result = new int[aLimit];
		Arrays.fill(result, 1);
		for ( int i = 2; i < aLimit; i++ ) {
			for ( int j = i; j < aLimit; j += i ) {
				result[j] += i;
			}
		}		
		
		// Set the divisor sum of non-Duffinian numbers to 0.
		result[1] = 0; //  1 is not a Duffinian number.
		for ( int n = 2; n < aLimit; n++ ) {
			int resultN = result[n];
		    if ( resultN == n + 1 || gcd(n, resultN) != 1 ) {
		    	// n is prime, or it is not relatively prime to its divisor sum.
		    	result[n] = 0;
		    }
		}	
		return result;
	}
	
	private static int gcd(int aOne, int aTwo) {
		if ( aTwo == 0 ) {
			return aOne;
		}
		return gcd(aTwo, aOne % aTwo);
	}

}
Output:
The first 50 Duffinian numbers:
   4   8   9  16  21  25  27  32  35  36  39  49  50  55  57  63  64  65  75  77  81  85  93  98 100
 111 115 119 121 125 128 129 133 143 144 155 161 169 171 175 183 185 187 189 201 203 205 209 215 217

The first 16 Duffinian triplets:
          (63, 64, 65)       (323, 324, 325)       (511, 512, 513)       (721, 722, 723)
       (899, 900, 901)    (1443, 1444, 1445)    (2303, 2304, 2305)    (2449, 2450, 2451)
    (3599, 3600, 3601)    (3871, 3872, 3873)    (5183, 5184, 5185)    (5617, 5618, 5619)
    (6049, 6050, 6051)    (6399, 6400, 6401)    (8449, 8450, 8451) (10081, 10082, 10083)

jq

Works with: jq

The solution presented here follows the Wren and similar entries on this page in hard-coding the size of the prime sieve used to produce the answer; while this approach helps speed things up, it does assume some foreknowledge.

Preliminaries

def count(s): reduce s as $x (0; .+1);

def isSquare:
  (sqrt|floor) as $sqrt
  | . == ($sqrt | .*.);

# Input:  a positive integer
# Output: an array, $a, of length .+1 such that
#         $a[$i] is $i if $i is prime, and false otherwise.
def primeSieve:
  # erase(i) sets .[i*j] to false for integral j > 1
  def erase($i):
    if .[$i] then
      reduce (range(2*$i; length; $i)) as $j (.; .[$j] = false) 
    else .
    end;
  (. + 1) as $n
  | (($n|sqrt) / 2) as $s
  | [null, null, range(2; $n)]
  | reduce (2, 1 + (2 * range(1; $s))) as $i (.; erase($i));

def gcd(a; b):
  # subfunction expects [a,b] as input
  # i.e. a ~ .[0] and b ~ .[1]
  def rgcd: if .[1] == 0 then .[0]
         else [.[1], .[0] % .[1]] | rgcd
         end;
  [a,b] | rgcd;

# divisors as an unsorted stream (without calling sqrt)
def divisors:
  if . == 1 then 1
  else . as $n
  | label $out
  | range(1; $n) as $i
  | ($i * $i) as $i2
  | if $i2 > $n then break $out
    else if $i2 == $n
         then $i
         elif ($n % $i) == 0
         then $i, ($n/$i)
         else empty
	 end
    end
  end;

The Task

# emit an array such that .[i] is i if i is Duffinian and false otherwise
def duffinianArray($limit):
   ($limit | primeSieve |  map(not))
   | .[1] = false
   | reduce range(2; $limit) as $i (.;
       if (.[$i]|not) then .
       else if ($i % 2) == 0 and ($i|isSquare|not) and (($i/2)|isSquare|not)
            then .[$i] = false 
            else sum($i|divisors) as $sigmaSum
            | if gcd($sigmaSum; $i) != 1
	      then .[$i] = false
	      else .
	      end
	    end
        end  );

# Input: duffinianArray($limit)
# Output: an array of the corresponding Duffinians
def duffinians:
 . as $d
 | reduce range(1;length) as $i ([]; if $d[$i] then . + [$i] else . end);

# Input: duffinians
# Output: stream of triplets
def triplets:
  . as $d
  | range (2; length) as $i
  | select( $d[$i] and $d[$i-1] and $d[$i-2] )
  | [$i-2, $i-1, $i];

def withCount(s; $msg):
  foreach (s,null) as $x (0; .+1;
    if $x == null then "\($msg) \(.-1)" else $x end );

def lpad($len): tostring | ($len - length) as $l | (" " * $l)[:$l] + .;

# 167039 is the minimum integer that is sufficient to produce 50 triplets
duffinianArray(167039)
| "First 50 Duffinian numbers:",
  (duffinians[0:50] | _nwise(10) | map(lpad(4)) | join(" ") ),
  "\nFirst 50 Duffinian triplets:",
  withCount(limit(50;triplets); "\nNumber of triplets: ")
Output:
First 50 Duffinian numbers:
   4    8    9   16   21   25   27   32   35   36
  39   49   50   55   57   63   64   65   75   77
  81   85   93   98  100  111  115  119  121  125
 128  129  133  143  144  155  161  169  171  175
 183  185  187  189  201  203  205  209  215  217

First 50 Duffinian triplets:
[63,64,65]
[323,324,325]
[511,512,513]
[721,722,723]
[899,900,901]
[1443,1444,1445]
[2303,2304,2305]
[2449,2450,2451]
[3599,3600,3601]
[3871,3872,3873]
[5183,5184,5185]
[5617,5618,5619]
[6049,6050,6051]
[6399,6400,6401]
[8449,8450,8451]
[10081,10082,10083]
[10403,10404,10405]
[11663,11664,11665]
[12481,12482,12483]
[13447,13448,13449]
[13777,13778,13779]
[15841,15842,15843]
[17423,17424,17425]
[19043,19044,19045]
[26911,26912,26913]
[30275,30276,30277]
[36863,36864,36865]
[42631,42632,42633]
[46655,46656,46657]
[47523,47524,47525]
[53137,53138,53139]
[58563,58564,58565]
[72961,72962,72963]
[76175,76176,76177]
[79523,79524,79525]
[84099,84100,84101]
[86527,86528,86529]
[94177,94178,94179]
[108899,108900,108901]
[121103,121104,121105]
[125315,125316,125317]
[128017,128018,128019]
[129599,129600,129601]
[137287,137288,137289]
[144399,144400,144401]
[144721,144722,144723]
[154567,154568,154569]
[158403,158404,158405]
[166463,166464,166465]
[167041,167042,167043]

Number of triplets:  50

Julia

using Primes

function σ(n)
    f = [one(n)]
    for (p,e) in factor(n)
        f = reduce(vcat, [f*p^j for j in 1:e], init=f)
    end
    return sum(f)
end

isDuffinian(n) = !isprime(n) && gcd(n, σ(n)) == 1

function testDuffinians()
    println("First 50 Duffinian numbers:")
    foreach(p -> print(rpad(p[2], 4), p[1] % 25 == 0 ? "\n" : ""),
        enumerate(filter(isDuffinian, 2:217)))
    n, found = 2, 0
    println("\nFifteen Duffinian triplets:")
    while found < 15
        if isDuffinian(n) && isDuffinian(n + 1) && isDuffinian(n + 2)
            println(lpad(n, 6), lpad(n +1, 6), lpad(n + 2, 6))
            found += 1
        end
        n += 1
    end
end

testDuffinians()
Output:
First 50 Duffinian numbers:
4   8   9   16  21  25  27  32  35  36  39  49  50  55  57  63  64  65  75  77  81  85  93  98  100 
111 115 119 121 125 128 129 133 143 144 155 161 169 171 175 183 185 187 189 201 203 205 209 215 217

Fifteen Duffinian triplets:
    63    64    65
   323   324   325
   511   512   513
   721   722   723
   899   900   901
  1443  1444  1445
  2303  2304  2305
  2449  2450  2451
  3599  3600  3601
  3871  3872  3873
  5183  5184  5185
  5617  5618  5619
  6049  6050  6051
  6399  6400  6401
  8449  8450  8451

MAD

            NORMAL MODE IS INTEGER
            DIMENSION SIGMA(10000),OUTROW(10)

            INTERNAL FUNCTION(AA,BB)
            ENTRY TO GCD.
            A = AA
            B = BB
STEP        WHENEVER A.E.B, FUNCTION RETURN A
            WHENEVER A.G.B, A = A-B
            WHENEVER A.L.B, B = B-A
            TRANSFER TO STEP
            END OF FUNCTION

            INTERNAL FUNCTION(N)
            ENTRY TO DUFF.
            SIG = SIGMA(N)
            FUNCTION RETURN SIG.G.N+1 .AND. GCD.(N,SIG).E.1
            END OF FUNCTION

            INTERNAL FUNCTION(N)
            ENTRY TO TRIP.
            FUNCTION RETURN DUFF.(N) .AND.
          0      DUFF.(N+1) .AND. DUFF.(N+2)
            END OF FUNCTION

            THROUGH SZERO, FOR I=1, 1, I.G.10000
SZERO       SIGMA(I) = 0
            THROUGH SCALC, FOR I=1, 1, I.G.10000
            THROUGH SCALC, FOR J=I, I, J.G.10000
SCALC       SIGMA(J) = SIGMA(J) + I

            PRINT COMMENT $ FIRST 50 DUFFINIAN NUMBERS$
            CAND = 0
            THROUGH DUFROW, FOR R=0, 1, R.GE.5
            THROUGH DUFCOL, FOR C=0, 1, C.GE.10
SCHDUF      THROUGH SCHDUF, FOR CAND=CAND+1, 1, DUFF.(CAND)
DUFCOL      OUTROW(C) = CAND
DUFROW      PRINT FORMAT ROWFMT,OUTROW(0),OUTROW(1),OUTROW(2),
          0           OUTROW(3),OUTROW(4),OUTROW(5),OUTROW(6),
          1           OUTROW(7),OUTROW(8),OUTROW(9)

            PRINT COMMENT $ $
            PRINT COMMENT $ FIRST 15 DUFFINIAN TRIPLETS$
            CAND = 0
            THROUGH DUFTRI, FOR S=0, 1, S.GE.15
SCHTRP      THROUGH SCHTRP, FOR CAND=CAND+1, 1, TRIP.(CAND)
DUFTRI      PRINT FORMAT TRIFMT,CAND,CAND+1,CAND+2

            VECTOR VALUES ROWFMT = $10(I5)*$
            VECTOR VALUES TRIFMT = $3(I7)*$
            END OF PROGRAM
Output:
FIRST 50 DUFFINIAN NUMBERS
    4    8    9   16   21   25   27   32   35   36
   39   49   50   55   57   63   64   65   75   77
   81   85   93   98  100  111  115  119  121  125
  128  129  133  143  144  155  161  169  171  175
  183  185  187  189  201  203  205  209  215  217

FIRST 15 DUFFINIAN TRIPLETS
     63     64     65
    323    324    325
    511    512    513
    721    722    723
    899    900    901
   1443   1444   1445
   2303   2304   2305
   2449   2450   2451
   3599   3600   3601
   3871   3872   3873
   5183   5184   5185
   5617   5618   5619
   6049   6050   6051
   6399   6400   6401
   8449   8450   8451

Mathematica /Wolfram Language

ClearAll[DuffianQ]
DuffianQ[n_Integer] := CompositeQ[n] \[And] CoprimeQ[DivisorSigma[1, n], n]
dns = Select[DuffianQ][Range[1000000]];
Take[dns, UpTo[50]]
triplets = ToString[dns[[#]]] <> "\[LongDash]" <> ToString[dns[[# + 2]]] & /@ SequencePosition[Differences[dns], {1, 1}][[All, 1]]
Multicolumn[triplets, {Automatic, 5}, Appearance -> "Horizontal"]
Output:

First 50 Duffinian numbers and Duffinian triplets below a million:

{4, 8, 9, 16, 21, 25, 27, 32, 35, 36, 39, 49, 50, 55, 57, 63, 64, 65, 75, 77, 81, 85, 93, 98, 100, 111, 115, 119, 121, 125, 128, 129, 133, 143, 144, 155, 161, 169, 171, 175, 183, 185, 187, 189, 201, 203, 205, 209, 215, 217}

63-65	323-325	511-513	721-723	899-901
1443-1445	2303-2305	2449-2451	3599-3601	3871-3873
5183-5185	5617-5619	6049-6051	6399-6401	8449-8451
10081-10083	10403-10405	11663-11665	12481-12483	13447-13449
13777-13779	15841-15843	17423-17425	19043-19045	26911-26913
30275-30277	36863-36865	42631-42633	46655-46657	47523-47525
53137-53139	58563-58565	72961-72963	76175-76177	79523-79525
84099-84101	86527-86529	94177-94179	108899-108901	121103-121105
125315-125317	128017-128019	129599-129601	137287-137289	144399-144401
144721-144723	154567-154569	158403-158405	166463-166465	167041-167043
175231-175233	177607-177609	181475-181477	186623-186625	188497-188499
197191-197193	199711-199713	202499-202501	211249-211251	230399-230401
231199-231201	232561-232563	236195-236197	242063-242065	243601-243603
248003-248005	260099-260101	260641-260643	272483-272485	274575-274577
285155-285157	291599-291601	293763-293765	300303-300305	301087-301089
318095-318097	344449-344451	354481-354483	359551-359553	359999-360001
367235-367237	374543-374545	403201-403203	406801-406803	417697-417699
419903-419905	423199-423201	435599-435601	468511-468513	470449-470451
488071-488073	504099-504101	506017-506019	518399-518401	521283-521285
522241-522243	529983-529985	547057-547059	585361-585363	589823-589825
617795-617797	640711-640713	647521-647523	656099-656101	659343-659345
675683-675685	682111-682113	685583-685585	688899-688901	700927-700929
703297-703299	710431-710433	725903-725905	746641-746643	751537-751539
756899-756901	791281-791283	798847-798849	809999-810001	814087-814089
834631-834633	837217-837219	842401-842403	842723-842725	857475-857477
860671-860673	910115-910117	913951-913953	963271-963273	968255-968257
991231-991233

Maxima

/* 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);
Output:
[4,8,9,16,21,25,27,32,35,36,39,49,50,55,57,63,64,65,75,77,81,85,93,98,100,111,115,119,121,125,128,129,133,143,144,155,161,169,171,175,183,185,187,189,201,203,205,209,215,217]

[[63,64,65],[323,324,325],[511,512,513],[721,722,723],[899,900,901],[1443,1444,1445],[2303,2304,2305],[2449,2450,2451],[3599,3600,3601],[3871,3872,3873],[5183,5184,5185],[5617,5618,5619],[6049,6050,6051],[6399,6400,6401],[8449,8450,8451]]

Modula-2

MODULE DuffinianNumbers;
FROM InOut IMPORT WriteCard, WriteString, WriteLn;

CONST
    MaxSigma = 10000;

VAR
    seen, cur: CARDINAL;
    sigma: ARRAY [1..MaxSigma] OF CARDINAL;

PROCEDURE CalculateSigmaTable;
    VAR i, j: CARDINAL;
BEGIN
    FOR i := 1 TO MaxSigma DO
        sigma[i] := 0
    END;
    FOR i := 1 TO MaxSigma DO
        j := i;
        WHILE j <= MaxSigma DO
            INC(sigma[j], i);
            INC(j, i);
        END
    END
END CalculateSigmaTable;

PROCEDURE GCD(a, b: CARDINAL): CARDINAL;
    VAR c: CARDINAL;
BEGIN
    WHILE b # 0 DO
        c := a MOD b;
        a := b;
        b := c
    END;
    RETURN a
END GCD;

PROCEDURE IsDuffinian(n: CARDINAL): BOOLEAN;
BEGIN
    RETURN (sigma[n] > n+1) AND (GCD(n, sigma[n]) = 1)
END IsDuffinian;

PROCEDURE IsDuffinianTriple(n: CARDINAL): BOOLEAN;
BEGIN
    RETURN IsDuffinian(n) AND IsDuffinian(n+1) AND IsDuffinian(n+2)
END IsDuffinianTriple;

BEGIN
    CalculateSigmaTable;
    WriteString("First 50 Duffinian numbers:");
    WriteLn;
    cur := 0;
    FOR seen := 1 TO 50 DO
        REPEAT INC(cur) UNTIL IsDuffinian(cur);
        WriteCard(cur, 4);
        IF seen MOD 10 = 0 THEN WriteLn END
    END;

    WriteLn;
    WriteString("First 15 Duffinian triples:");
    WriteLn;
    cur := 0;
    FOR seen := 1 TO 15 DO
        REPEAT INC(cur) UNTIL IsDuffinianTriple(cur);
        WriteCard(cur, 6);
        WriteCard(cur+1, 6);
        WriteCard(cur+2, 6);
        WriteLn
    END
END DuffinianNumbers.
Output:
First 50 Duffinian numbers:
   4   8   9  16  21  25  27  32  35  36
  39  49  50  55  57  63  64  65  75  77
  81  85  93  98 100 111 115 119 121 125
 128 129 133 143 144 155 161 169 171 175
 183 185 187 189 201 203 205 209 215 217

First 15 Duffinian triples:
    63    64    65
   323   324   325
   511   512   513
   721   722   723
   899   900   901
  1443  1444  1445
  2303  2304  2305
  2449  2450  2451
  3599  3600  3601
  3871  3872  3873
  5183  5184  5185
  5617  5618  5619
  6049  6050  6051
  6399  6400  6401
  8449  8450  8451

Nim

Translation of: ALGOL 68
import std/[algorithm, math, strformat]

const MaxNumber = 500_000

# Construct a table of the divisor counts.
var ds: array[1..MaxNumber, int]
ds.fill 1
for i in 2..MaxNumber:
  for j in countup(i, MaxNumber, i):
    ds[j] += i

# Set the divisor counts of non-Duffinian numbers to 0.
ds[1] = 0   # 1 is not Duffinian.
for n in 2..MaxNumber:
  let nds = ds[n]
  if nds == n + 1 or gcd(n, nds) != 1:
    # "n" is prime or is not relatively prime to its divisor sum.
    ds[n] = 0

# Show the first 50 Duffinian numbers.
echo "First 50 Duffinian numbers:"
var dcount = 0
var n = 1
while dcount < 50:
  if ds[n] != 0:
    stdout.write &" {n:3}"
    inc dcount
    if dcount mod 25 == 0:
      echo()
  inc n
echo()

# Show the Duffinian triplets below MaxNumber.
echo &"The Duffinian triplets up to {MaxNumber}:"
dcount = 0
for n in 3..MaxNumber:
  if ds[n - 2] != 0 and ds[n - 1] != 0 and ds[n] != 0:
    inc dcount
    stdout.write &" {(n - 2, n - 1, n): ^24}"
    stdout.write if dcount mod 4 == 0: '\n' else: ' '
echo()
Output:

The output is identical to that of the Algol 68 program, but the formatting is different.

First 50 Duffinian numbers:
   4   8   9  16  21  25  27  32  35  36  39  49  50  55  57  63  64  65  75  77  81  85  93  98 100
 111 115 119 121 125 128 129 133 143 144 155 161 169 171 175 183 185 187 189 201 203 205 209 215 217

The Duffinian triplets up to 500000:
       (63, 64, 65)            (323, 324, 325)           (511, 512, 513)           (721, 722, 723)     
     (899, 900, 901)          (1443, 1444, 1445)        (2303, 2304, 2305)        (2449, 2450, 2451)   
    (3599, 3600, 3601)        (3871, 3872, 3873)        (5183, 5184, 5185)        (5617, 5618, 5619)   
    (6049, 6050, 6051)        (6399, 6400, 6401)        (8449, 8450, 8451)      (10081, 10082, 10083)  
  (10403, 10404, 10405)     (11663, 11664, 11665)     (12481, 12482, 12483)     (13447, 13448, 13449)  
  (13777, 13778, 13779)     (15841, 15842, 15843)     (17423, 17424, 17425)     (19043, 19044, 19045)  
  (26911, 26912, 26913)     (30275, 30276, 30277)     (36863, 36864, 36865)     (42631, 42632, 42633)  
  (46655, 46656, 46657)     (47523, 47524, 47525)     (53137, 53138, 53139)     (58563, 58564, 58565)  
  (72961, 72962, 72963)     (76175, 76176, 76177)     (79523, 79524, 79525)     (84099, 84100, 84101)  
  (86527, 86528, 86529)     (94177, 94178, 94179)    (108899, 108900, 108901)  (121103, 121104, 121105)
 (125315, 125316, 125317)  (128017, 128018, 128019)  (129599, 129600, 129601)  (137287, 137288, 137289)
 (144399, 144400, 144401)  (144721, 144722, 144723)  (154567, 154568, 154569)  (158403, 158404, 158405)
 (166463, 166464, 166465)  (167041, 167042, 167043)  (175231, 175232, 175233)  (177607, 177608, 177609)
 (181475, 181476, 181477)  (186623, 186624, 186625)  (188497, 188498, 188499)  (197191, 197192, 197193)
 (199711, 199712, 199713)  (202499, 202500, 202501)  (211249, 211250, 211251)  (230399, 230400, 230401)
 (231199, 231200, 231201)  (232561, 232562, 232563)  (236195, 236196, 236197)  (242063, 242064, 242065)
 (243601, 243602, 243603)  (248003, 248004, 248005)  (260099, 260100, 260101)  (260641, 260642, 260643)
 (272483, 272484, 272485)  (274575, 274576, 274577)  (285155, 285156, 285157)  (291599, 291600, 291601)
 (293763, 293764, 293765)  (300303, 300304, 300305)  (301087, 301088, 301089)  (318095, 318096, 318097)
 (344449, 344450, 344451)  (354481, 354482, 354483)  (359551, 359552, 359553)  (359999, 360000, 360001)
 (367235, 367236, 367237)  (374543, 374544, 374545)  (403201, 403202, 403203)  (406801, 406802, 406803)
 (417697, 417698, 417699)  (419903, 419904, 419905)  (423199, 423200, 423201)  (435599, 435600, 435601)
 (468511, 468512, 468513)  (470449, 470450, 470451)  (488071, 488072, 488073) 

PARI/GP

Translation of: Julia
isDuffinian(n) = (!isprime(n)) && (gcd(n, sigma(n)) == 1);

testDuffinians()=
{
    print("First 50 Duffinian numbers:");
    count = 0; n = 2;
    while(count < 50, 
        if (isDuffinian(n),
            print1(n, " ");
            count++;
        );
        n++;
    );

    print("\n\nFifteen Duffinian triplets:");
    count = 0; n = 2;
    while (count < 15,
        if (isDuffinian(n) && isDuffinian(n + 1) && isDuffinian(n + 2),      
            print(n, " ", n + 1, " ", n + 2);
            count++;
        );
        n++;
    );
}

testDuffinians();
Output:
First 50 Duffinian numbers:
4 8 9 16 21 25 27 32 35 36 39 49 50 55 57 63 64 65 75 77 81 85 93 98 100 111 115 119 121 125 128 129 133 143 144 155 161 169 171 175 183 185 187 189 201 203 205 209 215 217 

Fifteen Duffinian triplets:
63 64 65
323 324 325
511 512 513
721 722 723
899 900 901
1443 1444 1445
2303 2304 2305
2449 2450 2451
3599 3600 3601
3871 3872 3873
5183 5184 5185
5617 5618 5619
6049 6050 6051
6399 6400 6401
8449 8450 8451

PascalABC.NET

uses school;

function is_duffinian(x: integer) := 
    (gcd(x, divisors(x).Sum) = 1) and 
    (divisors(x).ToArray.Length > 2);

begin
  var count := 0;
  var i := 0;
  writeln('First 50 Duffinian numbers:');
  while count < 50 do
  begin
    if is_duffinian(i) then 
    begin
      write(i:4);
      count += 1;
      if count mod 10 = 0 then writeln;
    end;
    i += 1
  end;
  
  count := 0;
  i := 0;
  writeln;
  writeln('First 15 Duffinian triplets:');
  while count < 15 do
  begin
    if is_duffinian(i) and is_duffinian(i + 1) and is_duffinian(i + 2) then 
    begin
      writeln(i:6, i + 1:6, i + 2:6);
      count += 1;
      i += 3;
    end;
    i += 1;
  end;
end.
Output:
First 50 Duffinian numbers:
   4   8   9  16  21  25  27  32  35  36
  39  49  50  55  57  63  64  65  75  77
  81  85  93  98 100 111 115 119 121 125
 128 129 133 143 144 155 161 169 171 175
 183 185 187 189 201 203 205 209 215 217

First 15 Duffinian triplets:
    63    64    65
   323   324   325
   511   512   513
   721   722   723
   899   900   901
  1443  1444  1445
  2303  2304  2305
  2449  2450  2451
  3599  3600  3601
  3871  3872  3873
  5183  5184  5185
  5617  5618  5619
  6049  6050  6051
  6399  6400  6401
  8449  8450  8451

Perl

Library: ntheory
use strict;
use warnings;
use feature <say state>;
use List::Util 'max';
use ntheory qw<divisor_sum is_prime gcd>;

sub table { my $t = shift() * (my $c = 1 + max map {length} @_); ( sprintf( ('%'.$c.'s')x@_, @_) ) =~ s/.{1,$t}\K/\n/gr }

sub duffinian {
    my($n) = @_;
    state $c = 1; state @D;
    do { push @D, $c if ! is_prime ++$c and 1 == gcd($c,divisor_sum($c)) } until @D > $n;
    $D[$n];
}

say "First 50 Duffinian numbers:";
say table 10, map { duffinian $_-1 } 1..50;

my(@d3,@triples) = (4, 8, 9); my $n = 3;
while (@triples < 39) {
    push @triples, '('.join(', ',@d3).')' if $d3[1] == 1+$d3[0] and $d3[2] == 2+$d3[0];
    shift @d3 and push @d3, duffinian ++$n;
}

say 'First 39 Duffinian triplets:';
say table 3, @triples;
Output:
First 50 Duffinian numbers:
   4   8   9  16  21  25  27  32  35  36
  39  49  50  55  57  63  64  65  75  77
  81  85  93  98 100 111 115 119 121 125
 128 129 133 143 144 155 161 169 171 175
 183 185 187 189 201 203 205 209 215 217

First 39 Duffinian triplets:
             (63, 64, 65)          (323, 324, 325)          (511, 512, 513)
          (721, 722, 723)          (899, 900, 901)       (1443, 1444, 1445)
       (2303, 2304, 2305)       (2449, 2450, 2451)       (3599, 3600, 3601)
       (3871, 3872, 3873)       (5183, 5184, 5185)       (5617, 5618, 5619)
       (6049, 6050, 6051)       (6399, 6400, 6401)       (8449, 8450, 8451)
    (10081, 10082, 10083)    (10403, 10404, 10405)    (11663, 11664, 11665)
    (12481, 12482, 12483)    (13447, 13448, 13449)    (13777, 13778, 13779)
    (15841, 15842, 15843)    (17423, 17424, 17425)    (19043, 19044, 19045)
    (26911, 26912, 26913)    (30275, 30276, 30277)    (36863, 36864, 36865)
    (42631, 42632, 42633)    (46655, 46656, 46657)    (47523, 47524, 47525)
    (53137, 53138, 53139)    (58563, 58564, 58565)    (72961, 72962, 72963)
    (76175, 76176, 76177)    (79523, 79524, 79525)    (84099, 84100, 84101)
    (86527, 86528, 86529)    (94177, 94178, 94179) (108899, 108900, 108901)

Phix

with javascript_semantics
sequence duffinian = {false}
integer n = 2, count = 0, triplet = 0, triple_count = 0
while triple_count<50 do
    bool bDuff = not is_prime(n) and gcd(n,sum(factors(n,1)))=1
    duffinian &= bDuff
    if bDuff then
        count += 1
        if count=50 then
            sequence s50 = apply(true,sprintf,{{"%3d"},find_all(true,duffinian)})
            printf(1,"First 50 Duffinian numbers:\n%s\n",join_by(s50,1,25," "))
        end if
        triplet += 1
        triple_count += (triplet>=3)
    else
        triplet = 0
    end if
    n += 1
end while
sequence s = apply(true,sq_add,{match_all({true,true,true},duffinian),{{0,1,2}}}),
         p = apply(true,pad_tail,{apply(true,sprintf,{{"[%d,%d,%d]"},s}),24})
printf(1,"First 50 Duffinian triplets:\n%s\n",{join_by(p,1,4," ")})
Output:
First 50 Duffinian numbers:
  4   8   9  16  21  25  27  32  35  36  39  49  50  55  57  63  64  65  75  77  81  85  93  98 100
111 115 119 121 125 128 129 133 143 144 155 161 169 171 175 183 185 187 189 201 203 205 209 215 217

First 50 Duffinian triplets:
[63,64,65]               [323,324,325]            [511,512,513]            [721,722,723]
[899,900,901]            [1443,1444,1445]         [2303,2304,2305]         [2449,2450,2451]
[3599,3600,3601]         [3871,3872,3873]         [5183,5184,5185]         [5617,5618,5619]
[6049,6050,6051]         [6399,6400,6401]         [8449,8450,8451]         [10081,10082,10083]
[10403,10404,10405]      [11663,11664,11665]      [12481,12482,12483]      [13447,13448,13449]
[13777,13778,13779]      [15841,15842,15843]      [17423,17424,17425]      [19043,19044,19045]
[26911,26912,26913]      [30275,30276,30277]      [36863,36864,36865]      [42631,42632,42633]
[46655,46656,46657]      [47523,47524,47525]      [53137,53138,53139]      [58563,58564,58565]
[72961,72962,72963]      [76175,76176,76177]      [79523,79524,79525]      [84099,84100,84101]
[86527,86528,86529]      [94177,94178,94179]      [108899,108900,108901]   [121103,121104,121105]
[125315,125316,125317]   [128017,128018,128019]   [129599,129600,129601]   [137287,137288,137289]
[144399,144400,144401]   [144721,144722,144723]   [154567,154568,154569]   [158403,158404,158405]
[166463,166464,166465]   [167041,167042,167043]

PL/I

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;
Output:
First 50 Duffinian numbers:
    4    8    9   16   21   25   27   32   35   36
   39   49   50   55   57   63   64   65   75   77
   81   85   93   98  100  111  115  119  121  125
  128  129  133  143  144  155  161  169  171  175
  183  185  187  189  201  203  205  209  215  217


First 15 Duffinian triplets:
     63     64     65
    323    324    325
    511    512    513
    721    722    723
    899    900    901
   1443   1444   1445
   2303   2304   2305
   2449   2450   2451
   3599   3600   3601
   3871   3872   3873
   5183   5184   5185
   5617   5618   5619
   6049   6050   6051
   6399   6400   6401
   8449   8450   8451

PL/M

100H:
BDOS: PROCEDURE (F,A); DECLARE F BYTE, A ADDRESS; GO TO 5; END BDOS;
EXIT: PROCEDURE; GO TO 0; END EXIT;
PR$CHAR: PROCEDURE (C); DECLARE C BYTE; CALL BDOS(2,C); END PR$CHAR;
PRINT: PROCEDURE (S); DECLARE S ADDRESS; CALL BDOS(9,S); END PRINT;

PR$NUM: PROCEDURE (N, WIDTH);
    DECLARE N ADDRESS, WIDTH BYTE;
    DECLARE S (6) BYTE INITIAL ('.....$');
    DECLARE P ADDRESS, DG BASED P BYTE;
    P = .S(5);
DIGIT:
    P = P - 1;
    DG = '0' + N MOD 10;
    IF WIDTH > 0 THEN WIDTH = WIDTH - 1;
    IF (N := N / 10) > 0 THEN GO TO DIGIT;
    CALL PRINT(P);
    DO WHILE WIDTH > 0;
        CALL PR$CHAR(' ');
        WIDTH = WIDTH - 1;
    END;
END PR$NUM;

DECLARE MAX$SIGMA LITERALLY '10$001';
DECLARE SIGMA (MAX$SIGMA) ADDRESS;
CALC$SIGMA: PROCEDURE;
    DECLARE (I, J) ADDRESS;
    DO I = 1 TO MAX$SIGMA-1;
        SIGMA(I) = 0;
    END;
    DO I = 1 TO MAX$SIGMA-1;
        DO J = I TO MAX$SIGMA-1 BY I;
            SIGMA(J) = SIGMA(J) + I;
        END;
    END;
END CALC$SIGMA;

GCD: PROCEDURE (X, Y) ADDRESS;
    DECLARE (X, Y, Z) ADDRESS;
    DO WHILE Y > 0;
        Z = X MOD Y;
        X = Y;
        Y = Z;
    END;
    RETURN X;
END GCD;

DUFF: PROCEDURE (N) BYTE;
    DECLARE N ADDRESS;
    RETURN SIGMA(N) > N+1 AND GCD(N, SIGMA(N)) = 1;
END DUFF;

DUFF$TRIPLE: PROCEDURE (N) BYTE;
    DECLARE N ADDRESS;
    RETURN DUFF(N) AND DUFF(N+1) AND DUFF(N+2);
END DUFF$TRIPLE;

DECLARE N ADDRESS, I BYTE;

CALL CALC$SIGMA;
CALL PRINT(.('FIRST 50 DUFFINIAN NUMBERS:',13,10,'$'));
N = 0;
DO I = 1 TO 50;
    DO WHILE NOT DUFF(N := N+1); END;
    CALL PR$NUM(N, 4);
    IF I MOD 10 = 0 THEN CALL PRINT(.(13,10,'$'));
END;

CALL PRINT(.(13,10,'FIRST 15 DUFFINIAN TRIPLES:',13,10,'$'));
N = 0;
DO I = 1 TO 15;
    DO WHILE NOT DUFF$TRIPLE(N := N+1); END;
    CALL PR$NUM(N, 6);
    CALL PR$NUM(N+1, 6);
    CALL PR$NUM(N+2, 6);
    CALL PRINT(.(13,10,'$'));
END;

CALL EXIT;
EOF
Output:
FIRST 50 DUFFINIAN NUMBERS:
4   8   9   16  21  25  27  32  35  36
39  49  50  55  57  63  64  65  75  77
81  85  93  98  100 111 115 119 121 125
128 129 133 143 144 155 161 169 171 175
183 185 187 189 201 203 205 209 215 217

FIRST 15 DUFFINIAN TRIPLES:
63    64    65
323   324   325
511   512   513
721   722   723
899   900   901
1443  1444  1445
2303  2304  2305
2449  2450  2451
3599  3600  3601
3871  3872  3873
5183  5184  5185
5617  5618  5619
6049  6050  6051
6399  6400  6401
8449  8450  8451

Python

# duffinian.py by xing216

def factors(n):
    factors = []
    for i in range(1, n + 1):
       if n % i == 0:
           factors.append(i)
    return factors
def gcd(a, b):
    while b != 0:
        a, b = b, a % b
    return a
is_relively_prime = lambda a, b: gcd(a, b) == 1
sigma_sum = lambda x: sum(factors(x))
is_duffinian = lambda x: is_relively_prime(x, sigma_sum(x)) and len(factors(x)) > 2
count = 0
i = 0
while count < 50:
    if is_duffinian(i):
        print(i, end=' ')
        count += 1
    i+=1
count2 = 0
j = 0
while count2 < 20:
    if is_duffinian(j) and is_duffinian(j+1) and is_duffinian(j+2):
        print(f"({j},{j+1},{j+2})", end=' ')
        count2 += 1
        j+=3
    j+=1

Quackery

factors is defined at Factors of an integer#Quackery.

gcd is defined at Greatest common divisor#Quackery.

  [ dup factors
    dup size 3 < iff
      [ 2drop false ] done
    0 swap witheach +
    gcd 1 = ]              is duffinian ( n --> b )

  [] 0
  [ dup duffinian if
      [ tuck join swap ]
    1+
    over size 50 = until ]
  drop
  [] swap
  witheach
    [ number$ nested join ]
  60 wrap$
  cr cr
  0 temp put
  [] 0
  [ dup duffinian iff
      [ 1 temp tally ]
    else
      [ 0 temp replace ]
    temp share 2 > if
      [ tuck 2 -
        join swap ]
    1+
    over size 15 = until ]
  drop
  [] swap
  witheach
    [ dup 1+ dup 1+
      join join
      nested join ]
  witheach [ echo cr ]
Output:
4 8 9 16 21 25 27 32 35 36 39 49 50 55 57 63 64 65 75 77 81
85 93 98 100 111 115 119 121 125 128 129 133 143 144 155 161
169 171 175 183 185 187 189 201 203 205 209 215 217

[ 63 64 65 ]
[ 323 324 325 ]
[ 511 512 513 ]
[ 721 722 723 ]
[ 899 900 901 ]
[ 1443 1444 1445 ]
[ 2303 2304 2305 ]
[ 2449 2450 2451 ]
[ 3599 3600 3601 ]
[ 3871 3872 3873 ]
[ 5183 5184 5185 ]
[ 5617 5618 5619 ]
[ 6049 6050 6051 ]
[ 6399 6400 6401 ]
[ 8449 8450 8451 ]

Raku

use Prime::Factor;

my @duffinians = lazy (3..*).hyper.grep: { !.is-prime && $_ gcd .&divisors.sum == 1 };

put "First 50 Duffinian numbers:\n" ~
@duffinians[^50].batch(10)».fmt("%3d").join: "\n";

put "\nFirst 40 Duffinian triplets:\n" ~
    ((^∞).grep: -> $n { (@duffinians[$n] + 1 == @duffinians[$n + 1]) && (@duffinians[$n] + 2 == @duffinians[$n + 2]) })[^40]\
    .map( { "({@duffinians[$_ .. $_+2].join: ', '})" } ).batch(4)».fmt("%-24s").join: "\n";
Output:
First 50 Duffinian numbers:
  4   8   9  16  21  25  27  32  35  36
 39  49  50  55  57  63  64  65  75  77
 81  85  93  98 100 111 115 119 121 125
128 129 133 143 144 155 161 169 171 175
183 185 187 189 201 203 205 209 215 217

First 40 Duffinian triplets:
(63, 64, 65)             (323, 324, 325)          (511, 512, 513)          (721, 722, 723)         
(899, 900, 901)          (1443, 1444, 1445)       (2303, 2304, 2305)       (2449, 2450, 2451)      
(3599, 3600, 3601)       (3871, 3872, 3873)       (5183, 5184, 5185)       (5617, 5618, 5619)      
(6049, 6050, 6051)       (6399, 6400, 6401)       (8449, 8450, 8451)       (10081, 10082, 10083)   
(10403, 10404, 10405)    (11663, 11664, 11665)    (12481, 12482, 12483)    (13447, 13448, 13449)   
(13777, 13778, 13779)    (15841, 15842, 15843)    (17423, 17424, 17425)    (19043, 19044, 19045)   
(26911, 26912, 26913)    (30275, 30276, 30277)    (36863, 36864, 36865)    (42631, 42632, 42633)   
(46655, 46656, 46657)    (47523, 47524, 47525)    (53137, 53138, 53139)    (58563, 58564, 58565)   
(72961, 72962, 72963)    (76175, 76176, 76177)    (79523, 79524, 79525)    (84099, 84100, 84101)   
(86527, 86528, 86529)    (94177, 94178, 94179)    (108899, 108900, 108901) (121103, 121104, 121105)

RPL

∑DIV, which returns the sum of the divisors of a given number, is defined at Sum_of_divisors#RPL.

GCD, which returns the GCD of 2 given numbers, is defined at Greatest_common_divisor#RPL.

Works with: Halcyon Calc version 4.2.8
RPL code Comment
≪ 
   DUP ∑DIV
   IF DUP2 1 - == THEN DROP2 0 ELSE GCD 1 == END 
≫ ‘DUFF?’ STO

≪ { } 2 
   WHILE OVER SIZE 50 < REPEAT 
     IF DUP DUFF? THEN SWAP OVER + SWAP END 
     1 + 
   END DROP 
≫ ‘TASK1’ STO

≪ { } 4 → duff n 
  ≪ 0 0 0 
     WHILE duff SIZE 15 ≤ REPEAT 
       ROT DROP n DUFF?
       IF 3 DUPN + + 3 == THEN 
           n 2 - n 1 - n 3 →ARRY 
           duff SWAP + 'duff' STO END 
        n 1 + 'n' STO END 
     3 DROPN duff
≫ ≫ ‘TASK2' STO
DUFF? ( n -- boolean ) 
get σ
if composite then check gcd(n,σ)


TASK1 ( -- { duff1..duff50 } ) 
loop from n=2 until 50 items in list
   if n is Duffinian then append to list
   n += 1
clean stack


TASK2 ( -- { [duff_triplets] } ) 
put  3 'false' boolean values in stack
loop from n=4 until 15 items in list
   update stack with n Duffinian status
   if first 3 stack levels are all 1
      create triplet
      append it to list
   n += 1
clean stack, display list

Output:
2: { 4 8 9 16 21 25 27 32 35 36 39 49 50 55 57 63 64 65 75 77 81 85 93 98 100 111 115 119 121 125 128 129 133 143 144 155 161 169 171 175 183 185 187 189 201 203 205 209 215 217 }
1: { [ 63 64 65 ] [ 323 324 325 ] [ 511 512 513 ] [ 721 722 723 ] [ 899 900 901 ] [ 1443 1444 1445 ] [ 2303 2304 2305 ] [ 2449 2450 2451 ] [ 3599 3600 3601 ] [ 3871 3872 3873 ] [ 5183 5184 5185 ] [ 5617 5618 5619 ] [ 6049 6050 6051 ] [ 6399 6400 6401 ] [ 8449 8450 8451 ] }

Ruby

require "prime"

class Integer

  def proper_divisors(prim_div = prime_division)
    return [] if self == 1
    primes = prim_div.flat_map{|prime, freq| [prime] * freq}
    (1...primes.size).each_with_object([1]) do |n, res|
      primes.combination(n).map{|combi| res << combi.inject(:*)}
    end.flatten.uniq
  end

  def duffinian?
    pd = prime_division
    return false if pd.sum(&:last) < 2
    gcd(proper_divisors(pd).sum + self) == 1
  end

end

n = 50
puts "The first #{n} Duffinian numbers:"
(1..).lazy.select(&:duffinian?).first(n).each_slice(10) do |slice|
  puts "%4d" * slice.size % slice
end

puts "\nThe first #{n} Duffinian triplets:"
(1..).each_cons(3).lazy.select{|slice| slice.all?(&:duffinian?)}.first(n).each do |group|
  puts "%8d" * group.size % group
end
Output:
The first 50 Duffinian numbers:
   4   8   9  16  21  25  27  32  35  36
  39  49  50  55  57  63  64  65  75  77
  81  85  93  98 100 111 115 119 121 125
 128 129 133 143 144 155 161 169 171 175
 183 185 187 189 201 203 205 209 215 217

The first 50 Duffinian triplets:
      63      64      65
     323     324     325
     511     512     513
     721     722     723
     899     900     901
    1443    1444    1445
    2303    2304    2305
    2449    2450    2451
    3599    3600    3601
    3871    3872    3873
    5183    5184    5185
    5617    5618    5619
    6049    6050    6051
    6399    6400    6401
    8449    8450    8451
   10081   10082   10083
   10403   10404   10405
   11663   11664   11665
   12481   12482   12483
   13447   13448   13449
   13777   13778   13779
   15841   15842   15843
   17423   17424   17425
   19043   19044   19045
   26911   26912   26913
   30275   30276   30277
   36863   36864   36865
   42631   42632   42633
   46655   46656   46657
   47523   47524   47525
   53137   53138   53139
   58563   58564   58565
   72961   72962   72963
   76175   76176   76177
   79523   79524   79525
   84099   84100   84101
   86527   86528   86529
   94177   94178   94179
  108899  108900  108901
  121103  121104  121105
  125315  125316  125317
  128017  128018  128019
  129599  129600  129601
  137287  137288  137289
  144399  144400  144401
  144721  144722  144723
  154567  154568  154569
  158403  158404  158405
  166463  166464  166465
  167041  167042  167043

Sidef

func is_duffinian(n) {
    n.is_composite && n.is_coprime(n.sigma)
}

say "First 50 Duffinian numbers:"
say 50.by(is_duffinian)

say "\nFirst 15 Duffinian triplets:"
15.by{|n| ^3 -> all {|k| is_duffinian(n+k) } }.each {|n|
    printf("(%s, %s, %s)\n", n, n+1, n+2)
}
Output:
First 50 Duffinian numbers:
[4, 8, 9, 16, 21, 25, 27, 32, 35, 36, 39, 49, 50, 55, 57, 63, 64, 65, 75, 77, 81, 85, 93, 98, 100, 111, 115, 119, 121, 125, 128, 129, 133, 143, 144, 155, 161, 169, 171, 175, 183, 185, 187, 189, 201, 203, 205, 209, 215, 217]

First 15 Duffinian triplets:
(63, 64, 65)
(323, 324, 325)
(511, 512, 513)
(721, 722, 723)
(899, 900, 901)
(1443, 1444, 1445)
(2303, 2304, 2305)
(2449, 2450, 2451)
(3599, 3600, 3601)
(3871, 3872, 3873)
(5183, 5184, 5185)
(5617, 5618, 5619)
(6049, 6050, 6051)
(6399, 6400, 6401)
(8449, 8450, 8451)

Wren

Library: Wren-math
Library: Wren-fmt
import "./math" for Int
import "./fmt" for Fmt

var limit = 200000  // say
var d = Int.primeSieve(limit-1, false)
d[1] = false
for (i in 2...limit) {
    if (!d[i]) continue
    if (i % 2 == 0 && !Int.isSquare(i) && !Int.isSquare(i/2)) {
        d[i] = false
        continue
    }
    var sigmaSum = Int.divisorSum(i)
    if (Int.gcd(sigmaSum, i) != 1) d[i] = false
}

var duff = (1...d.count).where { |i| d[i] }.toList
System.print("First 50 Duffinian numbers:")
Fmt.tprint("$3d", duff[0..49], 10)

var triplets = []
for (i in 2...limit) {
    if (d[i] && d[i-1] && d[i-2]) triplets.add([i-2, i-1, i])
}
System.print("\nFirst 50 Duffinian triplets:")
Fmt.tprint("$-25n", triplets[0..49], 4)
Output:
First 50 Duffinian numbers:
  4   8   9  16  21  25  27  32  35  36
 39  49  50  55  57  63  64  65  75  77
 81  85  93  98 100 111 115 119 121 125
128 129 133 143 144 155 161 169 171 175
183 185 187 189 201 203 205 209 215 217

First 50 Duffinian triplets:
[63, 64, 65]              [323, 324, 325]           [511, 512, 513]           [721, 722, 723]          
[899, 900, 901]           [1443, 1444, 1445]        [2303, 2304, 2305]        [2449, 2450, 2451]       
[3599, 3600, 3601]        [3871, 3872, 3873]        [5183, 5184, 5185]        [5617, 5618, 5619]       
[6049, 6050, 6051]        [6399, 6400, 6401]        [8449, 8450, 8451]        [10081, 10082, 10083]    
[10403, 10404, 10405]     [11663, 11664, 11665]     [12481, 12482, 12483]     [13447, 13448, 13449]    
[13777, 13778, 13779]     [15841, 15842, 15843]     [17423, 17424, 17425]     [19043, 19044, 19045]    
[26911, 26912, 26913]     [30275, 30276, 30277]     [36863, 36864, 36865]     [42631, 42632, 42633]    
[46655, 46656, 46657]     [47523, 47524, 47525]     [53137, 53138, 53139]     [58563, 58564, 58565]    
[72961, 72962, 72963]     [76175, 76176, 76177]     [79523, 79524, 79525]     [84099, 84100, 84101]    
[86527, 86528, 86529]     [94177, 94178, 94179]     [108899, 108900, 108901]  [121103, 121104, 121105] 
[125315, 125316, 125317]  [128017, 128018, 128019]  [129599, 129600, 129601]  [137287, 137288, 137289] 
[144399, 144400, 144401]  [144721, 144722, 144723]  [154567, 154568, 154569]  [158403, 158404, 158405] 
[166463, 166464, 166465]  [167041, 167042, 167043] 

XPL0

func IsPrime(N);        \Return 'true' if N is prime
int  N, I;
[if N <= 2 then return N = 2;
if (N&1) = 0 then \even >2\ return false;
for I:= 3 to sqrt(N) do
    [if rem(N/I) = 0 then return false;
    I:= I+1;
    ];
return true;
];

func    SumDiv(Num);    \Return sum of proper divisors of Num
int     Num, Div, Sum, Quot;
[Div:= 2;
Sum:= 0;
loop    [Quot:= Num/Div;
        if Div > Quot then quit;
        if rem(0) = 0 then
            [Sum:= Sum + Div;
            if Div # Quot then Sum:= Sum + Quot;
            ];
        Div:= Div+1;
        ];
return Sum+1;
];

func GCD(A, B);         \Return greatest common divisor of A and B
int  A, B;
[while A#B do
  if A>B then A:= A-B
         else B:= B-A;
return A;
];

func Duff(N);           \Return 'true' if N is a Duffinian number
int  N;
[if IsPrime(N) then return false;
return GCD(SumDiv(N), N) = 1;
];

int C, N;
[Format(4, 0);
C:= 0;  N:= 4;
loop    [if Duff(N) then
            [RlOut(0, float(N));
            C:= C+1;
            if C >= 50 then quit;
            if rem(C/20) = 0 then CrLf(0);
            ];
        N:= N+1;
        ];
CrLf(0);  CrLf(0);
Format(5, 0);
C:= 0;  N:= 4;
loop    [if Duff(N) & Duff(N+1) & Duff(N+2) then
            [RlOut(0, float(N));  RlOut(0, float(N+1));  RlOut(0, float(N+2));
            CrLf(0);
            C:= C+1;
            if C >= 15 then quit;
            ];
        N:= N+1;
        ];
]
Output:
   4   8   9  16  21  25  27  32  35  36  39  49  50  55  57  63  64  65  75  77
  81  85  93  98 100 111 115 119 121 125 128 129 133 143 144 155 161 169 171 175
 183 185 187 189 201 203 205 209 215 217

   63   64   65
  323  324  325
  511  512  513
  721  722  723
  899  900  901
 1443 1444 1445
 2303 2304 2305
 2449 2450 2451
 3599 3600 3601
 3871 3872 3873
 5183 5184 5185
 5617 5618 5619
 6049 6050 6051
 6399 6400 6401
 8449 8450 8451