Mercurial > hg > camir-aes2014

annotate toolboxes/FullBNT-1.0.7/KPMtools/ind2subvKPM.c @ 0:e9a9cd732c1e tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
first hg version after svn
author wolffd
date Tue, 10 Feb 2015 15:05:51 +0000
parents
children
rev   line source
wolffd@0 1 /* C mex version of ind2subv.m in misc directory */
wolffd@0 2 /* 2 input, 1 output */
wolffd@0 3 /* siz, ndx */
wolffd@0 4 /* sub */
wolffd@0 5
wolffd@0 6 #include "mex.h"
wolffd@0 7 #include <math.h>
wolffd@0 8
wolffd@0 9 void rbinary(int num, int n, double *rbits){
wolffd@0 10 int i, mask;
wolffd@0 11 num = num - 1;
wolffd@0 12
wolffd@0 13 mask = 1 << (n-1); /* mask = 00100...0 , where the 1 is in column n (rightmost = col 1) */
wolffd@0 14 for (i = 0; i < n; i++) {
wolffd@0 15 rbits[n-i-1] = ((num & mask) == 0) ? 1 : 2;
wolffd@0 16 num <<= 1;
wolffd@0 17 }
wolffd@0 18 }
wolffd@0 19
wolffd@0 20 void ind_subv(int num, const double *sizes, int n, double *rbits){
wolffd@0 21 int i;
wolffd@0 22 int *cumprod;
wolffd@0 23
wolffd@0 24 cumprod = malloc(n * sizeof(int));
wolffd@0 25 num = num - 1;
wolffd@0 26 cumprod[0] = 1;
wolffd@0 27 for (i = 0; i < n-1; i++)
wolffd@0 28 cumprod[i+1] = cumprod[i] * (int)sizes[i];
wolffd@0 29 for (i = n-1; i >= 0; i--) {
wolffd@0 30 rbits[i] = ((int)floor(num / cumprod[i])) + 1;
wolffd@0 31 num = num % cumprod[i];
wolffd@0 32 }
wolffd@0 33 free(cumprod);
wolffd@0 34 }
wolffd@0 35
wolffd@0 36
wolffd@0 37 void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){
wolffd@0 38 int i, j, k, nCol, nRow, nnRow, binary, count, temp, temp1, start;
wolffd@0 39 double *pSize, *pNdx, *pr;
wolffd@0 40 double ndx;
wolffd@0 41 int *subv, *cumprod, *templai;
wolffd@0 42
wolffd@0 43 pSize = mxGetPr(prhs[0]);
wolffd@0 44 pNdx = mxGetPr(prhs[1]);
wolffd@0 45 nCol = mxGetNumberOfElements(prhs[0]);
wolffd@0 46 nnRow = mxGetNumberOfElements(prhs[1]);
wolffd@0 47
wolffd@0 48 nRow = 1;
wolffd@0 49 for(i=0; i<nCol; i++){
wolffd@0 50 nRow *= (int)pSize[i];
wolffd@0 51 }
wolffd@0 52
wolffd@0 53 if(nCol == 0){
wolffd@0 54 plhs[0] = mxDuplicateArray(prhs[1]);
wolffd@0 55 return;
wolffd@0 56 }
wolffd@0 57
wolffd@0 58 binary = 2;
wolffd@0 59 for (i = 0; i < nCol; i++){
wolffd@0 60 if (pSize[i] > 2.0){
wolffd@0 61 binary = 0;
wolffd@0 62 break;
wolffd@0 63 }
wolffd@0 64 else if((int)pSize[i] == 1){
wolffd@0 65 binary = 1;
wolffd@0 66 }
wolffd@0 67 }
wolffd@0 68
wolffd@0 69 if(nnRow == 1){
wolffd@0 70 ndx = mxGetScalar(prhs[1]);
wolffd@0 71 plhs[0] = mxCreateDoubleMatrix(1, nCol, mxREAL);
wolffd@0 72 pr = mxGetPr(plhs[0]);
wolffd@0 73 if(binary == 2)rbinary((int)ndx, nCol, pr);
wolffd@0 74 else ind_subv((int)ndx, pSize, nCol, pr);
wolffd@0 75 return;
wolffd@0 76 }
wolffd@0 77
wolffd@0 78 plhs[0] = mxCreateDoubleMatrix(nnRow, nCol, mxREAL);
wolffd@0 79 pr = mxGetPr(plhs[0]);
wolffd@0 80
wolffd@0 81 subv = malloc(nRow * nCol * sizeof(int));
wolffd@0 82
wolffd@0 83 if (binary == 2) {
wolffd@0 84 for(j=0; j<nCol; j++){
wolffd@0 85 temp = (1 << j);
wolffd@0 86 temp1 = j * nRow;
wolffd@0 87 for(i=0; i<nRow; i++){
wolffd@0 88 subv[temp1 + i] = ((i & temp) == 0) ? 1 : 2;
wolffd@0 89 }
wolffd@0 90 }
wolffd@0 91 }
wolffd@0 92 else if(binary == 1){
wolffd@0 93 cumprod = (int *)malloc(nCol * sizeof(int));
wolffd@0 94 templai = (int *)malloc(nCol * sizeof(int));
wolffd@0 95 cumprod[0] = 1;
wolffd@0 96 templai[0] = nRow;
wolffd@0 97 for(i=1; i<nCol; i++){
wolffd@0 98 k = (int)pSize[i-1] - 1;
wolffd@0 99 cumprod[i] = cumprod[i-1] << k;
wolffd@0 100 templai[i] = templai[i-1] >> k;
wolffd@0 101 }
wolffd@0 102 for(j=0; j<nCol; j++){
wolffd@0 103 temp1 = j * nRow;
wolffd@0 104 if(pSize[j] == 1.0){
wolffd@0 105 for(i=0; i<nRow; i++){
wolffd@0 106 subv[temp1 + i] = 1;
wolffd@0 107 }
wolffd@0 108 }
wolffd@0 109 else{
wolffd@0 110 temp = 1;
wolffd@0 111 count = 0;
wolffd@0 112 for(i=0; i<templai[j]; i++){
wolffd@0 113 if(temp > 2) temp = 1;
wolffd@0 114 for(k=0; k<cumprod[j]; k++){
wolffd@0 115 subv[temp1 + count] = temp;
wolffd@0 116 count++;
wolffd@0 117 }
wolffd@0 118 temp++;
wolffd@0 119 }
wolffd@0 120 }
wolffd@0 121 }
wolffd@0 122 free(templai);
wolffd@0 123 free(cumprod);
wolffd@0 124 }
wolffd@0 125 else {
wolffd@0 126 cumprod = (int *)malloc(nCol * sizeof(int));
wolffd@0 127 templai = (int *)malloc(nCol * sizeof(int));
wolffd@0 128 cumprod[0] = 1;
wolffd@0 129 for(i=1; i<nCol; i++){
wolffd@0 130 cumprod[i] = (int)(cumprod[i-1] * pSize[i-1]);
wolffd@0 131 }
wolffd@0 132 templai[0] = nRow;
wolffd@0 133 for(i=1; i<nCol; i++){
wolffd@0 134 templai[i] = nRow / cumprod[i];
wolffd@0 135 }
wolffd@0 136 for(j=0; j<nCol; j++){
wolffd@0 137 temp1 = j * nRow;
wolffd@0 138 if(pSize[j] == 1.0){
wolffd@0 139 for(i=0; i<nRow; i++){
wolffd@0 140 subv[temp1 + i] = 1;
wolffd@0 141 }
wolffd@0 142 }
wolffd@0 143 else{
wolffd@0 144 temp = 1;
wolffd@0 145 count = 0;
wolffd@0 146 for(i=0; i<templai[j]; i++){
wolffd@0 147 if(temp > (int)pSize[j]) temp = 1;
wolffd@0 148 for(k=0; k<cumprod[j]; k++){
wolffd@0 149 subv[temp1 + count] = temp;
wolffd@0 150 count++;
wolffd@0 151 }
wolffd@0 152 temp++;
wolffd@0 153 }
wolffd@0 154 }
wolffd@0 155 }
wolffd@0 156 free(cumprod);
wolffd@0 157 free(templai);
wolffd@0 158 }
wolffd@0 159
wolffd@0 160 count = 0;
wolffd@0 161 for(j=0; j<nCol; j++){
wolffd@0 162 temp1 = j * nRow;
wolffd@0 163 for(i=0; i<nnRow; i++){
wolffd@0 164 start = (int)pNdx[i] - 1;
wolffd@0 165 pr[count] = subv[temp1 + start];
wolffd@0 166 count++;
wolffd@0 167 }
wolffd@0 168 }
wolffd@0 169 free(subv);
wolffd@0 170 }
wolffd@0 171
wolffd@0 172
wolffd@0 173
wolffd@0 174
wolffd@0 175
wolffd@0 176