C++

Forum

 
x86 to x64 UINT_MAX and/or INT_MAX issue

ballzac (18)
I am working with code that was largely written written by someone else.

One of the functions adds poisson noise to an image. Basically, it takes the image (wf) and the noise level (nl), and the counts (count = 1.0/(nl*nl)), randomly picks a number from the poisson distribution with mean "count", and alters wf accordingly, also adjusting the value of "nl" based on the error between the original value of "wf" and the altered value.

I don't understand the details of how the poisson distribution is generated in this code, but I know exactly what the function should be doing.

It was working fine on 32 bit windows (I am using VC++2010 if it makes a difference), but I changed to 64 bit to make use of the extra memory, and am getting errors.

This is the relevant code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
 
void PoidevWF(vector<vector<double> > &wf, double &nl, const double count)
{
	const size_t ny = wf.size();
	const size_t nx = wf[0].size();

	//generate a negative random number use as seed to ran1 in poidev
	unsigned int number;
	rand_s(&number);
	int idum = -static_cast<int>((double)number*(double)INT_MAX/(double)UINT_MAX);
	double error = 0.0;
	double val;
	double sum = 0.0;
	for(size_t j=0; j<ny; j++)
	{
		for(size_t i=0; i<nx; i++)
		{
			val = wf[j][i];
			if(count > DBL_MIN) 
			{
				wf[j][i] = poidev(count*val,idum)/count;
			}
			error += (wf[j][i] -val)*(wf[j][i] -val);
			sum += val*val;
		}
	}
	nl = sqrt(error/sum);
}

double poidev(const double xm, int &idum)
{
	#define PI 3.141592653589793
	static double sq,alxm,g,oldm=(-1.0);
	double em,t,y;

	if (xm < 12.0) {
		if (xm != oldm) {
			oldm=xm;
			g=exp(-xm);
		}
		em = -1;
		t=1.0;
		do {
			em += 1.0;
			t *= ran1(idum);
		} while (t > g);
	} else {
		if (xm != oldm) {
			oldm=xm;
			sq=sqrt(2.0*xm);
			alxm=log(xm);
			g=xm*alxm-gammln(xm+1.0);
		}
		do {
			do {
				y=tan(PI*ran1(idum));
				em=sq*y+xm;
			} while (em < 0.0);
			em=floor(em);
			t=0.9*(1.0+y*y)*exp(em*alxm-gammln(em+1.0)-g);
		} while (ran1(idum) > t);
	}
	return em;
	#undef PI
}

double ran1(int &idum)
{
	#define M1 259200
	#define IA1 7141
	#define IC1 54773
	#define RM1 (1.0/M1)
	#define M2 134456
	#define IA2 8121
	#define IC2 28411
	#define RM2 (1.0/M2)
	#define M3 243000
	#define IA3 4561
	#define IC3 51349
	//
	static long ix1,ix2,ix3;
	static double r[98];
	double temp;
	static int iff=0;
	int j;
	void nrerror();

	if (idum < 0 || iff == 0) {
		iff=1;
		ix1=(IC1-(idum)) % M1;
		ix1=(IA1*ix1+IC1) % M1;
		ix2=ix1 % M2;
		ix1=(IA1*ix1+IC1) % M1;
		ix3=ix1 % M3;
		for (j=1;j<=97;j++) {
			ix1=(IA1*ix1+IC1) % M1;
			ix2=(IA2*ix2+IC2) % M2;
			r[j]=(ix1+ix2*RM2)*RM1;
		}
		idum=1;
	}
	ix1=(IA1*ix1+IC1) % M1;
	ix2=(IA2*ix2+IC2) % M2;
	ix3=(IA3*ix3+IC3) % M3;
	j=1 + (97*ix3)/M3;
	/*extern vector<int> jCount;
	jCount[j-1]=jCount[j-1]+1;
	*/
	//
	if (j > 97 || j < 1) 
	{
		if(j==0 || j==-97)
		cerr << "RAN1: This cannot happen." << endl;
	}
	temp=r[j];
	r[j]= (ix1+ix2*RM2)*RM1;
	return temp;

	#undef M1
	#undef IA1
	#undef IC1
	#undef RM1
	#undef M2
	#undef IA2
	#undef IC2
	#undef RM2
	#undef M3
	#undef IA3
	#undef IC3
}

double gammln(double xx)
{
	double x,tmp,ser;
	static double cof[6]={76.18009173,-86.50532033,24.01409822,
		-1.231739516,0.120858003e-2,-0.536382e-5};
	int j;

	x=xx-1.0;
	tmp=x+5.5;
	tmp -= (x+0.5)*log(tmp);
	ser=1.0;
	for (j=0;j<=5;j++) {
		x += 1.0;
		ser += cof[j]/x;
	}
	return -tmp+log(2.50662827465*ser);
}



The problem is that when I switched to 64 bit, I started getting the output
RAN1: This cannot happen.

Looking through the code, I saw that this occurs here:
1
2
3
4
5
 
if (j > 97 || j < 1) 
	{
		if(j==0 || j==-97)
		cerr << "RAN1: This cannot happen." << endl;
	}


It turns out that occasionally I get a negative number for j. This happens only very occasionally, but when it does, it happens dozens of times in a row. The negative values for j were greater than -96, so I changed this:
j=1 + (97*ix3)/M3;
to this
j=1 + abs((97*ix3)/M3);
and the code ran without error. But then when I looked at my results, it appears that those values for j are producing erroneous results, rarely, but randomly. Obviously my solution did not work.

The only thing I can put the negative j values down to is that UINT_MAX and INT_MAX are used in the creation of the seed. I have tried searching around to find how these numbers are different on different operating systems, but have not had much luck.

Ultimately, I might need to find another function for the poisson distribution somewhere, but at the moment I would prefer to leave the code mostly as it is and just change something about those values if need be. If anyone has any idea how to solve this problem I would greatly appreciate the advice :) Cheers
Peter87 (11234)
Could the problem be that sizeof(long) is 4 on x86 and 8 on x64?

I don't know if that's the case or if that would make a difference here but you could try change long to int and see if that fixes it.
ballzac (18)
I hadn't thought of that. I will try it, thanks. One of the problems though is that, because it happens randomly, the code can run for several hours before the problem occurs. This makes it very difficult to test.
ballzac (18)
Unfortunately that didn't fix it. :(


I just realised I made a mistake in copying the code here. When I removed some lines that I had put in to help me troubleshoot the problem, I accidentally left one.

This
1
2
3
4
5
 
if (j > 97 || j < 1) 
	{
		if(j==0 || j==-97)
		cerr << "RAN1: This cannot happen." << endl;
	}

should actually read
1
2
3
4
 
if (j > 97 || j < 1) 
	{
			cerr << "RAN1: This cannot happen." << endl;
	}


Any other ideas as to what could be done to fix the problem?
JLBorges (13770)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
 
#include <random>

// randomly picks a number from the poisson distribution with mean "mean"
int rand_poisson( double mean )
{
    // if there is a random_device
    // static std::random_device rdev() ;
    // static std::mt19937 rng( rdev() ) ;

    // else
    static std::mt19937 rng( std::time(nullptr) ) ;

    std::poisson_distribution<int> distribution(mean) ;
    return distribution(rng) ;
}
ballzac (18)
Thanks for that. It seems to be working fine. Just a quick question though. I get this:

 warning C4244: 'argument' : conversion from 'time_t' to 'unsigned long', possible loss of data


is that going to cause any problems?
JLBorges (13770)
> is that going to cause any problems?

In this particular instance, no - you just want the low order bits of current time to seed the random number generator.

It appears that on your implementation, std::time_t is a 64 bit value; this allows representation of time values beyond 2038 AD. If you used a std::time_t to actually represent a point in time, this truncation to 32 bits would be significant for periods after 2038 AD.
ballzac (18)
Cool. Thanks heaps for your help :)
Peter87 (11234)
To get rid of that warning you can use a cast
static std::mt19937 rng(static_cast<std::uint_fast32_t>(std::time(nullptr)));
Topic archived. No new replies allowed.