float xrnd(unsigned int idx){
  {
    idx%=s.rsize;
    s.ml=rm[idx];
    s.da=rs[idx];
  }
  union{
    unsigned int tmp;
    float xrd;
  };
  do{
    s.da = s.ml * (unsigned long long) s.ax + s.dx;
    tmp = s.ax >> 9;
  } while(tmp==0);
  tmp |= 0x3f800000;
  {
    rs[idx]=s.da;
  }
  return xrd-1;
}

float xrnd(state & s){
  union{
    unsigned int tmp;
    float xrd;
  };
  do{
    s.da = s.ml * (unsigned long long) s.ax + s.dx;
    tmp = s.ax >> 9;
  } while(tmp==0);
  tmp |= 0x3f800000;
  return xrd-1;
}

void rotate(float cs, float si, float *n, state & s){
  float u[3];
  {
    float p1[3], p2[3];
    {
      int i0=0;
      float n0=n[0];
      for(int i=1; i<3; i++) if(fabsf(n[i])>fabsf(n0)){ i0=i; n0=n[i]; }
      int i1=(i0+1)%3, i2=(i0+2)%3;
      float n1=n[i1], n2=n[i2], nq=n0*n0;
      float r1=1/sqrtf(nq+n1*n1), r2=1/sqrtf(nq+n2*n2);

      p1[i0]=-n1*r1; p1[i1]=n0*r1; p1[i2]=0;
      p2[i0]=-n2*r2; p2[i1]=0; p2[i2]=n0*r2;
    }
    {
      float q1[3], q2[3], r1=0, r2=0;
      for(int i=0; i<3; i++){
	float a1=p1[i]-p2[i]; q1[i]=a1; r1+=a1*a1;
	float a2=p1[i]+p2[i]; q2[i]=a2; r2+=a2*a2;
      }

      r1=1/sqrtf(r1); r2=1/sqrtf(r2);

      for(int i=0; i<3; i++){
	p1[i]=q1[i]*r1;
	p2[i]=q2[i]*r2;
      }
    }

    {
      float zi=xrnd(s); zi*=2*M_PI;
      float px=cosf(zi), py=sinf(zi);
      for(int i=0; i<3; i++) u[i]=px*p1[i]+py*p2[i];
    }
  }

  {
    float r=0;
    for(int i=0; i<3; i++){
      float a=n[i]*cs+u[i]*si;
      n[i]=a; r+=a*a;
    }
    {
      r=1/sqrtf(r);
      for(int i=0; i<3; i++) n[i]*=r;
    }
  }
}

void propagate(unsigned int idx){
  {
    idx%=s.rsize;
    s.ml=rm[idx];
    s.da=rs[idx];
  }

  { // initialize photon
    t=xt;
    for(int i=0; i<3; i++){
      r[i]=xr[i]; n[i]=xn[i];
    }

    if(type==1) rotate(w->coschr, w->sinchr, n, s);
    else{
      static const float rms=9.2*sqrtf(3)*M_PI/180;
      float xi=xrnd(s);
      xi=rms*(2*xi-1);
      float nz=sinf(xi), np=cosf(xi);

      xi=xrnd(s);
      xi*=2*M_PI;
      n[0]=np*cosf(xi);
      n[1]=np*sinf(xi);
      n[2]=nz;
    }
  }

  float TOT=-logf(xrnd(s));
  while(true){
    float SCA=-logf(xrnd(s));

    float sca;
    { // get distance for overburden x
      float & x = SCA;
      float & z = r[2];
      float & nz = n[2];

      int i=lroundf((z-hmin)*rdh), j;
      if(i<0) i=0; else if(i>=size) i=size-1;

      float h=hmin+i*dh;
      float ai, ah;

      if(nz>0){
	ah=h+hdh;
	ai=(nz*x-(ah-z)*w->z[i].sca)*rdh;
	for(j=i; j<size-1 && ai>0; ah+=dh) ai-=w->z[++j].sca;
      }
      else{
	ah=h-hdh;
	ai=(nz*x-(ah-z)*w->z[i].sca)*rdh;
	for(j=i; j>0 && ai<0; ah-=dh) ai+=w->z[--j].sca;
      }

      if(i==j) sca=x/w->z[i].sca;
      else sca=(ai*dh/w->z[j].sca+ah-z)/nz;
    }

    float tot;
    { // get distance for overburden x
      float & x = TOT;
      float & z = r[2];
      float & nz = n[2];

      int i=lroundf((z-hmin)*rdh), j;
      if(i<0) i=0; else if(i>=size) i=size-1;

      float h=hmin+i*dh;
      float ai, ah;

      if(nz>0){
	ah=h+hdh;
	ai=(nz*x-(ah-z)*w->z[i].abs)*rdh;
	for(j=i; j<size-1 && ai>0; ah+=dh) ai-=w->z[++j].abs;
      }
      else{
	ah=h-hdh;
	ai=(nz*x-(ah-z)*w->z[i].abs)*rdh;
	for(j=i; j>0 && ai<0; ah-=dh) ai+=w->z[--j].abs;
      }

      if(i==j) tot=x/w->z[i].abs;
      else tot=(ai*dh/w->z[j].abs+ah-z)/nz;
    }

    if(tot<sca) sca=tot;

    name om;
    { // sphere

#ifdef ACCL0
      float dx[3], dn[3];
      for(int m=0; m<3; m++){
	float ini=r[m];
	float fin=ini+sca*n[m];
	if(n[m]>0){
	  dn[m]=ini;
	  dx[m]=fin;
	}
	else{
	  dn[m]=fin;
	  dx[m]=ini;
	}
	dn[m]-=R;
	dx[m]+=R;
      }
#endif

#ifdef ACCL2
      int xl[3], xh[3];
      for(int m=0; m<3; m++){
	xl[m]=min(max((int) lroundf((dn[m]-cl[m])*crst[m]), 0), cn[m]);
	xh[m]=max(min((int) lroundf((dx[m]-cl[m])*crst[m]), cn[m]-1), -1);
      }
      for(int i=xl[0]; i<=xh[0]; i++) for(int j=xl[1]; j<=xh[1]; j++) for(int k=xl[2]; k<=xh[2]; k++) for(int l=cidx[i][j][k]; carr[l]!=-1; l++){
	const DOM & dom=oms[carr[l]];
#ifdef ZZZ
      }
#endif
#else
      for(int l=0; l<gsize; l++){
	const DOM & dom=oms[l];
#endif
#ifdef ACCL1
	bool pass=false;
	for(int i=0; i<3; i++) if(dom.r[i]<dn[i] || dom.r[i]>dx[i]){ pass=true; break; }
	if(pass) continue;
#endif

	{
	  float b=0, c=0;
	  for(int i=0; i<3; i++){
	    float dr=dom.r[i]-r[i];
	    b+=n[i]*dr; c+=dr*dr;
	  }
	  c-=R*R;
	  float D=b*b-c, h=-1;
	  if(D>0) h=b-sqrtf(D);

	  if(h>0 && h<sca){ om=dom.n; sca=h; }
	}
      }
    }

    float xa;
    { // get overburden for distance a
      float & a = sca;
      float & z = r[2];
      float & nz = n[2];

      int i=lroundf((z-hmin)*rdh), j;
      if(i<0) i=0; else if(i>=size) i=size-1;

      float h=hmin+i*dh;
      float y=z+nz*a;

      j=lroundf((y-hmin)*rdh);
      if(j<0) j=0; else if(j>=size) j=size-1;

      if(i==j) xa=a*w->z[i].abs; 
      else{
	float g=hmin+j*dh;
	if(nz>0){
	  xa=(y-(g-hdh))*w->z[j].abs;
	  while(--j>i) xa+=hdh*w->z[j].abs;
	  xa+=(h+hdh-z)*w->z[i].abs;
	}
	else{
	  xa=(y-(g+hdh))*w->z[j].abs;
	  while(++j<i) xa-=hdh*w->z[j].abs;
	  xa+=(h-hdh-z)*w->z[i].abs;
	}
	xa/=nz;
      }
    }

    TOT-=xa;

    { // advance
      float & a= sca;
      t+=a*w->ocm;
      for(int i=0; i<3; i++) r[i]+=n[i]*a;
    }

    if(om.k!=0){
      hit h; h.n=om; h.t=t;

      float sum;
      { // hole ice
	const float s[11]={0.32752, 0.66696, 0.35389, -1.4107, -1.7738, 4.3302, 5.7151, -5.5262, -7.9000, 2.0155, 3.3965};
	float & x= n[2];
	float y=1;
	sum=s[0];
	for(int i=1; i<11; i++){ y*=x; sum+=s[i]*y; }
      }

      if(xrnd(s)<sum) if(hidx<hnum) hits[hidx++]=h;
      break;
    }

    if(TOT<xx) break;

    {
      float xi=xrnd(s);
      xi=2*xi-1;
      if(g!=0){
	float g2=g*g;
	float ga=(1-g2)/(1+g*xi);
	xi=(1+g2-ga*ga)/(2*g);
      }
      if(xi>1) xi=1; else if(xi<-1) xi=-1;

      float si=sqrtf(1-xi*xi);
      rotate(xi, si, n, s);
    }
  }

  {
    rs[idx]=s.da;
  }

}