Matlab 2017a mex of mandelbrot_step.c

In the book "Experiments with Matlab" of Cleve Moler (free ebook)

Cleve Moler : Experiments with matlab

with source code

Experiments with matlab - source code

In chapter "Mandelbrot set" there is a file mandelbrot_step.c which can be compiled by mex.

Now in R2017a this fails even with -compatibleArrayDims from 64 bit api changes

Matlab crashes on it.

What should change in

#include <math.h>
#include "mex.h"
void mexFunction( int nlhs, mxArray *plhs[], 
      int nrhs, const mxArray *prhs[] )

/* function [z,kz] = mandelbrot_step(z,kz,z0,d);
 * Take one step of the Mandelbrot iteration.
 * Complex arithmetic:
 *    z = z.^2 + z0
 *    kz(abs(z) < 2) == d
 * Real arithmetic:
 *    x <-> real(z);
 *    y <-> imag(z);
 *    u <-> real(z0);
 *    v <-> imag(z0);
 *    [x,y] = [x.^2-y.^2+u, 2*x.*y+v];
 *    kz(x.^2+y.^2 < 4) = d;
 */

{ 
    double *x,*y,*u,*v,t; 
    unsigned short *kz,d;
    int j,n; 

      x = mxGetPr(prhs[0]); 
      y = mxGetPi(prhs[0]);
      kz = (unsigned short *) mxGetData(prhs[1]); 
      u = mxGetPr(prhs[2]); 
      v = mxGetPi(prhs[2]);
      d = (unsigned short) mxGetScalar(prhs[3]);
      plhs[0] = prhs[0];
      plhs[1] = prhs[1];

      n = mxGetN(prhs[0]);
      for (j=0; j<n*n; j++) {
          if (kz[j] == d-1) {
              t = x[j];
              x[j] = x[j]*x[j] - y[j]*y[j] + u[j];
              y[j] = 2*t*y[j] + v[j];
              if (x[j]*x[j] + y[j]*y[j] < 4) {
                  kz[j] = d;
              }
          }
      }
      return;
  }

My mex compiler is Microsoft Visual Studio 2015 Community Edition.

Answer 1

[Jan]

Please try this:

unsigned short ==> uint16_T
int            ==> size_t

Writing into the input array and replying it as output is, well, bold. The officially supported way is (James, please correct me on demand):

plhs[0] = mxDuplicateArray(prhs[0]);
plhs[1] = mxDuplicateArray(prhs[1]);
x  = mxGetPr(plhs[0]); 
y  = mxGetPi(plhs[0]);
kz = (uint16_T *) mxGetData(plhs[1]); 

u = mxGetPr(prhs[2]); 
v = mxGetPi(prhs[2]);
d = (unsigned short) mxGetScalar(prhs[3]);

[Bert RAM Aerts]

I changed int into mwSize, as this is advised for arrays bigger than 2^32-1.

Indeed you are perfectly right that plhs[0]=prhs[0] is a tricky statement.

With your corrections it now works just fine.

Thanks!