Gillespie algorithm -- where's the glitch?

[JaroThe1st]

Hello All,

I am very new to c++ and I've been trying to write a simple stochastic algorithm (aka Gillespie algorithm). I know what the answer should be, but c++ is giving something not quite right; instead of 100 and 1000, I'm getting 96 and 871. For the life of me I can't find a problem. If any of you could help me out here, that would be much appreciated. Again, I'm a beginner, so please don't laugh at the code. Here it is:

#include <iostream>
#include <cmath>
#include <cstdlib>
#include <ctime>

using namespace std;

int main()
{

srand(time(NULL));

double z1;
double z2;
double r=10.0;
double k=0.1;
double K=1.0;
double q=0.1;
double xx=0;
double yy=0;
double x=0;
double y=0;
double t;
double r1;
double r2;
double r3;
double r4;
double p1;
double p2;
double p3;
double p4;
double a0;
double tau;

int zz1;
int zz2;
int M=1000;
int m=0;

while(m<M){

x=0;
y=0;
t=0.0;

while(t<=500){

zz1 = rand()%32767+1;
zz2 = rand();

z1=zz1/32767.0;
z2=zz2/32767.0;

r1=r;
r2=k*x;
r3=K*x;
r4=q*y;

p1=r1;
p2=p1+r2;
p3=p2+r3;
p4=p3+r4;

a0=z2*p4;
tau = log(1.0/z1)/a0;

if(p1>a0){
x = x + 1;
}

if(p1<=a0&&p2>a0){
x = x - 1;
}

if(p2<=a0&&p3>a0){
y = y + 1;
}

if(p3<=a0&&p4>a0){
y = y - 1;
}

t = t + tau;

}

xx = xx+x;
yy = yy+y;

m++;
}

cout<< xx/M << endl;
cout<< yy/M << endl;

}

Thanks a lot!

[chicofeo]

PLEASE USE CODE TAGS (the <> formatting button) when posting code.
It makes it easier to read your code and also easier to respond to your
post.
http://www.cplusplus.com/articles/jEywvCM9/
Hint: You can edit your post, highlight your code and press the <>
formatting button.