Duffinian numbers: Difference between revisions

=={{header|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'''

<syntaxhighlight lang=jq>

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

def isSquare:

(sqrt|floor) as $sqrt

| . == ($sqrt | .*.);

# return an array, $a, of length . or slightly greater, 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; (1 + length) / i) as $j (.; .[i * $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;

</syntaxhighlight>

'''The Task'''

<syntaxhighlight lang=jq>

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

</syntaxhighlight>

{{output}}

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

</pre>