/**
* These subroutines implement the actual iterators.
*
* The real combinatorics are done in-place on a private array of indices
* that is guaranteed to hold integers. We cannot assume they are IVs though,
* because in a few places in the Perl side there's some simple arithmetic
* that is enough to give NVs in 5.6.x.
*
* Once the next tuple has been computed the corresponding slice of data is
* copied in the Perl side. I tried to slice data here in C but it was in
* fact slightly slower. I think we would need to pass aliases to gain
* some more speed.
*
* All the subroutines return -1 when the sequence has been exhausted.
*/
#include "EXTERN.h"
#include "perl.h"
#include "XSUB.h"
#define SETIV(av, i, n) (sv_setiv(*av_fetch(av, i, 0), n))
#define GETIV(av, i) (SvIV(*av_fetch(av, i, 0)))
#define INCR(av, i) (SETIV(av, i, 1 + GETIV(av, i)))
#define GETAV(avptr) ((AV*) SvRV(avptr))
/**
* Swap the ith and jth elements in av.
*
* Assumes av contains IVs.
*/
void __swap(AV* av, int i, int j)
{
IV tmp = GETIV(av, i);
SETIV(av, i, GETIV(av, j));
SETIV(av, j, tmp);
}
/**
* This implementation emulates what we do by hand. It is faster than
* Algorithm T from [2], which gives another lexicographic ordering.
*/
int __next_combination(SV* tuple_avptr, int max_n)
{
AV* tuple = GETAV(tuple_avptr);
int i, j;
IV n;
I32 offset, len_tuple;
SV* e;
len_tuple = av_len(tuple);
offset = max_n - len_tuple;
for (i = len_tuple; i >= 0; --i) {
e = *av_fetch(tuple, i, 0);
n = SvIV(e);
if (n < i + offset) {
sv_setiv(e, ++n);
for (j = i+1; j <= len_tuple; ++j)
SETIV(tuple, j, ++n);
return i;
}
}
return -1;
}
/**
* This provisional implementation emulates what we do by hand.
*/
int __next_combination_with_repetition(SV* tuple_avptr, int max_n)
{
AV* tuple = GETAV(tuple_avptr);
int i, j;
IV n;
I32 len_tuple;
len_tuple = av_len(tuple);
for (i = len_tuple; i >= 0; --i) {
n = GETIV(tuple, i);
if (n < max_n) {
++n;
for (j = i; j <= len_tuple; ++j)
SETIV(tuple, j, n);
return i;
}
}
return -1;
}
/**
* This provisional implementation emulates what we do by hand, keeping
* and array of booleans (used) to keep track of the indices in use.
* That is, used[n] == 1 if and only if tuple[i] == n for some i.
*
*/
int __next_variation(SV* tuple_avptr, SV* used_avptr, int max_n)
{
AV* tuple = GETAV(tuple_avptr);
AV* used = GETAV(used_avptr);
int i, j;
I32 len_tuple;
SV* e;
IV n;
len_tuple = av_len(tuple);
for (i = len_tuple; i >= 0; --i) {
/* from right to left, find the first position that can be incremented */
e = *av_fetch(tuple, i, 0);
n = SvIV(e);
SETIV(used, n, 0);
while (++n <= max_n) {
if (!GETIV(used, n)) {
/* if we get here we nececessarily exit the subrutine, so forget about the outer while and for */
sv_setiv(e, n);
SETIV(used, n, 1);
for (j = i+1; j <= len_tuple; ++j) {
/* from there to the right, fill the tuple with the lowest available numbers */
n = -1;
while (++n <= max_n) {
if (!GETIV(used, n)) {
SETIV(tuple, j, n);
SETIV(used, n, 1);
break;
}
}
}
return i;
}
}
}
return -1;
}
/**
* This provisional implementation emulates what we do by hand.
*/
int __next_variation_with_repetition(SV* tuple_avptr, int max_n)
{
AV* tuple = GETAV(tuple_avptr);
int i;
I32 len_tuple;
SV* e;
len_tuple = av_len(tuple);
for (i = len_tuple; i >= 0; --i) {
e = *av_fetch(tuple, i, 0);
if (SvIV(e) < max_n) {
sv_setiv(e, 1 + SvIV(e));
return i;
}
sv_setiv(e, 0);
}
return -1;
}
/**
* Algorithm H (Loopless reflected mixed-radix Gray generation), from [1].
*
* [Initialize.] and [Visit.] are done in the Perl side.
*/
int __next_variation_with_repetition_gray_code(SV* tuple_avptr, SV* f_avptr, SV* o_avptr, int max_m)
{
AV* tuple = GETAV(tuple_avptr);
AV* f = GETAV(f_avptr);
AV* o = GETAV(o_avptr);
I32 n;
IV j, aj;
n = av_len(tuple) + 1;
/* [Choose j.] */
j = GETIV(f, 0);
SETIV(f, 0, 0);
/* [Change coordinate j.] */
if (j == n)
return -1;
else
SETIV(tuple, j, GETIV(tuple, j) + GETIV(o, j));
/* [Reflect?] */
aj = GETIV(tuple, j);
if (aj == 0 || aj == max_m) {
SETIV(o, j, -GETIV(o, j));
SETIV(f, j, GETIV(f, j+1));
SETIV(f, j+1, j+1);
}
return j;
}
/**
* Algorithm L (Lexicographic permutation generation), adapted from [1].
* I used "h" instead of the letter "l" for the sake of readability.
*
* This algorithm goes back at least to the 18th century, and has been rediscovered
* ever since.
*/
int __next_permutation(SV* tuple_avptr)
{
AV* tuple = GETAV(tuple_avptr);
I32 max_n, j, h, k;
IV aj;
max_n = av_len(tuple);
/* [Find j.] Find the element a(j) behind the longest decreasing tail. */
for (j = max_n-1; j >= 0 && GETIV(tuple, j) > GETIV(tuple, j+1); --j)
;
if (j == -1)
return -1;
/* [Increase a(j).] Find the rightmost element a(h) greater than a(j) and swap them. */
aj = GETIV(tuple, j);
for (h = max_n; aj > GETIV(tuple, h); --h)
;
__swap(tuple, j, h);
/* [Reverse a(j+1)...a(max_n)] Reverse the tail. */
for (k = j+1, h = max_n; k < h; ++k, --h)
__swap(tuple, k, h);
/* Done. */
return 1;
}
int __next_permutation_heap(SV* a_avptr, SV* c_avptr)
{
AV* a = GETAV(a_avptr);
AV* c = GETAV(c_avptr);
int k;
I32 n;
IV ck;
n = av_len(a) + 1;
for (k = 1, ck = GETIV(c, k); ck == k; ++k, ck = GETIV(c, k))
SETIV(c, k, 0);
if (k == n)
return -1;
++ck;
SETIV(c, k, ck);
k % 2 == 0 ? __swap(a, k, 0) : __swap(a, k, ck-1);
return k;
}
/**
* The only algorithms I have found by now are either recursive, or a
* naive wrapper around permutations() that loops over all of them and
* discards the ones with fixed-points.
*
* We take here a mixed-approach, which consists on starting with the
* algorithm in __next_permutation() and tweak a couple of places that
* allow us to skip a significant number of permutations sometimes.
*
* Benchmarking shows this subroutine makes derangements() more than
* two and a half times faster than permutations() for n = 8.
*/
int __next_derangement(SV* tuple_avptr)
{
AV* tuple = GETAV(tuple_avptr);
I32 max_n, min_j, j, h, k;
IV aj;
max_n = av_len(tuple);
min_j = max_n;
THERE_IS_A_FIXED_POINT:
/* Find the element a(j) behind the longest decreasing tail. */
for (j = max_n-1; j >= 0 && GETIV(tuple, j) > GETIV(tuple, j+1); --j)
;
if (j == -1)
return -1;
if (min_j > j)
min_j = j;
/* Find the rightmost element a(h) greater than a(j) and swap them. */
aj = GETIV(tuple, j);
for (h = max_n; aj > GETIV(tuple, h); --h)
;
__swap(tuple, j, h);
/* If a(h) was j leave the tail in decreasing order and try again. */
if (GETIV(tuple, j) == j)
goto THERE_IS_A_FIXED_POINT;
/* I tried an alternative approach that would in theory avoid the
generation of some permutations with fixed-points: keeping track of
the leftmost fixed-point, and reversing the elements to its right.
But benchmarks up to n = 11 showed no difference whatsoever.
Thus, I left this version, which is simpler.
That n = 11 does not mean there was a difference for n = 12, it
means I stopped benchmarking at n = 11. */
/* Otherwise reverse the tail and return if there's no fixed point. */
for (k = j+1, h = max_n; k < h; ++k, --h)
__swap(tuple, k, h);
for (k = max_n; k > min_j; --k)
if (GETIV(tuple, k) == k)
goto THERE_IS_A_FIXED_POINT;
return 1;
}
/*
* This is a transcription of algorithm 3 from [3].
*
* It is a classical approach based on restricted growth strings, which are
* introduced in the paper.
*/
int __next_partition(SV* k_avptr, SV* M_avptr)
{
AV* k = GETAV(k_avptr); /* follows notation in [3] */
AV* M = GETAV(M_avptr); /* follows notation in [3] */
int i, j;
IV mi;
I32 len_k;
len_k = av_len(k);
for (i = len_k; i > 0; --i) {
if (GETIV(k, i) <= GETIV(M, i-1)) {
INCR(k, i);
if (GETIV(k, i) > GETIV(M, i))
SETIV(M, i, GETIV(k, i));
mi = GETIV(M, i);
for (j = i+1; j <= len_k; ++j) {
SETIV(k, j, 0);
SETIV(M, j, mi);
}
return i;
}
}
return -1;
}
/*
* This is a transcription of algorithm 8 from [3].
*
* It is an adaptation of the previous one.
*/
int __next_partition_of_size_p(SV* k_avptr, SV* M_avptr, int p)
{
AV* k = GETAV(k_avptr); /* follows notation in [3] */
AV* M = GETAV(M_avptr); /* follows notation in [3] */
int i, j;
IV mi, x;
I32 len_k, n_minus_p;
len_k = av_len(k);
for (i = len_k; i > 0; --i) {
if (GETIV(k, i) < p-1 && GETIV(k, i) <= GETIV(M, i-1)) {
INCR(k, i);
if (GETIV(k, i) > GETIV(M, i))
SETIV(M, i, GETIV(k, i));
n_minus_p = len_k + 1 - p;
mi = GETIV(M, i);
x = n_minus_p + mi;
for (j = i+1; j <= x; ++j) {
SETIV(k, j, 0);
SETIV(M, j, mi);
}
for (j = x+1; j <= len_k; ++j) {
SETIV(k, j, j - n_minus_p);
SETIV(M, j, j - n_minus_p);
}
return i;
}
}
return -1;
}
/*
* This subroutine has been copied from List::PowerSet.
*
* It uses a vector of bits "odometer" to indicate which elements to include
* in each iteration. The odometer runs and eventually exhausts all possible
* combinations of 0s and 1s.
*/
AV* __next_subset(SV* data_avptr, SV* odometer_avptr)
{
AV* data = GETAV(data_avptr);
AV* odometer = GETAV(odometer_avptr);
I32 len_data = av_len(data);
AV* subset = newAV();
IV adjust = 1;
int i;
IV n;
for (i = 0; i <= len_data; ++i) {
n = GETIV(odometer, i);
if (n) {
av_push(subset, newSVsv(*av_fetch(data, i, 0)));
}
if (adjust) {
adjust = 1 - n;
SETIV(odometer, i, adjust);
}
}
return (AV*) sv_2mortal((SV*) subset);
}
/** -------------------------------------------------------------------
*
* XS stuff starts here.
*
*/
MODULE = Algorithm::Combinatorics PACKAGE = Algorithm::Combinatorics
PROTOTYPES: DISABLE
int
__next_combination(tuple_avptr, max_n)
SV* tuple_avptr
int max_n
int
__next_combination_with_repetition(tuple_avptr, max_n)
SV* tuple_avptr
int max_n
int
__next_variation(tuple_avptr, used_avptr, max_n)
SV* tuple_avptr
SV* used_avptr
int max_n
int
__next_variation_with_repetition(tuple_avptr, max_n)
SV* tuple_avptr
int max_n
int
__next_variation_with_repetition_gray_code(tuple_avptr, f_avptr, o_avptr, max_m)
SV* tuple_avptr
SV* f_avptr
SV* o_avptr
int max_m
int
__next_permutation(tuple_avptr)
SV* tuple_avptr
int
__next_permutation_heap(a_avptr, c_avptr)
SV* a_avptr
SV* c_avptr
int
__next_derangement(tuple_avptr)
SV* tuple_avptr
int
__next_partition(k_avptr, M_avptr)
SV* k_avptr
SV* M_avptr
int
__next_partition_of_size_p(k_avptr, M_avptr, p)
SV* k_avptr
SV* M_avptr
int p
AV*
__next_subset(data_avptr, odometer_avptr)
SV* data_avptr
SV* odometer_avptr