Help convert a C code to AutoIt

[andybiochem]

Hi,

I want to convert the following code to AutoIt to generate random numbers based on the normal distribution.

An Improved Ziggurat Method to Generate Normal Random Samples

JURGEN A DOORNIK

University of Oxford

expandcollapse popup

static double DRanNormalTail(double dMin, int iNegative)
{
    double x, y;
    do
    { x = log(DRanU()) / dMin;
    y = log(DRanU());
    } while (-2 * y < x * x);

return iNegative ? x - dMin : dMin - x;
}


#define ZIGNOR_C 128 /* number of blocks */
#define ZIGNOR_R 3.442619855899 /* start of the right tail */
/* (R * phi(R) + Pr(X>=R)) * sqrt(2\pi) */
#define ZIGNOR_V 9.91256303526217e-3
/* s_adZigX holds coordinates, such that each rectangle has*/
/* same area; s_adZigR holds s_adZigX[i + 1] / s_adZigX[i] */
static double s_adZigX[ZIGNOR_C + 1], s_adZigR[ZIGNOR_C];
static void zigNorInit(int iC, double dR, double dV)

{
    int i; double f;
    f = exp(-0.5 * dR * dR);
    s_adZigX[0] = dV / f; /* [0] is bottom block: V / f(R) */
    s_adZigX[1] = dR;
    s_adZigX[iC] = 0;
    for (i = 2; i < iC; ++i)
    {
        s_adZigX[i] = sqrt(-2 * log(dV / s_adZigX[i - 1] + f));
        f = exp(-0.5 * s_adZigX[i] * s_adZigX[i]);
    }
    for (i = 0; i < iC; ++i)
    s_adZigR[i] = s_adZigX[i + 1] / s_adZigX[i];
}

double DRanNormalZig(void)

{
    unsigned int i;
    double x, u, f0, f1;
    for (;;)
    {
        u = 2 * DRanU() - 1;
        i = IRanU() & 0x7F;
        /* first try the rectangular boxes */
        if (fabs(u) < s_adZigR[i])
        return u * s_adZigX[i];
        /* bottom box: sample from the tail */
        if (i == 0)
        return DRanNormalTail(ZIGNOR_R, u < 0);
        /* is this a sample from the wedges? */
        x = u * s_adZigX[i];
        f0 = exp(-0.5 * (s_adZigX[i] * s_adZigX[i] - x * x) );
        f1 = exp(-0.5 * (s_adZigX[i+1] * s_adZigX[i+1] - x * x) );
        Improved Ziggurat Method · 9
        if (f1 + DRanU() * (f0 - f1) < 1.0)
        return x;
    }
}

Considering I know nothing about C, an someone help me with the following questions...

1) what is the use of the ? and : in this:

return iNegative ? x - dMin : dMin - x;

So far for the first DRanNormalTail func I've got:

Func DRanNormalTail($dMin,$iNegative)
    Do
        $x = Log(Random()) / $dMin
        $y = Log(Random())
    Until Not((-2 * $y) < ($x^2))
    
    Return;what goes here?
EndFunc

2) Is text surrounded by /* */ a comment???

3) is this declaring an array?:

static double s_adZigX[ZIGNOR_C + 1]

4) what is the equivalent of a static void in AI?

5) what is the effect of using a math operator in a function call?

return DRanNormalTail(ZIGNOR_R, u < 0)

I think I could do this on my own, so would prefer it if people could just give me pointers rather than convert it straight out for me!!

Thanks for any help you can give!!!

p.s. why can't ALL languages be as straight-forward as autoit?!?!? ha ha

[martin]

1) what is the use of the ? and : in this:

return iNegative ? x - dMin : dMin - x;

So far for the first DRanNormalTail func I've got:

Func DRanNormalTail($dMin,$iNegative)
    Do
        $x = Log(Random()) / $dMin
        $y = Log(Random())
    Until Not((-2 * $y) < ($x^2))
    
    Return;what goes here?
EndFunc

In c

return statement?x:y

is a realy neat way (IMO) of doing this in AutoIt

If statement then
 return X
else
 return Y
endif

2) Is text surrounded by /* */ a comment???

Yes

3) is this declaring an array?:

static double s_adZigX[ZIGNOR_C + 1]

Yes

4) what is the equivalent of a static void in AI?

Global

well almost. Static variables declared inside a function will only be initialized the first time the function is called, so the next time the function is called it will have its last value. It isn't really global because it is still local to the function and can't be seen outside. So if you convert to AutoIt you need to decare it before the function is ever called and not again, and you need to make sure that there are no other Global variables with the same name. Then it will behave the same. You could instead add an extra parameter which is used ByRef in the function and then have the variable declared in the calling function.

A static function is not something I understand at the moment.

5) what is the effect of using a math operator in a function call?

return DRanNormalTail(ZIGNOR_R, u < 0)

Same as in AutoIt, u < 0 returns true or false so no need to convert that bit for AutoIt.

p.s. why can't ALL languages be as straight-forward as autoit?

IMO no languages are straightforward until you learn to use them and then they're obvious. But yes, AutoIt is definitely easier to learn than most. Edited August 9, 2008 by martin

[ame1011]

Problem is that dranu() and iranu() aren't included.

I did a little search and this is what's missing:

CODE

/*! @file

*==========================================================================

* This code is Copyright © 2005, Jurgen A. Doornik.

* Permission to use this code for non-commercial purposes

* is hereby given, provided proper reference is made to:

* Doornik, J.A. (2005), "An Improved Ziggurat Method to Generate Normal

* Random Samples", mimeo, Nuffield College, University of Oxford,

* and www.doornik.com/research.

* or the published version when available.

* This reference is still required when using modified versions of the code.

* This notice should be maintained in modified versions of the code.

* No warranty is given regarding the correctness of this code.

*==========================================================================

*

* @author Jurgen A. Doornik (Creator)

*

* $Author: ejt $

* $Name: $

* $Revision: 1.2 $

* $State: Exp $

* $Date: 2007/01/30 23:17:12 $

*/

#include <limits.h>

#include <math.h>

#include <stdio.h>

#include <string.h>

#include "zigrandom.h"

#ifdef __cplusplus

extern "C" {

#endif

/*---------------------------- GetInitialSeeds -----------------------------*/

void GetInitialSeeds(unsigned int auiSeed[], int cSeed,

unsigned int uiSeed, unsigned int uiMin)

{

int i;

unsigned int s = uiSeed; /* may be 0 */

for (i = 0; i < cSeed; )

{ /* see Knuth p.106, Table 1(16) and Numerical Recipes p.284 (ranqd1)*/

s = 1664525 * s + 1013904223;

if (s <= uiMin)

continue;

auiSeed = s;

++i;

}

}

/*-------------------------- END GetInitialSeeds ---------------------------*/

/*------------------------ George Marsaglia MWC ----------------------------*/

#define MWC_R 256

#define MWC_A LIT_UINT64(809430660)

#define MWC_AI 809430660

#define MWC_C 362436

static unsigned int s_uiStateMWC = MWC_R - 1;

static unsigned int s_uiCarryMWC = MWC_C;

static unsigned int s_auiStateMWC[MWC_R];

void RanSetSeed_MWC8222(int *piSeed, int cSeed)

{

s_uiStateMWC = MWC_R - 1;

s_uiCarryMWC = MWC_C;

if (cSeed == MWC_R)

{

int i;

for (i = 0; i < MWC_R; ++i)

{

s_auiStateMWC = (unsigned int)piSeed;

}

}

else

{

GetInitialSeeds(s_auiStateMWC, MWC_R, piSeed && cSeed ? piSeed[0] : 0, 0);

}

}

unsigned int IRan_MWC8222(void)

{

UINT64 t;

s_uiStateMWC = (s_uiStateMWC + 1) & (MWC_R - 1);

t = MWC_A * s_auiStateMWC[s_uiStateMWC] + s_uiCarryMWC;

s_uiCarryMWC = (unsigned int)(t >> 32);

s_auiStateMWC[s_uiStateMWC] = (unsigned int)t;

return (unsigned int)t;

}

double DRan_MWC8222(void)

{

UINT64 t;

s_uiStateMWC = (s_uiStateMWC + 1) & (MWC_R - 1);

t = MWC_A * s_auiStateMWC[s_uiStateMWC] + s_uiCarryMWC;

s_uiCarryMWC = (unsigned int)(t >> 32);

s_auiStateMWC[s_uiStateMWC] = (unsigned int)t;

return RANDBL_32new(t);

}

void VecIRan_MWC8222(unsigned int *auiRan, int cRan)

{

UINT64 t;

unsigned int carry = s_uiCarryMWC, status = s_uiStateMWC;

for (; cRan > 0; --cRan, ++auiRan)

{

status = (status + 1) & (MWC_R - 1);

t = MWC_A * s_auiStateMWC[status] + carry;

*auiRan = s_auiStateMWC[status] = (unsigned int)t;

carry = (unsigned int)(t >> 32);

}

s_uiCarryMWC = carry;

s_uiStateMWC = status;

}

void VecDRan_MWC8222(double *adRan, int cRan)

{

UINT64 t;

unsigned int carry = s_uiCarryMWC, status = s_uiStateMWC;

for (; cRan > 0; --cRan, ++adRan)

{

status = (status + 1) & (MWC_R - 1);

t = MWC_A * s_auiStateMWC[status] + carry;

s_auiStateMWC[status] = (unsigned int)t;

*adRan = RANDBL_32new(t);

carry = (unsigned int)(t >> 32);

}

s_uiCarryMWC = carry;

s_uiStateMWC = status;

}

/*----------------------- END George Marsaglia MWC -------------------------*/

/*------------------- normal random number generators ----------------------*/

static int s_cNormalInStore = 0; /* > 0 if a normal is in store */

static DRANFUN s_fnDRanu = DRan_MWC8222;

static IRANFUN s_fnIRanu = IRan_MWC8222;

static IVECRANFUN s_fnVecIRanu = VecIRan_MWC8222;

static DVECRANFUN s_fnVecDRanu = VecDRan_MWC8222;

static RANSETSEEDFUN s_fnRanSetSeed = RanSetSeed_MWC8222;

double DRanU(void)

{

return (*s_fnDRanu)();

}

unsigned int IRanU(void)

{

return (*s_fnIRanu)();

}

void RanVecIntU(unsigned int *auiRan, int cRan)

{

(*s_fnVecIRanu)(auiRan, cRan);

}

void RanVecU(double *adRan, int cRan)

{

(*s_fnVecDRanu)(adRan, cRan);

}

//void RanVecU(double *adRan, int cRan)

//{

// int i, j, c, airan[256];

//

// for (; cRan > 0; cRan -= 256)

// {

// c = min(cRan, 256);

// (*s_fnVecIRanu)(airan, c);

// for (j = 0; j < c; ++j)

// *adRan = RANDBL_32new(airan[j]);

// }

//}

void RanSetSeed(int *piSeed, int cSeed)

{

s_cNormalInStore = 0;

(*s_fnRanSetSeed)(piSeed, cSeed);

}

void RanSetRan(const char *sRan)

{

s_cNormalInStore = 0;

if (strcmp(sRan, "MWC8222") == 0)

{

s_fnDRanu = DRan_MWC8222;

s_fnIRanu = IRan_MWC8222;

s_fnVecIRanu = VecIRan_MWC8222;

s_fnRanSetSeed = RanSetSeed_MWC8222;

}

else

{

s_fnDRanu = NULL;

s_fnIRanu = NULL;

s_fnVecIRanu = NULL;

s_fnRanSetSeed = NULL;

}

}

static unsigned int IRanUfromDRanU(void)

{

return (unsigned int)(UINT_MAX * (*s_fnDRanu)());

}

static double DRanUfromIRanU(void)

{

return RANDBL_32new( (*s_fnIRanu)() );

}

void RanSetRanExt(DRANFUN DRanFun, IRANFUN IRanFun, IVECRANFUN IVecRanFun,

DVECRANFUN DVecRanFun, RANSETSEEDFUN RanSetSeedFun)

{

s_fnDRanu = DRanFun ? DRanFun : DRanUfromIRanU;

s_fnIRanu = IRanFun ? IRanFun : IRanUfromDRanU;

s_fnVecIRanu = IVecRanFun;

s_fnVecDRanu = DVecRanFun;

s_fnRanSetSeed = RanSetSeedFun;

}

/*---------------- END uniform random number generators --------------------*/

/*----------------------------- Polar normal RNG ---------------------------*/

#define POLARBLOCK(u1, u2, d) \

do \

{ u1 = (*s_fnDRanu)(); u1 = 2 * u1 - 1;\

u2 = (*s_fnDRanu)(); u2 = 2 * u2 - 1;\

d = u1 * u1 + u2 * u2; \

} while (d >= 1); \

d = sqrt( (-2.0 / d) * log(d) ); \

u1 *= d; u2 *= d

static double s_dNormalInStore;

double DRanNormalPolar(void) /* Polar Marsaglia */

{

double d, u1;

if (s_cNormalInStore)

u1 = s_dNormalInStore, s_cNormalInStore = 0;

else

{

POLARBLOCK(u1, s_dNormalInStore, d);

s_cNormalInStore = 1;

}

return u1;

}

#define FPOLARBLOCK(u1, u2, d) \

do \

{ u1 = (float)((*s_fnDRanu)()); u1 = 2 * u1 - 1;\

u2 = (float)((*s_fnDRanu)()); u2 = 2 * u2 - 1;\

d = u1 * u1 + u2 * u2; \

} while (d >= 1); \

d = sqrt( (-2.0 / d) * log(d) ); \

u1 *= d; u2 *= d

static float s_fNormalInStore;

double FRanNormalPolar(void) /* Polar Marsaglia */

{

float d, u1;

if (s_cNormalInStore)

u1 = s_fNormalInStore, s_cNormalInStore = 0;

else

{

POLARBLOCK(u1, s_fNormalInStore, d);

s_cNormalInStore = 1;

}

return (double)u1;

}

/*--------------------------- END Polar normal RNG -------------------------*/

/*------------------------------ DRanQuanNormal -----------------------------*/

static double dProbN(double x, int fUpper)

{

double p; double y; int fnegative = 0;

if (x < 0)

x = -x, fnegative = 1, fUpper = !fUpper;

else if (x == 0)

return 0.5;

if ( !(x <= 8 || (fUpper && x <= 37) ) )

return (fUpper) ? 0 : 1;

y = x * x / 2;

if (x <= 1.28)

{

p = 0.5 - x * (0.398942280444 - 0.399903438504 * y /

(y + 5.75885480458 - 29.8213557808 /

(y + 2.62433121679 + 48.6959930692 /

(y + 5.92885724438))));

}

else

{

p = 0.398942280385 * exp(-y) /

(x - 3.8052e-8 + 1.00000615302 /

(x + 3.98064794e-4 + 1.98615381364 /

(x - 0.151679116635 + 5.29330324926 /

(x + 4.8385912808 - 15.1508972451 /

(x + 0.742380924027 + 30.789933034 /

(x + 3.99019417011))))));

}

return (fUpper) ? p : 1 - p;

}

double DProbNormal(double x)

{

return dProbN(x, 0);

}

double DRanQuanNormal(void)

{

return DProbNormal(DRanNormalPolar());

}

double FRanQuanNormal(void)

{

return DProbNormal(FRanNormalPolar());

}

/*----------------------------- END DRanQuanNormal -------------------------*/

#ifdef __cplusplus

} //extern "C"

#endif