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
 Numbers Aplenty - Duffinian numbers
 OEIS:A003624 - Duffinian numbers: composite numbers k relatively prime to sigma(k)




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