Question: In or out?

Every time a new set is defined, the question arises about whether or not an object is in or out, i.e. whether it is an element of the set or not. Just to answer this question I provided the XGCQ function, that works like the PrimeQ function. It gets the XGC parameters and the object you want to test for membership, and returns true or false whether the object is in or out.

XGCQ[ c, m, n ]
XGCQ[ c_Integer, m_Integer, n_Integer, option_Rule:(Method -> One) ] :=
  Module[ {opt},
    opt = Method /. option;
    Switch[ opt,
      One, XGCfastQ[ c, m, n ],
      All, XGCslowQ[ c, m, n ],
      _, "-Invalid Option-"
    ]
  ] /; Mod[ n, m ] == 0

It yields True if n belongs to XGC(c, m), and it yields False otherwise.
XGCQ[ c_Integer, m_Integer, n_Integer, Method -> All ] uses XGC (slow and intuitive), i.e. "All" the XGC[ c, m, n - (m - 1) (c + m) ] is computed and then it is verified whether n is a member or not.
XGCQ[ c_Integer, m_Integer, n_Integer, Method -> One ], the default, uses FastPartition (fast and complex), i.e. Euclid[ c, m, n - (m - 1) (c + m) ] is searched for only "One" partition of n using m of its terms.

As you can see, XGCQ is only a wrapper for the two functions XGCslowQ and XGCfastQ, which you can select by means of the Method parameter. By default, XGCQ uses XGCfastQ to get a fast answer, but, if you want and have time to waste, you can specify Method -> All, to get an answer that is much more "secure" because you can easily understand its algorithm.

 

A naive method

XGCslowQ uses a naive method to accomplish its task. It calculates all the XGC(c, m) set up to a sufficiently large element T and then checks the set for the membership of the object n. The value of T is n - (m - 1) (c + m) because n = x1 + x2 + ... + xm and Min[ xi ] = Min[ XGC[ c, m, infinity ] ] = c + m. So xk <= T = n - Sum[ Min[ xi ], {i, 1, m - 1} ] = n - Sum[ c + m, {i, 1, m - 1} ] = n - (m - 1) (c + m).
The problem with such a method is the amount of unnecessary calculations it involves, which directly lead to a big waste of time. On the contrary, it is so a simple method it will be liked by a lot of people.

XGCslowQ[ c, m, n ]
XGCslowQ[ c_Integer, m_Integer, n_Integer ] :=
  MemberQ[ XGC[ c, m, n - (m - 1) (c + m) ], n ]

 

A fast method

XGCfastQ uses a complicated method to get an answer as soon as possible. It looks for a partition of the object n, using exactly m elements of the Euclid(c, m) set calculated up to n - (m - 1) (c + m). If such a partition, whatever it be, is found then XGCfastQ returns True, otherwise returns False.

XGCfastQ[ c, m, n ]
XGCfastQ[ c_Integer, m_Integer, n_Integer ] :=
  FastPartition[ n, m, Euclid[ c, m, n - (m - 1) (c + m) ] ] != {}

The complicated part of this method is the search for a partition of n, if any exists.
In general, we have to perform an exhaustive search in a restriction of the search space, i.e. the space of all the m-permutations of the terms of the Euclid(c, m) set. Firstly, we can cut off all the permutations that are variations of the same m-combination (with repetitions) of elements of Euclid(c, m). For example, there is no need to test the permutations 132, 213, 231, 312, and 321 if we already failed on 123. To understand that this is still an enormous amount of combinations let's test the following function.

CombinationsCount[ n, k ]
CombinationsCount[ n_Integer, k_Integer ] := (n - 1 + k)! / ((n - 1)! k!)

It gives the length of Combinations[ n, k ].

 

Combinations

The set of k-combinations of n objects can be traversed downward and upward in an ordered way. Such a property is extremely useful to program an exhaustive search, where you have to test a combination, and, depending on the result, you need or not the adjacent combination. The following two functions (simple enumerators) are used to get the previous and the next combination of a given one.

PrevCombination[ comb ]
pPrevCombination[ comb_List ] :=
  Module[ {newComb = comb, i, j, combSize = Length[ comb ]},
    If[ newComb[[ 1 ]] > 1,
      newComb[[ 1 ]]--;
      Return[ newComb ]
    , If[ combSize > 1,
        For[ i = 2, (i <= combSize) && (newComb[[ i ]] == 1), i++ ];
        If[ i <= combSize,
          newComb[[ i ]]--;
          For[ j = 1, j < i, j++, 
            newComb[[ j ]] = newComb[[ i ]] 
          ];
          Return[ newComb ]
        , Return[ {} ]
        ]
      , Return[ {} ]
      ]
    ]
  ]

PrevCombination[ comb_List ] := pPrevCombination[ comb ]  /; OrderedQ[ comb ]

It orderly enumerates downward all the k-conbinations (k = Length[ comb ]) of some elements with repetitions. 1 is the minimum value valid for a term of comb.

NextCombination[comb, max]
pNextCombination[comb_List, max_Integer] :=
  Module[ {newComb = comb, i, j, combSize = Length[ comb ]},
    If[ (combSize > 1) && (newComb[[ 1 ]] < newComb[[ 2 ]]) || (combSize == 1) && (newComb[[ 1 ]] < max), 
      newComb[[ 1 ]]++;
      Return[ newComb ]
    , If[ combSize > 1, 
        For[ i = 2, (i < combSize) && (newComb[[ i+1 ]] == newComb[[ i ]]), i++ ];
        If[ newComb[[ i ]] < max, 
          newComb[[ i ]]++;
          For[ j = 1, j < i, j++, 
            newComb[[ j ]] = 1 
          ];
          Return[ newComb ]
        , Return[ {} ]
        ]
      , Return[ {} ]
      ]
    ]
  ]

NextCombination[comb_List, max_Integer] := pNextCombination[comb, max]  /; OrderedQ[ comb ] && IsSetLE[ comb, max ]

It orderly enumerates upward all the k-conbinations (k = Length[ comb ]) of max elements with repetitions. max is the maximum value valid for a term of comb.

Just to see what is a set of all the k-combinations of n elements, let's test the following function. As you will note, the set is sorted in ascending order. Tip: to see it, you have to revert each combination.

Combinations[ n, k ]
Combinations[ n_Integer, k_Integer ]:= 
  Module[ {comb = Table[ 1, {k} ], set = {}}, 
    While[ comb != {}, 
      AppendTo[ set, comb ];
      comb = pNextCombination[ comb, n ]
    ];
    set
  ]

It gives the set of the k-combinations of n items WITH REPETITIONS.

 

FastPartition

FastPartition is a twofold function. The first component uses a greedy method to get a partition of n in m parts of the set in a bit of time. Not only that method is fast, but it has a very high hit rate, too. Unfortunately, sometimes it fails, so takes the field the second component. It uses an exhaustive method that is highly optimized in dynamically reducing the search space as much as possible.

FastPartition[ n, m, set ]
FastPartition[ n_Integer, m_Integer, set_List ] :=
  Module[ {extra, extraIndex, greedyComb, prevComb, nextComb, minCombSet, maxCombSet},

(* greedy method *)
    greedyComb = FastLoading[ n, m, set ];
    If[ n > m Last[ set ],
      Return[ {} ]
    ];
    extra = n - Apply[ Plus, set[[ greedyComb ]] ];
    If[ extra == 0,
      Return[ set[[ greedyComb ]] ]
    ];

(* exhaustive method *)
(* enumerate all the m-combinations of n items WITH REPETITIONS, 
 * going from greedy result down to bottom and up to top
 *)
    minCombSet = {Table[ 1, {m - 1} ]};
    prevComb = pPrevCombinationGTmin[ Rest[ greedyComb ], minCombSet ];

    maxCombSet = {Table[ Length[ set ], {m - 1} ]};
    nextComb = pNextCombinationLTmax[ Rest[ greedyComb ], Length[ set ], maxCombSet ];

    While[ prevComb != {} || nextComb != {},

      If[ prevComb != {},
        extra = n - Apply[ Plus, set[[ prevComb ]] ];
        extraIndex = BinSearchEQ[ extra, set ];
        If[ extraIndex != 0,
          Return[ set[[ Flatten[ {extraIndex, prevComb} ] ]] ]
        ];
        If[ extra > Last[ set ], (* then prevComb is a new limit *)
          PrependTo[ minCombSet, prevComb ]
        ];
        prevComb = pPrevCombinationGTmin[ prevComb, minCombSet ]
      ];

      If[ nextComb != {},
        extra = n - Apply[ Plus, set[[ nextComb ]] ];
        extraIndex = BinSearchEQ[ extra, set ];
        If[ extraIndex != 0,
          Return[ set[[ Flatten[ {extraIndex, nextComb} ] ]] ]
        ];
        If[ extra < First[ set ], (* then nextComb is a new limit *)
          PrependTo[ maxCombSet, nextComb ]
        ];
        nextComb = pNextCombinationLTmax[ nextComb, Length[ set ], maxCombSet ]
      ];

    ];
    Return[ {} ]
  ]

It gives a m-partition of n with numbers of the set.

 

FastPartition: Exhaustive Search

The search space is greatly reduced (nearly 40%, 30%, and 20% in the examples above) by the use of the "combinations versus permutations" concept, but that reduction is still not sufficient to find an exhaustive answer in a human lifetime. For example, if you could test one billion combinations per second (10^9/s) you will need more than 3239 years to be sure that a partition of n by means of m17with100terms does not exist.

Among all the m-combinations there is a huge lot of them that are not partitions of n. So it would mean a great save of time if we could not explore all those combinations that a priori are not partitions of n.

For example, if we take the Combinations[ 5, 3 ] set, then there is no need to test the 334, 244, 344, 444, 235, 335, 245, 345, 445, 255, 355, 455, and 555 combinations if we already under failed on 234.

Under/Over Failing
When we search for a partition n = term1 + term2 + ... termm - 1 + termm, where all the terms belong to a given set, we can search for a combination (term1, term2, ... termm - 1) of the terms of the set such that n - (term1 + term2 ... + termm - 1) = restm.
Now, if the restm is less than the minimum (greater than the maximum) of the set, then we have under (over) failed on that combination.

Back to the example, by means of the under/over failing concept (and the fact that the Euclid set is ordered), you can easily understand why you can save the test of all those combinations. Any of them is such that always the 3 terms are greater than or equal to the terms 2, 3, and 4, respectively. So, we can say a priori that we will under fail on any of those combinations.

During the search process, in general we under/over fail on a lot of different combinations. So it is a good idea to maintain a pool of max/min combinations to use as thresholds for the search process. To traverse the set of combinations using those thresholds, I provided the following two functions.

PrevCombinationGTmin[ comb, minCombSet ]
PrevCombinationGTmin[ comb_List, minCombSet_List ] :=
  Module[ {newComb = PrevCombination[ comb ]},
    While[ newComb != {} && IsSetLE2[ newComb, minCombSet ],
      newComb = PrevCombination[ newComb ];
    ];
    newComb
  ]

It orderly enumerates downward all the k-conbinations (k = Length[ comb ]) of some elements with repetitions, skipping each newComb such that IsSetLE2[ newComb, minCombSet ].

NextCombinationLTmax[ comb, max, maxCombSet ]
NextCombinationLTmax[ comb_List, max_Integer, maxCombSet_List ] :=
  Module[ {newComb = NextCombination[ comb, max ]},
    While[ newComb != {} && IsSetGE2[ newComb, maxCombSet ],
      newComb = NextCombination[ newComb, max ];
    ];
    newComb
  ]

It orderly enumerates upward all the k-conbinations (k = Length[ comb ]) of max elements with repetitions, skipping each newComb such that IsSetGE2[ newComb, maxCombSet ].

 

FastPartition: Greedy Search

The speed of the FastPartition function is mainly due to the greedy search. To clearly explain how it finds a simple partition, I gave it the form of a train loading scenario.

Problem
  • A railway engine is going to pull a train of wagons.
    • The power of the railway engine is enginePower.
    • The number of wagons is wagonsCount.
    • Each wagon carries one of the loadsWeights.
  • The train has to be loaded as much as possible.
  • FastLoading[ enginePower, wagonsCount, loadsWeights ]
    FastLoading[ enginePower_Integer, wagonsCount_Integer, loadsWeights_List ] :=
      Module[ {heaviestIndex, leftWagons, maxLoad,
          lightest = First[ loadsWeights ], 
          leftPower = enginePower, 
          loading = Table[ 0, {wagonsCount} ]},
        If[ enginePower < wagonsCount lightest,
          Return[ {} ]
        ];
        For[ leftWagons = wagonsCount, leftWagons > 1, leftWagons--,
          maxLoad = leftPower - (leftWagons - 1) lightest;
          heaviestIndex = pBinSearchLE[ maxLoad, loadsWeights ];
          loading[[ leftWagons ]] = heaviestIndex;
          leftPower -= loadsWeights[[ heaviestIndex ]]
        ];
        loading[[ 1 ]] = pBinSearchLE[ leftPower, loadsWeights ];
        Return[ loading ]
      ]  /; IsNatural[ {enginePower, wagonsCount, loadsWeights} ] && loadsWeights != {} && loadsWeights == Union[ loadsWeights ]
    

    It gives a (greedy) wagonsCount-combination of the loadsWeights such that the railway engine will waste minimal enginePower.