Re: Working with factors of triangular numbers.
- To: mathgroup at smc.vnet.net
- Subject: [mg79864] Re: Working with factors of triangular numbers.
- From: sashap <pavlyk at gmail.com>
- Date: Tue, 7 Aug 2007 01:31:22 -0400 (EDT)
- References: <200707030923.FAA17995@smc.vnet.net>
On Aug 5, 3:58 am, Carl Woll <ca... at wolfram.com> wrote:
> sashap wrote:
> >Carl's code is very good on numbers not involving
> >products of large powers of primes. Otherwise it
> >it is memory intensive and not as efficient:
>
> >In[109]:= Table[Timing[LargestPartition[{k, k}]], {k, 7, 15}]
>
> >Out[109]=
> >{{0.031,7},{0.141,7},{0.172,8},{0.187,9},{1.531,9},{1.875,10},
> >{12.016,10},{16.234,11},{141.391,11}}
>
> >One can use LinearProgramming to improve on the
> >situation.
>
> >Additional improvement comes from the following
> >observation. The longest factorization can be built
> >by taking the longest sequence of divisors of degree 1,
> >then of degree 2 and so on. The degree of the divisor
> >is defined as
>
> > deg[n_] := Total[ FactorInteger[n][[All,2]] ]
>
> >In other words, let n == p1^e1 *...* pn^en.
> >First take {p1, p2,.., pn}, then the longest
> >sequence of pairs with repetitions and so on.
>
> Nice work! However, I think there is a problem with your algorithm.
> There are many possible ways to a longest sequence of (distinct) pairs
> with repetitions. It is possible that some of these sequences will not
> produce the longest factorization. For instance, consider 2^10 * 3^3 *
> 5^2. After removing the degree 1 factors we are left with 2^9 * 3^2 * 5.
> Two possible length 3 sequences of pairs are:
>
> a) 2^2, 2*3, 2*5 leaving 2^5*3
>
> and
>
> b) 2^2, 2*3, 3*5 leaving 2^6
>
> Now, for case a) we can create two more unique factors, 2^3 and 2^2*3,
> while for case b) we can't create two more unique factors.
>
> So, I think one more thing is needed in your algorithm, a way of knowing
> which longest sequence of pairs to pick.
>
> Now, it turns out that the above example is handled correctly by your
> algorithm, that is, your algorithm by design or chance picks a longest
> sequence that leads to a longest factorization. I played around with a
> bunch of different examples, and came up with one example where your
> algorithm picks a longest sequence that doesn't lead to a longest
> factorization. This wasn't easy because my algorithm is soo slow (but I
> think, correct).
>
> In[594]:= counterexample = {2,3,4,5,6,7,8,10,12,14,15,35};
>
> In[595]:= Times @@ counterexample
>
> Out[595]= 35562240000
>
> In[596]:= Length[counterexample]
>
> Out[596]= 12
>
> In[597]:= LengthOfLongestDecomposition[Times @@ counterexample]
>
> Out[597]= 11
>
> My guess is that a simple heuristic probably exists which will allow you
> to pick a good longest sequence instead of a bad longest sequence.
> Perhaps something like choosing the longest sequence that keeps as many
> different types of prime factors available for the next stage? If so,
> then the biggest bottleneck will be just factoring the triangular numbers.
>
> Carl
>
Carl,
good points ! Indeed:
In[98]:= ffaux /@ Permutations[FactorInteger[chi][[All, 2]] ]
Out[98]= {12,12,11,12,12,12,12,12,12,11,11,11}
The heuristics you mention might be different sorting of exponents,
even though we can not guarantee it, that is
changing LengthOfLongestDecomposition to the following:
LengthOfLongestDecomposition[nn_Integer] :=
ffaux[Reverse@Sort@FactorInteger[nn][[All, 2]]]
Another approach is to write n == d1^s1 * ... * dk^sk and just use
LinearProgramming once. Here
s1, ..., sk are either 0 or 1, and d1, ..., dk are all distinct
divisors of the given integer n.
Clear[lldaux, lld];
lldaux[Lpow_] := (lldaux[Lpow] =
With[{Mtup = Tuples[ Range[0, #] & /@ Lpow]},
Total[
Quiet@LinearProgramming[ConstantArray[-1, Length[Mtup]],
Transpose@Mtup, Thread@{Lpow, 0},
Array[{0, 1} &, Length[Mtup]], Integers]] - 1])
lld[n_Integer] := lldaux[FactorInteger[n][[All, 2]] // Sort]
T2c[n_] := # (# + 1)/2 &@
NestWhile[# + 1 &, 2, lld[# (# + 1)/2 - 1] < n &]
The performance of T2c is competitive to Carl's code for small values
of the argument
and has an edge over it for large ones:
In[168]:= {Timing@T2[11], Timing@T2c[11]}
Out[168]= {{3.829,2080995841},{3.921,2080995841}}
In[169]:= {Timing@T2[12], Timing@T2c[12]}
Out[169]= {{26.532,68302634401},{22.797,68302634401}}
In[170]:= {Timing@T2[13], Timing@T2c[13]}
Out[170]= {{144.406,924972048001},{123.875,924972048001}}
Oleksandr Pavlyk and Maxim Rytin
> >In order to choose the longest sequence of factors of a given degree,
> >we set this up as an integer linear programming problem. We illustrate
> >this for degree 2. Let p1*p1, p1*p2, and so on be candidate factors of
> >degree 2. Then
>
> > (p1*p1)^s1 (p1*p2)^s2 ...
>
> >must divide the original number. Additionally, each sk is either 0 or
> >1.
> >This gives us the constraints for variables s1, s2, ... and
> >the objective function being maximized is s1+s2+....
>
> >The code implementing those ideas:
>
> >CombinationsWithRepetitions =
> > Compile[{{n, _Integer}, {k, _Integer}},
> > Module[{ans, i = 1, j = 1, l = 1},
> > ans = Array[0 &, {Binomial[n + k - 1, k], k}];
> > While[True, While[j <= k, ans[[i, j++]] = l];
> > If[i == Length@ans, Break[]];
> > While[ans[[i, --j]] == n,];
> > l = ans[[i++, j]] + 1;
> > ans[[i]] = ans[[i - 1]];];
> > ans]];
>
> >Clear[LengthOfLongestDecomposition, ffaux]
> >LengthOfLongestDecomposition[nn_] :=
> > ffaux[Sort@FactorInteger[nn][[All, 2]]]
> >ffaux[$Lpow_] :=
> > ffaux[$Lpow] =
> > Module[{Lpow = $Lpow, Ldiv, n, m, k = 1, Lind, ans = 0},
> > n = Length@Lpow;
> > While[Total@Lpow >= k,
> > Ldiv = CombinationsWithRepetitions[n, k++];
> > m = Length@Ldiv;
> > Lind =
> > Quiet@LinearProgramming[ConstantArray[-1, m],
> > Total[1 - Unitize[Ldiv - #], {2}] & /@ Range@n,
> > Thread@{Lpow, -1}, Array[{0, 1} &, m], Integers];
> > ans += Total@Lind;
> > Lpow -= BinCounts[Flatten@Pick[Ldiv, Lind, 1], {Range[n + 1]}]];
> > ans]
>
> >Compare the timing given earlier with the following
>
> >In[132]:= Table[Timing[ffaux[{k, k}]], {k, 7, 15}]
>
> >Out[132]= \
> >{{0.016,7},{1.78746*10^-14,7},{1.78746*10^-14,8},{0.015,9},{3.15303*\
> >10^-14,9},{0.016,10},{0.,10},{0.015,11},{0.,11}}
>
> >Clearly, the use of heavy machinery of LinearProgramming comes at the
> >cost of
> >noticeable overhead which is why the performance of this algorithm is
> >comparable
> >to that of Carl's code when measured over a sequence of triangular
> >numbers because
> >problematic numbers are rare.
>
> >T2b[n_] := # (# + 1)/2 &@
> > NestWhile[# + 1 &, 1,
> > LengthOfLongestDecomposition[# (# + 1)/2 - 1] < n &]
>
> >In[66]:= {Timing[T2[11]], Timing[T2b[11]]}
>
> >Out[66]= {{3.172,2080995841},{2.328,2080995841}}
>
> >In[67]:= {Timing[T2[12]], Timing[T2b[12]]}
>
> >Out[67]= {{26.312,68302634401},{18.235,68302634401}}
>
> >In[68]:= {Timing[T2[13]], Timing[T2b[13]]}
>
> >Out[68]= {{142.765,924972048001},{111.422,924972048001}}
>
> >Oleksandr Pavlyk and Maxim Rytin
>
> >On Jul 10, 5:37 am, Carl Woll <ca... at wolfram.com> 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 Partition=
able.
>
> >>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[li=
st]]];
> >> 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, s=
pace},
> >>> 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 =
>
> ...
>
> read more =BB