How to use Numerical Recipes?

[johnny92]

I am having trouble generating a normal distributed number using the Box-Muller Method in Numerical Recipes. I tried declaring in my main section a norm_dev as a Normaldev_BM type. Then printing, but no luck! Please help

#include "nr3.h" 

struct Ran {
	Ullong u,v,w;
	Ran(Ullong j) : v(4101842887655102017LL), w(1) {
		u = j ^ v; int64();
		v = u; int64();
		w = v; int64();
	}
	inline Ullong int64() {
		u = u * 2862933555777941757LL + 7046029254386353087LL;
		v ^= v >> 17; v ^= v << 31; v ^= v >> 8;
		w = 4294957665U*(w & 0xffffffff) + (w >> 32);
		Ullong x = u ^ (u << 21); x ^= x >> 35; x ^= x << 4;
		return (x + v) ^ w;
	}
	inline Doub doub() { return 5.42101086242752217E-20 * int64(); }
	inline Uint int32() { return (Uint)int64(); }
};

struct Normaldev_BM : Ran 
    {
	Doub mu,sig;
	Doub storedval;
	Normaldev_BM(Doub mmu, Doub ssig, Ullong i)
	: Ran(i), mu(mmu), sig(ssig), storedval(0.) {}
	Doub dev() 
    {
		Doub v1,v2,rsq,fac;
		if (storedval == 0.) 
        {
			do 
            {
				v1=2.0*doub()-1.0;
				v2=2.0*doub()-1.0;
				rsq=v1*v1+v2*v2;
			} 
            while (rsq >= 1.0 || rsq == 0.0);
			fac=sqrt(-2.0*log(rsq)/rsq);
			storedval = v1*fac;
			return mu + sig*v2*fac;
		} 
        else 
        {
			fac = storedval;
			storedval = 0.;
			return mu + sig*fac;
		}
	}
};

int main()
{
    Normaldev_BM norm_dev; 
    cout << norm_dev.Normaldev_BM(0,0.447,2);
    
    system("pause");
    return 0;
}

[MiiNiPaa]

You cannot call constructor manually.

Try:

Normaldev_BM norm_dev(0,0.447,2); 
std::cout << norm_dev.dev()

Doub What is this? Does it have overloaded << operator?

[johnny92]

Thank you! That worked like a charm, but now I seem to only be generating the same number over and over again. How do I go about changing this? So every time I print norm_dev.dev() i want a new normal deviate.


I think Doub is a typedef for double. I don't know why Numerical Recipes does this though.

[johnny92]

Every time I run my code I get the same 10 numbers.
I used this in my code:

int main()
{
    //Line below creates a new uniform number between 0 and 1 when I call ran.doub()
    Ran ran(time(0));    
    Normaldev_BM norm_dev(0,0.447,ran.doub()); 
    for(int i=0;i<10;i++)
    {
       cout << norm_dev.dev() << endl;;
    }
    system("pause");
    return 0;
}

I want to be able to generate 10 new numbers each time I run the code.

[johnny92]

I used the following code to solve my problem.

#include "nr3.h" 

struct Ran {
	Ullong u,v,w;
	Ran(Ullong j) : v(4101842887655102017LL), w(1) {
		u = j ^ v; int64();
		v = u; int64();
		w = v; int64();
	}
	inline Ullong int64() {
		u = u * 2862933555777941757LL + 7046029254386353087LL;
		v ^= v >> 17; v ^= v << 31; v ^= v >> 8;
		w = 4294957665U*(w & 0xffffffff) + (w >> 32);
		Ullong x = u ^ (u << 21); x ^= x >> 35; x ^= x << 4;
		return (x + v) ^ w;
	}
	inline Doub doub() { return 5.42101086242752217E-20 * int64(); }
	inline Uint int32() { return (Uint)int64(); }
};

struct Normaldev_BM : Ran 
    {
	Doub mu,sig;
	Doub storedval;
	Normaldev_BM(Doub mmu, Doub ssig, Ullong i)
	: Ran(i), mu(mmu), sig(ssig), storedval(0.) {}
	Doub dev() 
    {
		Doub v1,v2,rsq,fac;
		if (storedval == 0.) 
        {
			do 
            {
				v1=2.0*doub()-1.0;
				v2=2.0*doub()-1.0;
				rsq=v1*v1+v2*v2;
			} 
            while (rsq >= 1.0 || rsq == 0.0);
			fac=sqrt(-2.0*log(rsq)/rsq);
			storedval = v1*fac;
			return mu + sig*v2*fac;
		} 
        else 
        {
			fac = storedval;
			storedval = 0.;
			return mu + sig*fac;
		}
	}
};

int main()
{
    Normaldev_BM norm_dev(0,0.447,time(0)); 
    for(int i=0;i<10;i++)
    {
       cout << norm_dev.dev() << endl;;
    }
    system("pause");
    return 0;
}

[cire]

The third parameter for your Normaldev_BM constructor takes a value of type unsigned long long. Perhaps you should should feed it a value of that type instead of a double between 0 and 1?

int main()
{
    Ran ran(std::time(0)) ;
    Normaldev_BM norm_dev(0.0, 0.447, ran.int64()) ;

    for ( unsigned i=0; i<10; ++i )
        std::cout << norm_dev.dev() << '\n' ;
}