I'm working on modernizing Rosetta Code's infrastructure. Starting with communications. Please accept this time-limited open invite to RC's Slack.. --Michael Mol (talk) 20:59, 30 May 2020 (UTC)

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[edit]

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)

AppleScript[edit]

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[edit]

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]

C++[edit]

#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) 

Factor[edit]

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

Go[edit]

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[edit]

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

Julia[edit]

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

Perl[edit]

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[edit]

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]

Raku[edit]

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)

Sidef[edit]

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[edit]

Library: Wren-math
Library: Wren-seq
Library: Wren-fmt
import "./math" for Int, Nums
import "./seq" for Lst
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 = Nums.sum(Int.divisors(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:")
for (chunk in Lst.chunks(duff[0..49], 10)) Fmt.print("$3d", chunk)
 
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:")
for (chunk in Lst.chunks(triplets[0..49], 4)) Fmt.print("$-25s", chunk)
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[edit]

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