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)
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]
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)
Factor
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
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
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
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
Nim
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)
Perl
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]
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.
| 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 ] }
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
import "./math" for Int
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 = 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