Re: Re: Working with factors of triangular numbers.
- To: mathgroup at smc.vnet.net
- Subject: [mg78814] Re: [mg78766] Re: [mg78490] Working with factors of triangular numbers.
- From: Andrzej Kozlowski <akoz at mimuw.edu.pl>
- Date: Tue, 10 Jul 2007 06:31:47 -0400 (EDT)
- References: <200707030923.FAA17995@smc.vnet.net> <6D2F6E70-2462-4B69-A617-4DEC322D69BF@mimuw.edu.pl> <93E8C621-E93F-4761-A0FF-205BF84249CD@mimuw.edu.pl> <a851af150707070640g1f6a22e9ue072d4eba5c9262a@mail.gmail.com> <6F7B3048-4C96-48BC-AADF-2A62837DE6D3@mimuw.edu.pl> <a851af150707072125o69edf471o4b5be8932a924a32@mail.gmail.com> <EC981B8F-E91E-4CB7-A4C2-030EE8C5F319@mimuw.edu.pl> <B9DF7620-D7E2-49D7-B44E-619C8EB636BF@mimuw.edu.pl> <a851af150707080555r966885cl6a03f44f637ac3ef@mail.gmail.com> <a851af150707080608u4040f2a8y5368c781342d05c1@mail.gmail.com> <a851af150707080617s436ad717uc39e1e68cb974315@mail.gmail.com> <200707090533.BAA08008@smc.vnet.net> <4692812C.6050301@wolfram.com>
I am very happy to see a recognized Grand Master of The Art of
Lightning Fast Coding join this thread!
Your code is , as usual, spectaculary fast although I am not quite
convinced that your "algorithm" is actually faster than backtracking
although it is certainly faster than my implementation of it, based
on Combinatorica's Backtrack function. Usually (there have been a
quite many examples of that on this list in the past) the advantage
of "custom" backtracking code is that it is memory efficient and can
be compiled. In fact, I think, compiling is essential for
backtracking code to be really useful and it can't be achieved with
Combinatorica's Backtrack function. If this were done by one of the
other recognized Grand Masters (e.g. Fred Simons) then I think it
could provide a real challenge to your code.
But this is not the kind of challenge I wish to take up.
However, thinking again about this problem I noticed a horrible
elementary programming howler that I made in my implementation of
backtracking for this problem so I really feel oblidged to correct
it. What happened is that I failed to take account of commutativity
of multiplication (!!) so my code is actually computing all possible
factorizations that differ only in the order of factors. This, of
course, produces horrible performance. So, here is corrected version
of my backtracking code, which no longer has this fault.
<< Combinatorica`
FF[n_] := Module[{u = FactorInteger[n], s, k, partialQ, finalQ, space},
s = u[[All, 2]]; k = Length[u];
partialQ[l_List] /; Length[l] == 1 := True; partialQ[l_List] :=
And @@
Flatten[{OrderedQ[Take[l, -2]],
Last[l] == Array[0 &, k] || l[[-1]] != l[[-2]],
Thread[Total[l] <= s - 1]}];
finalQ[l_List] := Total[l] == s - 1;
space = Table[DeleteCases[Tuples[(Range[0, #1] & ) /@ (s - 1)],
Alternatives @@ IdentityMatrix[k], Infinity], {k}];
k +
Max[0, Length /@ DeleteCases[Backtrack[space, partialQ, finalQ,
All],
Array[0 & , k], Infinity]]]
T[n_] := Block[{k = 1, $Messages},
While[k++; FF[k*((k + 1)/2) - 1] < n, Null]; k*((k + 1)/2)]
with that we get:
T[6] // Timing
{0.851239, 13861}
This is still much slower than your code but a lot faster than my
previous implementation.
It is also useful to have a fucntion that actually computes the
factorizations. One can make one just by changing the last line in
the code for FF:
GG[n_] := Module[{u = FactorInteger[n], s, k, partialQ, finalQ, space},
s = u[[All, 2]]; k = Length[u];
partialQ[l_List] /; Length[l] == 1 := True; partialQ[l_List] :=
And @@
Flatten[{OrderedQ[Take[l, -2]],
Last[l] == Array[0 &, k] || l[[-1]] != l[[-2]],
Thread[Total[l] <= s - 1]}];
finalQ[l_List] := Total[l] == s - 1;
space = Table[DeleteCases[Tuples[(Range[0, #1] & ) /@ (s - 1)],
Alternatives @@ IdentityMatrix[k], Infinity], {k}];
DeleteCases[Backtrack[space, partialQ, finalQ, All], Array[0 &,
k],
Infinity]]
For example:
GG[665280]
{{{5, 2, 0, 0, 0}}, {{0, 2, 0, 0, 0}, {5, 0, 0, 0, 0}},
{{1, 1, 0, 0, 0}, {4, 1, 0, 0, 0}}, {{1, 2, 0, 0, 0}, {4, 0, 0,
0, 0}},
{{2, 0, 0, 0, 0}, {3, 2, 0, 0, 0}}, {{2, 1, 0, 0, 0}, {3, 1, 0,
0, 0}},
{{2, 2, 0, 0, 0}, {3, 0, 0, 0, 0}}, {{0, 2, 0, 0, 0}, {2, 0, 0,
0, 0},
{3, 0, 0, 0, 0}}, {{1, 1, 0, 0, 0}, {2, 0, 0, 0, 0}, {2, 1, 0,
0, 0}}}
Let me quickly explain what this means. We have:
ff = FactorInteger[665280]
{{2, 6}, {3, 3}, {5, 1}, {7, 1}, {11, 1}}
The prime factors are:
ff[[All, 1]]
{2, 3, 5, 7, 11}
In the above output of GG[665280], the list {{5, 2, 0, 0, 0}}
corresponds to the factor
Times @@ ({2, 3, 5, 7, 11}^{5, 2, 0, 0, 0}) = 288.
In a maximal factorization there are always going to be present
all the prime factors 2*3*5*7*11 and then some other factors and only
these other factors are found by GG so the above gives the
factorization of length 5+1 = 6 (5 primes + one other factor)
665280 == 2*3*5*7*11*288
True
The list
{{3, 0, 0, 0, 0}}, {{1, 1, 0, 0, 0}, {2, 0, 0, 0, 0}, {2, 1, 0, 0,
0}} gives us a factorization of maximal possible length, with the
prime factors {2, 3, 5, 7, 11} and the extra factors:
Times @@ ({2, 3, 5, 7, 11}^#) & /@ {{1, 1, 0, 0, 0}, {2, 0, 0, 0,
0}, {2, 1,
0, 0, 0}}
{6, 4, 12}
so the factorization has length 5+ 3:
665280 == 2*3*5*7*11*6*4*12
True
Al this is, of course, much slower than your code, so I am writing it
only for the sake of "education", although, as I wrote above, I still
think that a compiled code using the same algorithm could be quite
competitive - I have seen in the past compiled backtrackking code two
orders of magnitude faster than a version using the Backtrack function.
Andrzej Kozlowski
On 10 Jul 2007, at 03:40, Carl Woll wrote:
> Andrzej and Diana,
>
> Here is a faster algorithm. First, an outline:
>
> For a given number, the easy first step in writing it as a
> factorization of the most distinct factors is to use each prime as
> a factor. Then, the hard part is to take the remaining factors and
> partition them in a way to get the most distinct factors. However,
> note that this step only depends on the counts of each prime
> factor, and not on the values of the primes themselves. Hence, we
> can memoize the number of distinct factors possible for a set of
> counts. I do this as follows:
>
> Let f be the list of remaining factors, of length l.
>
> 1. Determine possible integer partitions of l, with the smallest
> partition allowed being 2. We also need to order these partitions
> by length, from largest to smallest:
>
> partitions = IntegerPartitions[l, All, Range[2, l]];
> partitions = Reverse@partitions[[ Ordering[Length/@partitions] ]];
>
> 2. Now, we need to test each partition to see if it's possible to
> fill in the partitions with the factors such that each partition is
> unique. Once we find such a partition, we are done, and we know how
> many distinct factors that number can be written as. This step I do
> by recursion, i.e., take the first member of a partition and fill
> it in with one of the possible subsets of factors, and then repeat
> with the remaining members of a partition. I do this with the
> function Partitionable.
>
> Here is the code:
>
> LargestPartition[list : {1 ..}] := Length[list]
> LargestPartition[list_List] := Length[list] + LargestTwoPartition
> [Reverse@Sort@list - 1 /. {a__, 0 ..} -> {a}]
>
> Clear[LargestTwoPartition]
> LargestTwoPartition[list_] := LargestTwoPartition[list] = Module
> [{set, partitions, res},
> set = CountToSet[list];
> partitions = IntegerPartitions[Total[list], All, Range[2, Total
> [list]]];
> partitions = Reverse@partitions[[Ordering[Length /@ partitions]]];
> res = Cases[partitions, x_ /; Partitionable[{}, set, x], 1, 1];
> If[res === {},
> 0,
> Length@First@res
> ]
> ]
>
> Partitionable[used_, unused_, part_] := Module[{first, rest},
> first = Complement[DistinctSubsets[unused, First@part], used];
> If[first === {}, Return[False]];
> rest = CountToSet /@ Transpose[
> SetToCount[unused, Max[unused]] -
> Transpose[SetToCount[#, Max[unused]] & /@ first]
> ];
> Block[{Partitionable},
> Or @@ MapThread[
> Partitionable[Join[used, {#1}], #2, Rest[part]] &,
> {first, rest}
> ]
> ]
> ]
>
> Partitionable[used_, rest_, {last_}] := ! MemberQ[used, rest]
>
> CountToSet[list_] := Flatten@MapIndexed[Table[#2, {#1}] &, list]
> SetToCount[list_, max_] := BinCounts[list, {1, max + 1}]
>
> DistinctSubsets[list_, len_] := Union[Subsets[list, {len}]]
>
> T2[n_] := Module[{k=1},
> While[k++; LargestPartition[FactorInteger[k*((k + 1)/2) - 1][[All,
> 2]]] < n];
> k*((k + 1)/2)
> ]
>
> Then:
>
> In[11]:= T2[6] // Timing
> Out[11]= {0.047,13861}
>
> In[12]:= T2[7] // Timing
> Out[12]= {0.031,115921}
>
> In[13]:= T2[8] // Timing
> Out[13]= {0.11,665281}
>
> In[14]:= T2[9] // Timing
> Out[14]= {0.625,18280081}
>
> In[15]:= T2[10] // Timing
> Out[15]= {1.157,75479041}
>
> In[16]:= T2[11] // Timing
> Out[16]= {5.875,2080995841}
>
> In[17]:= T2[12] // Timing
> Out[17]= {48.703,68302634401}
>
> For Diana, note that
>
> In[20]:= 481 482/2 - 1 == 115921 - 1 == 2 3 4 5 6 7 23
> Out[20]= True
>
> so 115921-1 can indeed be written as a product of 7 distinct factors.
>
> Partitionable can be improved, but another bottleneck for larger
> cases is evaluating FactorInteger. One possibility for improving
> FactorInteger speed is to note that
>
> In[21]:= n (n + 1)/2 - 1 == (n + 2) (n - 1)/2 // Expand
>
> Out[21]= True
>
> So, rather than applying FactorInteger to n(n+1)/2-1 you can
> instead somehow combine FactorInteger[n+2] and FactorInteger[n-1].
>
> At any rate, it takes a bit longer, but with enough patience one
> can also find:
>
> In[53]:= 1360126 1360127/2 - 1 == 924972048001 - 1 == 2 3 4 5 6 7
> 10 11 12 13 15 23 31
> Out[53]= True
>
> and
>
> In[56]:= 9830401 9830402/2 - 1 == 48318396825601 - 1 == 2 3 4 5 6
> 8 9 10 11 12 17 22 32 59
> Out[56]= True
>
> Carl Woll
> Wolfram Research
>
> Andrzej Kozlowski wrote:
>
>> Well, I stayed up longer than I wanted and I think I have now
>> fixed it, I hope for the final time. Here is the new FF:
>>
>> FFF[n_] := Module[{u = FactorInteger[n], s, k, partialQ, finalQ,
>> space},
>> s = u[[All,2]]; k = Length[u]; partialQ[l_List] :=
>> And @@ Flatten[{Last[l] == Array[0 & , k] ||
>> !MemberQ[Most[l], Last[l]], Thread[Total[l] <= s - 1]}];
>> finalQ[l_List] := And @@ Flatten[{Last[l] == Array[0 & , k] ||
>> !MemberQ[Most[l], Last[l]], Total[l] == s - 1}];
>> space = Table[DeleteCases[Tuples[(Range[0, #1] & ) /@ (s - 1)],
>> Alternatives @@ IdentityMatrix[k], Infinity], {k}];
>> k + Max[0, Length /@ DeleteCases[Backtrack[space, partialQ,
>> finalQ, All],
>> Array[0 & , k], Infinity]]]
>>
>> T[n_] := Block[{k = 1, $Messages},
>> While[k++; FF[k*((k + 1)/2) - 1] < n, Null]; k*((k + 1)/2)]
>>
>> This gives the right answers for the first 8 values and in
>> particular:
>>
>> T[8]
>> 665281
>>
>> I am pretty sure it now works fine! ;-))
>>
>> Andrzej
>>
>>
>> On 8 Jul 2007, at 22:17, Diana Mecum wrote:
>>
>>
>>> Andrzej,
>>>
>>> Please disregard my last e-mail.
>>>
>>> The eighth term of the sequence should be 665281
>>>
>>> 665280 = 11*7*5*3*9*2*4*8
>>>
>>> Diana
>>>
>>> On 7/8/07, Diana Mecum <diana.mecum at gmail.com> wrote: Andrzej,
>>>
>>> I appreciate all of your work. I have stopped working on this
>>> problem,
>>>
>>> The best first 8 terms I have at this point are:
>>>
>>> 3,15,55,253,1081,13861,138601,665281
>>>
>>> The first 8 terms I get with your algorithm are:
>>>
>>> 3,15,55,253,1081,13861,115921,1413721
>>>
>>> 115921 = 13 * 37 * 241, which does not fit the rule.
>>>
>>> I would be interested in any further information you would have,
>>> but also would understand your not wanting to take further time
>>> with this.
>>>
>>> Thanks again,
>>>
>>> Diana M.
>>>
> <snip>
- References:
- Working with factors of triangular numbers.
- From: Diana <diana.mecum@gmail.com>
- Re: Working with factors of triangular numbers.
- From: Andrzej Kozlowski <akoz@mimuw.edu.pl>
- Working with factors of triangular numbers.