Implementing Ryser's formula

[skitini]

Hi all,

I'm currently trying to write a piece of code that gives the number of Latin squares from a formula that was found. One part of the code requires to generate all [0,1] n x n matrices for the given number n, and I have managed to get a bit of code for this. One bit it needs to do for all these values however is then calculate the permanent of each of these matrices.

In doing some research, I managed to find some code that calculates the permanent, which does so by using Ryser's formula. Whilst I understand most of the code, I am unaware how to get or test this code, and I cannot get a vector to work.

For now, I just want to either declare a vector/ n x n matrix in the code and it calculate this permanent, or I enter a vector as a string and it does the same result.

The code is split into two parts (for reference you can check http://www.codeproject.com/Articles/21282/Compute-Permanent-of-a-Matrix-with-Ryser-s-Algorit)

inline int* dec2binarr(long n, int dim)
{
    // note: res[dim] will save the sum res[0]+...+res[dim-1]
    int* res = (int*)calloc(dim + 1, sizeof(int));   
    int pos = dim - 1;

    // note: this will crash if dim < log_2(n)...
    while (n > 0)
    {
        res[pos] = n % 2;
        res[dim] += res[pos];
        n = n / 2; // integer division        
        pos--;
    }

    return res;
}

and

long per(int* A, int n) // expects n by n matrix encoded as vector
{
    long sum = 0;
    long rowsumprod, rowsum;
    int* chi = new int[n + 1];    
    double C = (double)pow((double)2, n); 

    // loop all 2^n submatrices of A
    for (int k = 1; k < C; k++)
    {
        rowsumprod = 1;
        chi = dec2binarr(k, n); // characteristic vector        

        // loop columns of submatrix #k
        for (int m = 0; m < n; m++)
        {
            rowsum = 0;

            // loop rows and compute rowsum
            for (int p = 0; p < n; p++)
                rowsum += chi[p] * A[m * n + p];
        
            // update product of rowsums
            rowsumprod *= rowsum;    
        
            // (optional -- use for sparse matrices)
            // if (rowsumprod == 0) break;    
        }        
        
        sum += (long)pow((double)-1, n - chi[n]) * rowsumprod;        
    }    

    return sum;
} 

I have tried various things in the main, but I am still very new to C++ and don't know how to implement this particular code.

As some background, I am learning C++ as an extension to my maths course, and am just doing this to further my studies. This is a topic we are studying in our group project, and seemed like a perfect project for me to learn, and I learn better from having a task.

Thanks!

[keskiverto]

Horror code.

Something is allocated with new and promptly forgotten without deallocation. Then allocation with calloc (a separate system), which is never deallocated either.

We want that fixed, yes?

First, the chi. Array of bits. There is a std::bitset
http://www.cplusplus.com/reference/bitset/bitset/

However, the chi[n] is not a bit. It is the sum of bits. std::bitset has count() for that.

This is untested:

long per( const std::vector<int> & A, int n )
{
  if ( n*n != A.size() ) {
    std::cout << "A is not a n*n matrix\n";
    return 0;
  }

  unsigned long C = 1;
  C << n; // bitshift equals integer pow() for base2

  // loop all 2^n submatrices of A
  long sum = 0;
  for ( unsigned long k = 1; k < C; ++k ) {
    long rowsumprod = 1;

    const std::bitset<sizeof<unsigned long>> chi( k );
    const auto chilast = chi.size() - 1; // bitset indexing starts from right
    // only the last n elements in chi will be used

    // loop columns of submatrix #k
    for ( int m = 0; m < n; ++m ) {
      long rowsum = 0;
      // loop rows and compute rowsum
      for ( int p = 0; p < n; ++p ) {
         rowsum += chi[chilast - p] * A[m * n + p];
      }
      // update product of rowsums
      rowsumprod *= rowsum;

      // (optional -- use for sparse matrices)
      // if (rowsumprod == 0) break;
    }
    const auto t = n - chi.count(); // n - chi[n]
    const bool odd = t % 2;
    // if t is odd, then -1^t == -1
    // if t is even, then -1^t == 1
    sum += rowsumprod * ( odd ? -1 : 1 );
  }    
  return sum;
}