Duffinian numbers: Difference between revisions

From Rosetta Code
Content added Content deleted
(Added XPL0 example.)
(Add Factor)
Line 35: Line 35:





=={{header|Factor}}==
{{works with|Factor|0.99 2021-06-02}}
<lang 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.</lang>
{{out}}
<pre>
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
</pre>


=={{header|Julia}}==
=={{header|Julia}}==

Revision as of 18:01, 25 February 2022

Duffinian numbers is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.

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



Factor

Works with: Factor version 0.99 2021-06-02

<lang 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.</lang>

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

Julia

<lang 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()

</lang>

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

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]

Raku

<lang perl6>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";</lang>
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)

Wren

Library: Wren-math
Library: Wren-seq
Library: Wren-fmt

<lang ecmascript>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)</lang>

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

<lang 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;
       ];

]</lang>

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