Re: Working with factors of triangular numbers.
- To: mathgroup at smc.vnet.net
- Subject: [mg79784] Re: Working with factors of triangular numbers.
- From: sashap <pavlyk at gmail.com>
- Date: Sat, 4 Aug 2007 05:56:23 -0400 (EDT)
- References: <200707030923.FAA17995@smc.vnet.net>
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.
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 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.me... 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>
- Follow-Ups:
- Re: Re: Working with factors of triangular numbers.
- From: Carl Woll <carlw@wolfram.com>
- Re: Re: Working with factors of triangular numbers.