       Re: Re: Working with factors of triangular numbers.

• To: mathgroup at smc.vnet.net
• Subject: [mg79811] Re: [mg79784] Re: Working with factors of triangular numbers.
• From: Carl Woll <carlw at wolfram.com>
• Date: Sun, 5 Aug 2007 04:53:48 -0400 (EDT)
• References: <200707030923.FAA17995@smc.vnet.net> <200708040956.FAA05882@smc.vnet.net>

```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:= Table[Timing[LargestPartition[{k, k}]], {k, 7, 15}]
>
>Out=
>{{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:= counterexample =  {2,3,4,5,6,7,8,10,12,14,15,35};

In:= Times @@ counterexample

Out= 35562240000

In:= Length[counterexample]

Out= 12

In:= LengthOfLongestDecomposition[Times @@ counterexample]

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

>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:= Table[Timing[ffaux[{k, k}]], {k, 7, 15}]
>
>Out= \
>{{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:= {Timing[T2], Timing[T2b]}
>
>Out= {{3.172,2080995841},{2.328,2080995841}}
>
>In:= {Timing[T2], Timing[T2b]}
>
>Out= {{26.312,68302634401},{18.235,68302634401}}
>
>In:= {Timing[T2], Timing[T2b]}
>
>Out= {{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},
>>      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:= T2 // Timing
>>Out= {0.047,13861}
>>
>>In:= T2 // Timing
>>Out= {0.031,115921}
>>
>>In:= T2 // Timing
>>Out= {0.11,665281}
>>
>>In:= T2 // Timing
>>Out= {0.625,18280081}
>>
>>In:= T2 // Timing
>>Out= {1.157,75479041}
>>
>>In:= T2 // Timing
>>Out= {5.875,2080995841}
>>
>>In:= T2 // Timing
>>Out= {48.703,68302634401}
>>
>>For Diana, note that
>>
>>In:= 481 482/2 - 1 == 115921 - 1 == 2 3 4 5 6 7 23
>>Out= 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:= n (n + 1)/2 - 1 == (n + 2) (n - 1)/2 // Expand
>>
>>Out= 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:= 1360126 1360127/2 - 1 == 924972048001 - 1 ==  2 3 4 5 6 7 10 11
>>12 13 15 23 31
>>Out= True
>>
>>and
>>
>>In:= 9830401 9830402/2 - 1 == 48318396825601 - 1 ==  2 3 4 5 6 8 9
>>10 11 12 17 22 32 59
>>Out= 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
>>>665281
>>>
>>>
>>>I am pretty sure it now works fine!  ;-))
>>>
>>>
>>>Andrzej
>>>
>>>
>>>On 8 Jul 2007, at 22:17, Diana Mecum wrote:
>>>
>>>
>>>>Andrzej,
>>>>
>>>>
>>>>
>>>>
>>>>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>
>>
>>
>
>
>
>
>
>

```

• Prev by Date: Re: Documentation Center (v6): do-it-yourself Mathematica Book
• Next by Date: Re: Documentation Center (v6): do-it-yourself Mathematica
• Previous by thread: Re: Working with factors of triangular numbers.
• Next by thread: Re: Working with factors of triangular numbers.