       Re: Re: Working with factors of triangular numbers.

• To: mathgroup at smc.vnet.net
• Subject: [mg80019] Re: [mg79864] Re: Working with factors of triangular numbers.
• From: DrMajorBob <drmajorbob at bigfoot.com>
• Date: Fri, 10 Aug 2007 06:44:09 -0400 (EDT)
• References: <200707030923.FAA17995@smc.vnet.net> <24165820.1186478071750.JavaMail.root@m35>

```Here's a significantly faster code (for n=12,13 anyway).

Clear[lldaux5, lld5, t5, factorCounts]
factorCounts = factorCounts = factorCounts = {1};
factorCounts = {2};
factorCounts = {3};
factorCounts[k_Integer] :=
If[EvenQ@k, factorCounts[(k + 2)/2, k - 1],
factorCounts[(k - 1)/2, k + 2]]
factorCounts[j_Integer, k_Integer] :=
Switch[GCD[j, k],
1,
Join[FactorInteger@j, FactorInteger@k],
3, Replace[
Join[FactorInteger[j/3],
FactorInteger[k/3]], {3, s_} :> {3, s + 2}, {1}]][[All, -1]] //
Sort
lldaux5[{}] := 0
lldaux5[x_List] /; MemberQ[x, 1] :=
Count[x, 1] + lldaux5@DeleteCases[x, 1]
lldaux5[Lpow_] :=
lldaux5[Lpow] =
With[{Mtup = Tuples[Range[0, #] & /@ Lpow]},
Total[Quiet@
LinearProgramming[ConstantArray[-1, Length[Mtup]],
Array[{0, 1} &, Length[Mtup]], Integers]] - 1]
lld5[n_Integer] := lldaux5@factorCounts@n
t5[n_] := # (# + 1)/2 &@NestWhile[# + 1 &, 2, lld5[#] < n &]

T2c /@ Range // Timing
t5 /@ Range // Timing

{0.422, {3, 15, 55, 253, 1081, 13861, 115921, 665281, 18280081,
75479041}}

{1.046, {3, 15, 55, 325, 1081, 18145, 226801, 665281, 18280081,
75479041}}

T2c // Timing
t5 // Timing

{1.625, 2080995841}

{3.063, 2080995841}

Timing differences aren't reliable either way for n<=11, but beyond that
it's a different story:

T2c // Timing
t5 // Timing

{27., 68302634401}

{17.187, 68302634401}

T2c // Timing
t5 // Timing

{150.172, 924972048001}

{64.032, 924972048001}

t5 // Timing

{516.546, 48318396825601}

The trick (already used by someone else, I think)

lldaux5[x_List] /; MemberQ[x, 1] :=
Count[x, 1] + lldaux5@DeleteCases[x, 1]

means that if the factor counts includes n ones, the LP problem used by
T2c has 2^n times as many variables as the one used by t5. Yet this
yielded only a 10-15% improvement, by itself.

A bigger gain was achieved in factorCounts, based on one of Carl Woll's
ideas in

http://forums.wolfram.com/mathgroup/archive/2007/Jul/msg00411.html

factorCounts[k_Integer] :=
If[EvenQ@k, factorCounts[(k + 2)/2, k - 1],
factorCounts[(k - 1)/2, k + 2]]

Several fine adjustments were required to get any improvement at all.
Defining factorCounts for 1,2,3,4, and 7 removed steps from the code for
other cases, and it was crucial to take full advantage of the fact that
the GCD of the two factors could only be 1 or 3.

Bobby

On Tue, 07 Aug 2007 00:31:22 -0500, sashap <pavlyk at gmail.com> wrote:

> 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:= 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 ma
ny
>> 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:= ffaux /@ Permutations[FactorInteger[chi][[All, 2]] ]
>
> Out= {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]],
>         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:= {Timing@T2, Timing@T2c}
>
> Out= {{3.829,2080995841},{3.921,2080995841}}
>
> In:= {Timing@T2, Timing@T2c}
>
> Out= {{26.532,68302634401},{22.797,68302634401}}
>
> In:= {Timing@T2, Timing@T2c}
>
> Out= {{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:= 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
>> 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, Tota l[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 factor=
s.
>>
>> >>Partitionable can be improved, but another bottleneck for larger ca=
ses
>> >>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 instea=
d
>> >>somehow combine FactorInteger[n+2] and FactorInteger[n-1].
>>
>> >>At any rate, it takes a bit longer, but with enough patience one ca=
n
>> >>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 fixe=
d
>> >>>it, I hope for the final time. Here is the new FF:
>>
>> >>>FFF[n_] := Module[{u = FactorInteger[n], s, k, partialQ, final=
Q, 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 =
>>
>> ...
>>
>
>
>
>

-- =

DrMajorBob at bigfoot.com

```

• Prev by Date: Re: Integration of Singular function
• Next by Date: Re: Integrate with PrincipalValue->True
• Previous by thread: Re: Working with factors of triangular numbers.
• Next by thread: Re: Re: Working with factors of triangular numbers.