libcrbn/sampler/faure.hpp

#ifndef _FAURE_HPP_
#define _FAURE_HPP_

// Generate faure permutation matrix
int** faure(const int& iSize) {
        int** permutation;
        int   i,j;
        
        if(iSize < 3) return 0;
        
        permutation = new (int*)[iSize];

        permutation[0] = 0;

        permutation[1] = new int[2];
        permutation[1][0] = 0;
        permutation[1][1] = 1;

        for(i = 3; i < iSize; ++i) {
                permutation[i-1] = new int[i];

                /* odd */       
                if(i & 1) {
                        for(j = 0; j < ((i-1)/2); ++j) {
                                if((2 * permutation[i-2][j]) >= (i-1)) {
                                        permutation[i-1][j] = permutation[i-2][j] + 1;
                                }
                                else {
                                        permutation[i-1][j] = permutation[i-2][j];
                                }
                        }

                        permutation[i-1][(i-1)/2] = (i-1) / 2;

                        for(j = (i+1)/2; j < i; ++j) {
                                if((2 * permutation[i-2][j-1]) >= (i-1)) {
                                        permutation[i-1][j] = permutation[i-2][j-1] + 1;
                                }
                                else {
                                        permutation[i-1][j] = permutation[i-2][j-1];
                                }
                        }
                }
                /* even */
                else {
                        for(j = 0; j < (i/2); ++j) {
                                permutation[i-1][j] = 2 * permutation[(i/2)-1][j];
                        }

                        for(j = i/2; j < i; ++j) {
                                permutation[i-1][j] = permutation[i-1][j-(i/2)] + 1;
                        }
                }
        }
        
        return permutation;
}

#endif

---