/*
 This mex function, randomArray_qdMac, returns an array of numbers generated by the old Macintosh quickdraw function Random().
 INPUT:  [size], randomSeed <optional>
 OUTPUT: ND array of int16, with requested dimensions
 -by xaq pitkow, 19dec07

 Based on algorithm described by Mike Engber:
 random_qdMac() is a reconstruction of the old Quickdraw Random() used on the macintosh.
 Random() generates random numbers in the range [-32767,32767] using a 32-bit seed stored in the QuickDraw global randSeed.
 Each time it is called, it does the following:
 1. Multiply randSeed by 0x41A7 to form a 47-bit product.
 2. Take bits 31 through 46 (MSB) of the product and add them as an unsigned value to bits 0 (LSB) through 30 of the product to form the new seed.
 3. Return the low sixteen bits of the product as the random number, except when they are the value 0x8000, in which case return zero.
*/

#include "mex.h"
#include "matrix.h"
#include "time.h"

#define A       0x41A7
#define HI(x)   ((0xffff0000 & (x))>>16)
#define LO(x)   (0xffff & (x))

static unsigned long randSeed;
static unsigned char initialized = 0;
short Random_qdMac(void);


void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{   
    short *outArray;

    double *dimsDouble;
    mwSize *dimsInt;
    mxArray *dimsMXArray = NULL;
    short ndims;
	mwIndex index;
    long i,nelems=1;
    
    if (!initialized)
    {
        time_t t1;
        time(&t1);
        randSeed = (unsigned long) t1;
        initialized = 1;
    }
    if (nrhs==0)
        mexErrMsgTxt("randomArray_qdMac requires a size vector as the first argument.");
    
    if (nrhs>0)
    {
        ndims = (mwSize) mxGetNumberOfElements(prhs[0]); // first argument is size vector to create
        dimsMXArray = mxCreateNumericMatrix(ndims,1,mxUINT32_CLASS, mxREAL);
        dimsInt = (mwSize *)mxGetData(dimsMXArray);
        dimsDouble = mxGetPr(prhs[0]);
        for (index=0;index<ndims;index++)
        {
            dimsInt[index] = (mwSize)dimsDouble[index];
            nelems *= (long)dimsInt[index];
        }
        plhs[0] = mxCreateNumericArray(ndims, dimsInt, mxINT16_CLASS, mxREAL); // Allocate memory for real short ints
        outArray = (short *)mxGetData(plhs[0]); // point to the data to be filled
    }
    if (nrhs>1)
        randSeed = (unsigned long) *mxGetPr(prhs[1]);  // second argument is randomSeed 
        
    for (i=0;i<nelems;i++)
        outArray[i] = Random_qdMac();

    mxDestroyArray(dimsMXArray);
    return;
}



short Random_qdMac()
{
    short r;
    unsigned long temp;
    // ---------- 32-bit implementation
    temp = A*HI(randSeed) + HI(A * LO(randSeed));
    randSeed = ((temp & 0x7fff) << 16) + HI(temp<<1) + LO(A*randSeed);
    temp = LO(randSeed);
    if (temp == 0x8000)
        r = 0;
    else
        r = (short)temp;
    return r;
}