Factor-perfect numbers: Difference between revisions
m (→{{header|Phix}}: oops, comparing/predicting 9 against 8) |
(→{{header|Wren}}: Now uses Phix's approach for last part - considerable speed-up.) |
||
Line 308: | Line 308: | ||
{{libheader|Wren-math}} |
{{libheader|Wren-math}} |
||
{{libheader|Wren-fmt}} |
{{libheader|Wren-fmt}} |
||
However, using Phix's 'g' function for the last part which is considerably quicker than generating a bunch of lists and then counting them - overall time down from 76 to 7 seconds. Still not really worth attempting stretch goal though. |
|||
A bit slow on the last part - 76 seconds compared to 46 seconds for Python - so I haven't attempted the stretch goal. |
|||
<syntaxhighlight lang="ecmascript">import "./math" for Int |
<syntaxhighlight lang="ecmascript">import "./math" for Int, Nums |
||
import "./fmt" for Fmt |
import "./fmt" for Fmt |
||
Line 328: | Line 328: | ||
// Get factorization sequence count. |
// Get factorization sequence count. |
||
var countMultipleSequences = Fn.new { |n| |
var countMultipleSequences = Fn.new { |n| |
||
var f = Int.properDivisors(n) |
|||
f.removeAt(0) |
|||
var l = f.count |
|||
var d = List.filled(l, 0) |
|||
for (i in l-1..0) { |
|||
var k = f[i] |
|||
d[i] = 1 |
|||
var j = i + 1 |
|||
var x = 2 * k |
|||
while (x <= n) { |
|||
while (j < l && f[j] < x) j = j + 1 |
|||
if (j > l-1) break |
|||
if (f[j] == x) d[i] = d[i] + d[j] |
|||
x = x + k |
|||
} |
|||
} |
|||
return 1 + Nums.sum(d) |
|||
} |
} |
||
Line 357: | Line 373: | ||
System.print("\nOEIS A163272:") |
System.print("\nOEIS A163272:") |
||
var n = |
var n = 4 |
||
var fpns = [] |
var fpns = [0, 1] |
||
while (fpns.count < 7) { |
while (fpns.count < 7) { |
||
if ( |
if (countMultipleSequences.call(n) == n) fpns.add(n) |
||
n = n + |
n = n + 4 |
||
} |
} |
||
System.print(fpns)</syntaxhighlight> |
System.print(fpns)</syntaxhighlight> |
Revision as of 16:44, 7 October 2022
Consider the list of factors (divisors) of an integer, such as 12. The factors of 12 are [1, 2, 3, 4, 6, 12]. Consider all sorted sequences of the factors of n such that each succeeding number in such a sequnce is a multiple of its predecessor. So, for 6, we have the factors (divisors) [1, 2, 3, 6]. The 3 unique lists of sequential multiples starting with 1 and ending with 6 that can be derived from these factors are [1, 6], [1, 2, 6], and [1, 3, 6].
Another way to see these sequences is as an set of all the ordered factorizations of a number taken so that their product is that number (excluding 1 from the sequence). So, for 6, we would have [6], [2, 3], and [3, 2]. In this description of the sequences, we are looking at the numbers needed to multiply by, in order to generate the next element in the sequences previously listed in our first definition of the sequence type, as we described it in the preceding paragraph, above.
For example, for the factorization of 6, if the first type of sequence is [1, 6], this is generated by [6] since 1 * 6 = 6. Similarly, the first type of sequence [1, 2, 6] is generated by the second type of sequence [2, 3] because 1 * 2 = 2 and 2 * 3 = 6. Similarly, [1, 3, 6] is generated by [3, 2] because 1 * 3 = 3 and 3 * 2 = 6.
If we count the number of such sorted sequences of multiples, or ordered factorizations, and using that count find all integers `n` for which the count of such sequences equals `n`, we have re-created the sequence of the "factor-perfect" numbers (OEIS 163272).
By some convention, on its OEIS page, the factor-perfect number sequence starts with 0 rather than 1. As might be expected
with a sequence involving factorization and combinations, finding factor-perfect numbers becomes
more demanding on CPU time as the numbers become large.
- Task
- Show all 48 ordered sequences for each of the two methods for n = 48, which is the first non-trivial factor-perfect number.
- Write a program to calculate and show the first 7 numbers of the factor-perfect numbers.
- Stretch task
- Calculate and show more of the subsequent numbers in the sequence.
- see also
OEIS A163272 On the maximal order of numbers in the “factorisatio numerorum” problem Wikipedia: Enumerative Combinatorics
Julia
using Primes
""" Return the factors of n, including 1, n """
function factors(n::T)::Vector{T} where T <: Integer
sort(vec(map(prod, Iterators.product((p.^(0:m) for (p, m) in eachfactor(n))...))))
end
""" Uses the first definition and recursion to generate the sequences """
function more_multiples(to_seq, from_seq)
onemores = [[to_seq; i] for i in from_seq if i > to_seq[end] && i % to_seq[end] == 0]
isempty(onemores) && return Int[]
return append!(onemores, mapreduce(seq -> more_multiples(seq, from_seq), append!, onemores))
end
listing = sort!(push!(more_multiples([1], factors(48)[begin:end-1]), [1, 48]))
println("48 sequences using first definition:")
for (i, seq) in enumerate(listing)
print(rpad(seq, 20), i % 4 == 0 ? "\n" : "")
end
println("\n48 sequences using second definition:")
for (i, seq) in enumerate(listing)
seq[end] != 48 && push!(seq, 48)
seq2 = [seq[i] ÷ seq[i - 1] for i in 2:length(seq)]
print(rpad(seq2, 20), i % 4 == 0 ? "\n" : "")
end
""" Get factorization sequence count """
count_multiple_sequences(n) = length(more_multiples([1], factors(n)[begin:end-1])) + 1
println("\nOEIS A163272: ")
for n in 0:2_400_000
if n == 0 || count_multiple_sequences(n) == n
print(n, ", ")
end
end
- Output:
48 sequences using first definition: [1, 2] [1, 2, 4] [1, 2, 4, 8] [1, 2, 4, 8, 16] [1, 2, 4, 8, 24] [1, 2, 4, 12] [1, 2, 4, 12, 24] [1, 2, 4, 16] [1, 2, 4, 24] [1, 2, 6] [1, 2, 6, 12] [1, 2, 6, 12, 24] [1, 2, 6, 24] [1, 2, 8] [1, 2, 8, 16] [1, 2, 8, 24] [1, 2, 12] [1, 2, 12, 24] [1, 2, 16] [1, 2, 24] [1, 3] [1, 3, 6] [1, 3, 6, 12] [1, 3, 6, 12, 24] [1, 3, 6, 24] [1, 3, 12] [1, 3, 12, 24] [1, 3, 24] [1, 4] [1, 4, 8] [1, 4, 8, 16] [1, 4, 8, 24] [1, 4, 12] [1, 4, 12, 24] [1, 4, 16] [1, 4, 24] [1, 6] [1, 6, 12] [1, 6, 12, 24] [1, 6, 24] [1, 8] [1, 8, 16] [1, 8, 24] [1, 12] [1, 12, 24] [1, 16] [1, 24] [1, 48] 48 sequences using second definition: [2, 24] [2, 2, 12] [2, 2, 2, 6] [2, 2, 2, 2, 3] [2, 2, 2, 3, 2] [2, 2, 3, 4] [2, 2, 3, 2, 2] [2, 2, 4, 3] [2, 2, 6, 2] [2, 3, 8] [2, 3, 2, 4] [2, 3, 2, 2, 2] [2, 3, 4, 2] [2, 4, 6] [2, 4, 2, 3] [2, 4, 3, 2] [2, 6, 4] [2, 6, 2, 2] [2, 8, 3] [2, 12, 2] [3, 16] [3, 2, 8] [3, 2, 2, 4] [3, 2, 2, 2, 2] [3, 2, 4, 2] [3, 4, 4] [3, 4, 2, 2] [3, 8, 2] [4, 12] [4, 2, 6] [4, 2, 2, 3] [4, 2, 3, 2] [4, 3, 4] [4, 3, 2, 2] [4, 4, 3] [4, 6, 2] [6, 8] [6, 2, 4] [6, 2, 2, 2] [6, 4, 2] [8, 6] [8, 2, 3] [8, 3, 2] [12, 4] [12, 2, 2] [16, 3] [24, 2] [48] OEIS A163272: 0, 1, 48, 1280, 2496, 28672, 29808, 454656, 2342912,
Phix
You can run this online here.
with javascript_semantics function f(integer x) if x=1 then return {1} end if sequence res = {} for k=1 to x-1 do if remainder(x,k)=0 then for y in f(k) do res = append(res,y&x) end for end if end for res = sort(res) return res end function function m(sequence s, integer f, bool munge) sequence res = {} for x in s do x = deep_copy(x) if not munge then if length(x)>2 and x[$]=f then x = x[1..-2] end if res = append(res,x) else if x[$]!=f then x &= f end if for i=length(x) to 2 by -1 do x[i] /= x[i-1] end for res = append(res,x[2..$]) end if end for if not munge then res = sort(res) end if return res end function constant N = 48 sequence rN = f(N) function jbm(bool munge) rN = m(rN,N,munge) -- nb: true depends on false having happened return {length(rN),join_by(apply(rN,ppf),1,4," ",fmt:="%-16s")} end function ppOpt({pp_IntCh,false,pp_StrFmt,3}) printf(1,"%d sequences using first definition:\n%s\n",jbm(false)) printf(1,"%d sequences using second definition:\n%s\n",jbm(true)) function g(integer n) sequence f = factors(n) integer l = length(f) sequence d = repeat(0,l) for i=l to 1 by -1 do integer k = f[i] d[i] = 1 integer j = i+1 for x=2*k to n by k do while j<=l and f[j]<x do j+=1 end while if j>l then exit end if if f[j]=x then d[i] += d[j] end if end for end for return 1+sum(d) end function atom t = time(), t1 = t+1 integer n = 4 sequence res = {"0","1"} while length(res)<7 do if g(n)=n then res = append(res,sprintf("%d",n)) end if n += 4 if time()>t1 then progress("%d found, checking %d...\r",{length(res),n}) t1 = time()+1 end if end while progress("") printf(1,"Found %d: %s (%s)\n",{length(res),join(res," "),elapsed(time()-t)})
- Output:
Note that m(false) removes the final 48 where needed to match the Julia output, and m(true) puts it back on.
48 sequences using first definition: {1,2} {1,2,4} {1,2,4,8} {1,2,4,8,16} {1,2,4,8,24} {1,2,4,12} {1,2,4,12,24} {1,2,4,16} {1,2,4,24} {1,2,6} {1,2,6,12} {1,2,6,12,24} {1,2,6,24} {1,2,8} {1,2,8,16} {1,2,8,24} {1,2,12} {1,2,12,24} {1,2,16} {1,2,24} {1,3} {1,3,6} {1,3,6,12} {1,3,6,12,24} {1,3,6,24} {1,3,12} {1,3,12,24} {1,3,24} {1,4} {1,4,8} {1,4,8,16} {1,4,8,24} {1,4,12} {1,4,12,24} {1,4,16} {1,4,24} {1,6} {1,6,12} {1,6,12,24} {1,6,24} {1,8} {1,8,16} {1,8,24} {1,12} {1,12,24} {1,16} {1,24} {1,48} 48 sequences using second definition: {2,24} {2,2,12} {2,2,2,6} {2,2,2,2,3} {2,2,2,3,2} {2,2,3,4} {2,2,3,2,2} {2,2,4,3} {2,2,6,2} {2,3,8} {2,3,2,4} {2,3,2,2,2} {2,3,4,2} {2,4,6} {2,4,2,3} {2,4,3,2} {2,6,4} {2,6,2,2} {2,8,3} {2,12,2} {3,16} {3,2,8} {3,2,2,4} {3,2,2,2,2} {3,2,4,2} {3,4,4} {3,4,2,2} {3,8,2} {4,12} {4,2,6} {4,2,2,3} {4,2,3,2} {4,3,4} {4,3,2,2} {4,4,3} {4,6,2} {6,8} {6,2,4} {6,2,2,2} {6,4,2} {8,6} {8,2,3} {8,3,2} {12,4} {12,2,2} {16,3} {24,2} {48} Found 7: 0 1 48 1280 2496 28672 29808 (1.2s)
It takes about 6s under p2js, and 4½ minutes to get the next one on desktop/Phix, and I'm currently predicting Julia will take about the same time.
Python
''' Rosetta Code task Factor-perfect_numbers '''
from sympy import divisors
def more_multiples(to_seq, from_seq):
''' Uses the first definition and recursion to generate the sequences '''
onemores = [to_seq + [i]
for i in from_seq if i > to_seq[-1] and i % to_seq[-1] == 0]
if len(onemores) == 0:
return []
for i in range(len(onemores)):
for arr in more_multiples(onemores[i], from_seq):
onemores.append(arr)
return onemores
listing = [a + [48] for a in sorted(more_multiples([1], divisors(48)[1:-1]))] + [[1, 48]]
print('48 sequences using first definition:')
for j, seq in enumerate(listing):
print(f'{str(seq):22}', end='\n' if (j + 1) % 4 == 0 else '')
# Derive second definition's sequences
print('\n48 sequences using second definition:')
for k, seq in enumerate(listing):
seq2 = [seq[i] // seq[i - 1] for i in range(1, len(seq))]
print(f'{str(seq2):20}', end='\n' if (k + 1) % 4 == 0 else '')
def count_multiple_sequences(number):
''' Counts using the first definition, plus one extra for [1, n] '''
return len(more_multiples([1], divisors(number)[1:-1])) + 1
print("\nOEIS A163272: ", end='')
for num in range(500_000):
if num == 0 or count_multiple_sequences(num) == num:
print(num, end=', ')
- Output:
48 sequences using first definition: [1, 2, 48] [1, 2, 4, 48] [1, 2, 4, 8, 48] [1, 2, 4, 8, 16, 48] [1, 2, 4, 8, 24, 48] [1, 2, 4, 12, 48] [1, 2, 4, 12, 24, 48] [1, 2, 4, 16, 48] [1, 2, 4, 24, 48] [1, 2, 6, 48] [1, 2, 6, 12, 48] [1, 2, 6, 12, 24, 48] [1, 2, 6, 24, 48] [1, 2, 8, 48] [1, 2, 8, 16, 48] [1, 2, 8, 24, 48] [1, 2, 12, 48] [1, 2, 12, 24, 48] [1, 2, 16, 48] [1, 2, 24, 48] [1, 3, 48] [1, 3, 6, 48] [1, 3, 6, 12, 48] [1, 3, 6, 12, 24, 48] [1, 3, 6, 24, 48] [1, 3, 12, 48] [1, 3, 12, 24, 48] [1, 3, 24, 48] [1, 4, 48] [1, 4, 8, 48] [1, 4, 8, 16, 48] [1, 4, 8, 24, 48] [1, 4, 12, 48] [1, 4, 12, 24, 48] [1, 4, 16, 48] [1, 4, 24, 48] [1, 6, 48] [1, 6, 12, 48] [1, 6, 12, 24, 48] [1, 6, 24, 48] [1, 8, 48] [1, 8, 16, 48] [1, 8, 24, 48] [1, 12, 48] [1, 12, 24, 48] [1, 16, 48] [1, 24, 48] [1, 48] 48 sequences using second definition: [2, 24] [2, 2, 12] [2, 2, 2, 6] [2, 2, 2, 2, 3] [2, 2, 2, 3, 2] [2, 2, 3, 4] [2, 2, 3, 2, 2] [2, 2, 4, 3] [2, 2, 6, 2] [2, 3, 8] [2, 3, 2, 4] [2, 3, 2, 2, 2] [2, 3, 4, 2] [2, 4, 6] [2, 4, 2, 3] [2, 4, 3, 2] [2, 6, 4] [2, 6, 2, 2] [2, 8, 3] [2, 12, 2] [3, 16] [3, 2, 8] [3, 2, 2, 4] [3, 2, 2, 2, 2] [3, 2, 4, 2] [3, 4, 4] [3, 4, 2, 2] [3, 8, 2] [4, 12] [4, 2, 6] [4, 2, 2, 3] [4, 2, 3, 2] [4, 3, 4] [4, 3, 2, 2] [4, 4, 3] [4, 6, 2] [6, 8] [6, 2, 4] [6, 2, 2, 2] [6, 4, 2] [8, 6] [8, 2, 3] [8, 3, 2] [12, 4] [12, 2, 2] [16, 3] [24, 2] [48] OEIS A163272: 0, 1, 48, 1280, 2496, 28672, 29808, 454656,
Wren
However, using Phix's 'g' function for the last part which is considerably quicker than generating a bunch of lists and then counting them - overall time down from 76 to 7 seconds. Still not really worth attempting stretch goal though.
import "./math" for Int, Nums
import "./fmt" for Fmt
// Uses the first definition and recursion to generate the sequences.
var moreMultiples
moreMultiples = Fn.new { |toSeq, fromSeq|
var oneMores = []
for (i in fromSeq) {
if (i > toSeq[-1] && i%toSeq[-1] == 0) oneMores.add(toSeq + [i])
}
if (oneMores.isEmpty) return []
for (i in 0...oneMores.count) {
oneMores.addAll(moreMultiples.call(oneMores[i], fromSeq))
}
return oneMores
}
// Get factorization sequence count.
var countMultipleSequences = Fn.new { |n|
var f = Int.properDivisors(n)
f.removeAt(0)
var l = f.count
var d = List.filled(l, 0)
for (i in l-1..0) {
var k = f[i]
d[i] = 1
var j = i + 1
var x = 2 * k
while (x <= n) {
while (j < l && f[j] < x) j = j + 1
if (j > l-1) break
if (f[j] == x) d[i] = d[i] + d[j]
x = x + k
}
}
return 1 + Nums.sum(d)
}
var listing = moreMultiples.call([1], Int.properDivisors(48))
listing.add([1, 48])
listing.sort { |l1, l2|
var c1 = l1.count
var c2 = l2.count
for (i in 1...c1.min(c2)) {
if (l1[i] < l2[i]) return true
if (l1[i] > l2[i]) return false
}
if (c1 < c2) return true
return false
}
System.print("%(listing.count) sequences using first definition:")
Fmt.tprint("$-17n", listing, 4)
System.print("\n%(listing.count) sequences using second definition:")
var listing2 = []
for (i in 0...listing.count) {
var seq = listing[i]
if (seq[-1] != 48) seq.add(48)
var seq2 = (1...seq.count).map { |i| (seq[i]/seq[i-1]).floor }.toList
listing2.add(seq2)
}
Fmt.tprint("$-17n", listing2, 4)
System.print("\nOEIS A163272:")
var n = 4
var fpns = [0, 1]
while (fpns.count < 7) {
if (countMultipleSequences.call(n) == n) fpns.add(n)
n = n + 4
}
System.print(fpns)
- Output:
48 sequences using first definition: [1, 2] [1, 2, 4] [1, 2, 4, 8] [1, 2, 4, 8, 16] [1, 2, 4, 8, 24] [1, 2, 4, 12] [1, 2, 4, 12, 24] [1, 2, 4, 16] [1, 2, 4, 24] [1, 2, 6] [1, 2, 6, 12] [1, 2, 6, 12, 24] [1, 2, 6, 24] [1, 2, 8] [1, 2, 8, 16] [1, 2, 8, 24] [1, 2, 12] [1, 2, 12, 24] [1, 2, 16] [1, 2, 24] [1, 3] [1, 3, 6] [1, 3, 6, 12] [1, 3, 6, 12, 24] [1, 3, 6, 24] [1, 3, 12] [1, 3, 12, 24] [1, 3, 24] [1, 4] [1, 4, 8] [1, 4, 8, 16] [1, 4, 8, 24] [1, 4, 12] [1, 4, 12, 24] [1, 4, 16] [1, 4, 24] [1, 6] [1, 6, 12] [1, 6, 12, 24] [1, 6, 24] [1, 8] [1, 8, 16] [1, 8, 24] [1, 12] [1, 12, 24] [1, 16] [1, 24] [1, 48] 48 sequences using second definition: [2, 24] [2, 2, 12] [2, 2, 2, 6] [2, 2, 2, 2, 3] [2, 2, 2, 3, 2] [2, 2, 3, 4] [2, 2, 3, 2, 2] [2, 2, 4, 3] [2, 2, 6, 2] [2, 3, 8] [2, 3, 2, 4] [2, 3, 2, 2, 2] [2, 3, 4, 2] [2, 4, 6] [2, 4, 2, 3] [2, 4, 3, 2] [2, 6, 4] [2, 6, 2, 2] [2, 8, 3] [2, 12, 2] [3, 16] [3, 2, 8] [3, 2, 2, 4] [3, 2, 2, 2, 2] [3, 2, 4, 2] [3, 4, 4] [3, 4, 2, 2] [3, 8, 2] [4, 12] [4, 2, 6] [4, 2, 2, 3] [4, 2, 3, 2] [4, 3, 4] [4, 3, 2, 2] [4, 4, 3] [4, 6, 2] [6, 8] [6, 2, 4] [6, 2, 2, 2] [6, 4, 2] [8, 6] [8, 2, 3] [8, 3, 2] [12, 4] [12, 2, 2] [16, 3] [24, 2] [48] OEIS A163272: [0, 1, 48, 1280, 2496, 28672, 29808]