Using R to solve the Lucky 26 game

Question

I am trying to show my son how coding can be used to solve a problem posed by a game as well as seeing how R handles big data. The game in question is called "Lucky 26". In this game numbers (1-12 with no duplicates) are positioned on 12 points on a star of david (6 vertex, 6 intersections) and the 6 lines of 4 numbers must all add to 26. Of the approximately 479 million possibilities (12P12) there are apparently 144 solutions. I tried to code this in R as follows but memory is an issue it seems. I would greatly appreciate any advice to advance the answer if members have time. Thanking members in advance.

library(gtools)

x=c()
elements <- 12
for (i in 1:elements)
{ 
    x[i]<-i
}

soln=c()            

y<-permutations(n=elements,r=elements,v=x)  
j<-nrow(y)
for (i in 1:j) 
{
L1 <- y[i,1] + y[i,3] + y[i,6] + y[i,8]
L2 <- y[i,1] + y[i,4] + y[i,7] + y[i,11]
L3 <- y[i,8] + y[i,9] + y[i,10] + y[i,11]
L4 <- y[i,2] + y[i,3] + y[i,4] + y[i,5]
L5 <- y[i,2] + y[i,6] + y[i,9] + y[i,12]
L6 <- y[i,5] + y[i,7] + y[i,10] + y[i,12]
soln[i] <- (L1 == 26)&(L2 == 26)&(L3 == 26)&(L4 == 26)&(L5 == 26)&(L6 == 26) 
}

z<-which(soln)
z

Answer 1

For permutations, rcppalgos is great. Unfortunately, there are 479 million possibilities with 12 fields which means that takes up too much memory for most people:

library(RcppAlgos)
elements <- 12
permuteGeneral(elements, elements)
#> Error: cannot allocate vector of size 21.4 Gb

There are some alternatives.

1. Take a sample of the permutations. Meaning, only do 1 million instead of 479 million. To do this, you can use permuteSample(12, 12, n = 1e6). See @JosephWood's answer for a somewhat similar approach except he samples out to 479 million permutations ;)

2. Build a loop in rcpp to evaluate the permutation on creation. This saves memory because you would end up building the function to return only the correct results.

3. Approach the problem with a different algorithm. I will focus on this option.

New algorithm w/ constraints

Segments should be 26

We know that each line segment in the star above needs to add up to 26. We can add that constraint to generating our permutations - give us only combinations that add up to 26:

# only certain combinations will add to 26
lucky_combo <- comboGeneral(12, 4, comparisonFun = '==', constraintFun = 'sum', limitConstraints = 26L)

ABCD and EFGH groups

In the star above, I have colored three groups differently : ABCD, EFGH, and IJLK. The first two groups also have no points in common and are also on line segments of interest. Therefore, we can add another constraint: for combinations that add up to 26, we need to ensure ABCD and EFGH have no number overlap. IJLK will be assigned the remaining 4 numbers.

library(RcppAlgos)
lucky_combo <- comboGeneral(12, 4, comparisonFun = '==', constraintFun = 'sum', limitConstraints = 26L)
two_combo <- comboGeneral(nrow(lucky_combo), 2)

unique_combos <- !apply(cbind(lucky_combo[two_combo[, 1], ], lucky_combo[two_combo[, 2], ]), 1, anyDuplicated)

grp1 <- lucky_combo[two_combo[unique_combos, 1],]
grp2 <- lucky_combo[two_combo[unique_combos, 2],]
grp3 <- t(apply(cbind(grp1, grp2), 1, function(x) setdiff(1:12, x)))

Permute through the groups

We need to find all permutations of each group. That is, we only have combinations that add up to 26. For example, we need to take 1, 2, 11, 12 and make 1, 2, 12, 11; 1, 12, 2, 11; ....

#create group perms (i.e., we need all permutations of grp1, grp2, and grp3)
n <- 4
grp_perms <- permuteGeneral(n, n)
n_perm <- nrow(grp_perms)

# We create all of the permutations of grp1. Then we have to repeat grp1 permutations
# for all grp2 permutations and then we need to repeat one more time for grp3 permutations.
stars <- cbind(do.call(rbind, lapply(asplit(grp1, 1), function(x) matrix(x[grp_perms], ncol = n)))[rep(seq_len(sum(unique_combos) * n_perm), each = n_perm^2), ],
           do.call(rbind, lapply(asplit(grp2, 1), function(x) matrix(x[grp_perms], ncol = n)[rep(1:n_perm, n_perm), ]))[rep(seq_len(sum(unique_combos) * n_perm^2), each = n_perm), ],
           do.call(rbind, lapply(asplit(grp3, 1), function(x) matrix(x[grp_perms], ncol = n)[rep(1:n_perm, n_perm^2), ])))

colnames(stars) <- LETTERS[1:12]

Final Calculations

The last step is to do the math. I use lapply() and Reduce() here to do more functional programming - otherwise, a lot of code would be typed six times. See the original solution for a more thorough explanation of the math code.

# creating a list will simplify our math as we can use Reduce()
col_ind <- list(c('A', 'B', 'C', 'D'), #these two will always be 26
                c('E', 'F', 'G', 'H'),  #these two will always be 26
                c('I', 'C', 'J', 'H'), 
                c('D', 'J', 'G', 'K'),
                c('K', 'F', 'L', 'A'),
                c('E', 'L', 'B', 'I'))

# Determine which permutations result in a lucky star
L <- lapply(col_ind, function(cols) rowSums(stars[, cols]) == 26)
soln <- Reduce(`&`, L)

# A couple of ways to analyze the result
rbind(stars[which(soln),], stars[which(soln), c(1,8, 9, 10, 11, 6, 7, 2, 3, 4, 5, 12)])
table(Reduce('+', L)) * 2

      2       3       4       6 
2090304  493824   69120     960 

Swapping ABCD and EFGH

At the end of the code above, I took advantage that we can swap ABCD and EFGH to get the remaining permutations. Here is the code to confirm that yes, we can swap the two groups and be correct:

# swap grp1 and grp2
stars2 <- stars[, c('E', 'F', 'G', 'H', 'A', 'B', 'C', 'D', 'I', 'J', 'K', 'L')]

# do the calculations again
L2 <- lapply(col_ind, function(cols) rowSums(stars2[, cols]) == 26)
soln2 <- Reduce(`&`, L2)

identical(soln, soln2)
#[1] TRUE

#show that col_ind[1:2] always equal 26:
sapply(L, all)

[1]  TRUE  TRUE FALSE FALSE FALSE FALSE

Performance

In the end, we evaluated only 1.3 million of the 479 permutations and only only shuffled through 550 MB of RAM. It takes around 0.7s to run

# A tibble: 1 x 13
  expression   min median `itr/sec` mem_alloc `gc/sec` n_itr  n_gc
  <bch:expr> <bch> <bch:>     <dbl> <bch:byt>    <dbl> <int> <dbl>
1 new_algo   688ms  688ms      1.45     550MB     7.27     1     5

Answer 2

There are actually 960 solutions. Below we make use of Rcpp, RcppAlgos^*, and the parallel package to obtain the solution in just over 6 seconds using 4 cores. Even if you choose to use a single threaded approach with base R's lapply, the solution is returned in around 25 seconds.

First, we write a simple algorithm in C++ that checks a particular permutation. You will note that we use one array to store all six lines. This is for performance as we utilize cache memory more effectively than using 6 individual arrays. You will also have to keep in mind that C++ uses zero based indexing.

#include <Rcpp.h>
using namespace Rcpp;
// [[Rcpp::plugins(cpp11)]]

constexpr int index26[24] = {0, 2, 5, 7,
                             0, 3, 6, 10,
                             7, 8, 9, 10,
                             1, 2, 3, 4,
                             1, 5, 8, 11,
                             4, 6, 9, 11};

// [[Rcpp::export]]
IntegerVector DavidIndex(IntegerMatrix mat) {
    const int nRows = mat.nrow();
    std::vector<int> res;

    for (int i = 0; i < nRows; ++i) {
        int lucky = 0;

        for (int j = 0, s = 0, e = 4;
             j < 6 && j == lucky; ++j, s += 4, e += 4) {

            int sum = 0;

            for (int k = s; k < e; ++k)
                sum += mat(i, index26[k]);

            lucky += (sum == 26);
        }

        if (lucky == 6) res.push_back(i);
    }

    return wrap(res);
}

Now, using the lower and upper arguments in permuteGeneral, we can generate chunks of permutations and test these individually to keep memory in check. Below, I have chosen to test about 4.7 million permutations at a time. The output gives the lexicographical indices of the permutations of 12! such that the Lucky 26 condition is satisfied.

library(RcppAlgos)
## N.B. 4790016L evenly divides 12!, so there is no need to check
## the upper bound on the last iteration below

system.time(solution <- do.call(c, parallel::mclapply(seq(1L, factorial(12), 4790016L), function(x) {
    perms <- permuteGeneral(12, 12, lower = x, upper = x + 4790015)
    ind <- DavidIndex(perms)
    ind + x
}, mc.cores = 4)))

  user  system elapsed 
13.005   6.258   6.644

## Foregoing the parallel package and simply using lapply,
## we obtain the solution in about 25 seconds:
##   user  system elapsed 
## 18.495   6.221  24.729

Now, we verify using permuteSample and the argument sampleVec which allows you to generate specific permutations (e.g. if you pass 1, it will give you the first permutation (i.e. 1:12)).

system.time(Lucky26 <- permuteSample(12, 12, sampleVec=solution))
 user  system elapsed 
0.001   0.000   0.001

head(Lucky26)
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
[1,]    1    2    4   12    8   10    6   11    5     3     7     9
[2,]    1    2    6   10    8   12    4    7    3     5    11     9
[3,]    1    2    7   11    6    8    5   10    4     3     9    12
[4,]    1    2    7   12    5   10    4    8    3     6     9    11
[5,]    1    2    8    9    7   11    4    6    3     5    12    10
[6,]    1    2    8   10    6   12    4    5    3     7    11     9

tail(Lucky26)
       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
[955,]   12   11    5    3    7    1    9    8   10     6     2     4
[956,]   12   11    5    4    6    2    9    7   10     8     1     3
[957,]   12   11    6    1    8    3    9    5   10     7     4     2
[958,]   12   11    6    2    7    5    8    3    9    10     4     1
[959,]   12   11    7    3    5    1    9    6   10     8     2     4
[960,]   12   11    9    1    5    3    7    2    8    10     6     4

Finally, we verify our solution with base R rowSums:

all(rowSums(Lucky26[, c(1, 3, 6, 8]) == 26)
[1] TRUE

all(rowSums(Lucky26[, c(1, 4, 7, 11)]) == 26)
[1] TRUE

all(rowSums(Lucky26[, c(8, 9, 10, 11)]) == 26)
[1] TRUE

all(rowSums(Lucky26[, c(2, 3, 4, 5)]) == 26)
[1] TRUE

all(rowSums(Lucky26[, c(2, 6, 9, 12)]) == 26)
[1] TRUE

all(rowSums(Lucky26[, c(5, 7, 10, 12)]) == 26)
[1] TRUE

^* I am the author of RcppAlgos

Answer 3

Here's the solution for the little fella:

numbersToDrawnFrom = 1:12
bling=0

while(T==T){

  bling=bling+1
  x=sample(numbersToDrawnFrom,12,replace = F)

  A<-x[1]+x[2]+x[3]+x[4] == 26
  B<-x[4]+x[5]+x[6]+x[7] == 26
  C<-x[7] + x[8] + x[9] + x[1] == 26
  D<-x[10] + x[2] + x[9] + x[11] == 26
  E<-x[10] + x[3] + x[5] + x[12] == 26
  F1<-x[12] + x[6] + x[8] + x[11] == 26

  vectorTrue <- c(A,B,C,D,E,F1)

  if(min(vectorTrue)==1){break}
  if(bling == 1000000){break}

}

x
vectorTrue

Answer 4

Here is another approach. It's based on a MathWorks blog post by Cleve Moler, the author of the first MATLAB.

In the blog post, to save memory the author permutes only 10 elements, keeping the first element as the apex element and the 7th as the base element. Therefore, only 10! == 3628800 permutations need to be tested.
In the code below,

1. Generate the permutations of elements 1 to 10. There are a total of 10! == 3628800 of them.
2. Choose 11 as the apex element and keep it fixed. It really doesn't matter where the assignments start, the other elements will be in the right relative positions.
3. Then assign the 12th element to the 2nd position, 3rd position, etc, in a for loop.

This should produce most of the solutions, give or take rotations and reflections. But it doesn't guarantee that the solutions are unique. It is also reasonably fast.

elements <- 12
x <- seq_len(elements)
p <- gtools::permutations(n = elements - 2, r = elements - 2, v = x[1:10])  

i1 <- c(1, 3, 6, 8)
i2 <- c(1, 4, 7, 11)
i3 <- c(8, 9, 10, 11)
i4 <- c(2, 3, 4, 5)
i5 <- c(2, 6, 9, 12)
i6 <- c(5, 7, 10, 12)

result <- vector("list", elements - 1)
for(i in 0:10){
  if(i < 1){
    p2 <- cbind(11, 12, p)
  }else if(i == 10){
    p2 <- cbind(11, p, 12)
  }else{
    p2 <- cbind(11, p[, 1:i], 12, p[, (i + 1):10])
  }
  L1 <- rowSums(p2[, i1]) == 26
  L2 <- rowSums(p2[, i2]) == 26
  L3 <- rowSums(p2[, i3]) == 26
  L4 <- rowSums(p2[, i4]) == 26
  L5 <- rowSums(p2[, i5]) == 26
  L6 <- rowSums(p2[, i6]) == 26

  i_sol <- which(L1 & L2 & L3 & L4 & L5 & L6)
  result[[i + 1]] <- if(length(i_sol) > 0) p2[i_sol, ] else NA
}
result <- do.call(rbind, result)
dim(result)
#[1] 82 12

head(result)
#     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
#[1,]   11   12    1    3   10    5    8    9    7     6     4     2
#[2,]   11   12    1    3   10    8    5    6    4     9     7     2
#[3,]   11   12    1    7    6    4    3   10    2     9     5     8
#[4,]   11   12    3    2    9    8    6    4    5    10     7     1
#[5,]   11   12    3    5    6    2    9   10    8     7     1     4
#[6,]   11   12    3    6    5    4    2    8    1    10     7     9