How to access elements of GF array in C++?

Hello everyone,

There is a problem with getting elements from GF array in C++. For example, I create GF array:

   gf([1 2], 5) % actually equals to [gf(1,5) gf(2,5)]

Then I want to read, for example, the first element. But I can't, I can only read whole array. There is no gf type in MATLAB API, so mxGetData won't work. gf type has some properties:

   uint32_t x          % array of elements
   double   m          % field extention (I don't know why it's double)
   uint32_t prim_poly  % decimal representation of primitive polynomial

Array elements are stored in x, so it is possible to use mxGetProperty and convert each element back to gf. In this case, I'll have to call gf constructor a lot, but I'm trying to speed up my program by this C++ implementation.

Also, it seems that I have to do all numeric operations by calling them as MATLAB operators (like mexCallMATLAB(1,&hurr,2,&durr,"+"))

Is there any possible way to access one element of GF array explicitly?

MATLAB code:

   function [ f_x ] = lagrange( rts, func )

   rts_N = length(rts);
   func_N = length(func);

   f_x = 0;
   for i=1:func_N
       b_i = rts(i);
       buf = 1;
       for j=1:rts_N
           b_j = rts(j);
           if b_i ~= b_j
               buf = conv(buf, [1 -b_j]/(b_i - b_j));
           end
       end
       f_x = f_x + func(i)*buf;
   end

   end

C++ code (some bugs here):

   void lagrange(mxArray *rts, mxArray *func) {
        size_t rts_N = mxGetN(mxGetProperty(rts, 0, "x"));
  size_t func_N = mxGetN(mxGetProperty(func, 0, "x"));
  size_t func_dims[] = {1, func_N};
  size_t rts_dims[] = {1, rts_N};

        mxArray *m = mxGetProperty(rts, 0, "m");
  if (m == NULL) {mexPrintf("Error! m is NULL \n"); return;}

        mxArray *f_x, *prhs[2]; prhs[0] = mxCreateDoubleScalar(0); prhs[1] = m;
  if (prhs[0] == NULL) {mexPrintf("Error! create f_X: prhs[0]\n"); return;}
        if (prhs[1] == NULL) {mexPrintf("Error! create f_X: prhs[1]\n"); return;}
  
        if (mexCallMATLAB(1,&f_x,2,prhs,"gf") != NULL) {mexPrintf("Error! create f_x\n"); return;}

        uint32_t *rts_numbers = mxGetUint32s(mxGetProperty(rts, 0, "x"));
        if (rts_numbers == NULL) {mexPrintf("Error! create rts_numbers\n"); return;}

        uint32_t *func_numbers = mxGetUint32s(mxGetProperty(func, 0, "x"));
        if (rts_numbers == NULL) {mexPrintf("Error! create rts_numbers\n");  return;}

        mexPrintf("rts_numbers: %d %d %d\n", rts_numbers[0], rts_numbers[1], func_N);

        for (int i = 0; i < func_N; i++) {
    //mexCallMATLAB(0,NULL,1,&f_x,"disp");
    mxArray *b_i; prhs[0] = prhs[0] = mxCreateDoubleScalar(rts_numbers[i]); prhs[1] = m; 
    
    if (mexCallMATLAB(1,&b_i,2,prhs,"gf") != NULL) {mexPrintf("Error! create b_i\n"); return;}
    mexCallMATLAB(0,NULL,1,&b_i,"disp");
    
    mxArray *buf = mxCreateDoubleScalar(0);

                for (int j = 0; j < rts_N; j++) {
      mxArray *b_j; prhs[0] = mxCreateDoubleScalar(rts_numbers[j]); prhs[1] = m;
      if (mexCallMATLAB(1,&b_j,2,prhs,"gf") != NULL) {mexPrintf("Error! create b_i\n"); return;}

                        if (rts_numbers[i] != rts_numbers[j]) {
        mxArray *buf1, *buf2, *buf3; 

                                prhs[0] = b_i; prhs[1] = b_j; mexCallMATLAB(1,&buf2,2,prhs,"-");
                                if (buf2 == NULL) {mexPrintf("Error! b_i - b_j\n"); return;}

                                mxArray* arr = mxCreateDoubleMatrix((mwSize) 1, (mwSize) 2, mxREAL);
                                double* d_arr = new double[2]; d_arr[0] = 1; d_arr[1] = rts_numbers[j];
                                mxSetData(arr, d_arr);
        
                                prhs[0] = arr; prhs[1] = m; mexCallMATLAB(1, &buf1, 2, prhs, "gf");
                                if (buf1 == NULL) {mexPrintf("Error! [1 b_j] to gf\n"); return;}

                                prhs[0] = buf1; prhs[1] = buf2; mexCallMATLAB(1, &buf3, 2, prhs, "/");
                                if (buf3 == NULL) {mexPrintf("Error! [1 b_j]/(b_i - b_j)\n"); return;}
        
                                buf1 = buf; buf2 = buf3;
                                prhs[0] = buf1; prhs[1] = buf2; mexCallMATLAB(1, &buf, 2, prhs, "conv");
                                if (buf == NULL) {mexPrintf("Error! conv(buf, [1 b_j]/(b_i - b_j))\n"); return;}
      }
    }

                prhs[0] = mxCreateDoubleScalar(func_numbers[i]); prhs[1] = buf; 

                mxArray* f_x_buf;
                if (mexCallMATLAB(1, &f_x_buf, 2, prhs, "*") != NULL) {mexPrintf("Error! func_i*buf\n"); return;}

                prhs[0] = f_x; prhs[1] = f_x_buf;
                if (mexCallMATLAB(1, &f_x, 2, prhs, "+") != NULL) {mexPrintf("Error! f_x = f_x + func_i*buf;\n"); return;}
  }
        mexCallMATLAB(0,NULL,1,&f_x,"disp");
   }

Answer 1

[James Tursa]

"... but I'm trying to speed up my program by this C++ implementation ..."

The MATLAB API has you hamstrung in this effort. A gf variable is a type of classdef OOP object. The only things you can do with these in a mex routine via official API functions are mxGetProperty and mxSetProperty, which both work with deep copies, not pointers. Other than that you are stuck with mexCallMATLAB. All of this is going to slow you down. There are unofficial mxGetPropertyShared and mxPutPropertyShared which you might be able to use to avoid the deep copies, or look into these FEX submissions:

https://www.mathworks.com/matlabcentral/fileexchange/30672-mxgetpropertyptr-mxsetpropertysdc-c-mex-functions

If you can, I would advise doing all of the gf stuff outside the mex routine. E.g., pass in the individual gf properties to your mex routine, work with these as regular variables inside your mex routine, then pass the new properties out and build your resulting gf object at the m-code level.

[Nikita Zeulin]

James, thank you for youn answer and your code. It seems that I can't actually access single element of gf array, because the array itself is a gf object with corresponding integer values in x. You write that using these functions makes sense only if there is a "heavy" property you want to read, and I'm going to have only 10-20 integers.

The thing is that all operations must be done in finite fields, so the only way I see to make it work faster is to write this gf thing as a C++ class and overload all numeric operators.