#include <tc.h>

#include <math.h>
#include <stdio.h>
#include <stdlib.h>

typedef unsigned long long u64;
typedef long long s64;
typedef double d64;

typedef struct {
  u64 rec_count;
  u64 start;
  u64 size;
} pack_prm;

typedef struct {
  u64 sum[10];
  d64 Sx;
  d64 Sy;
} pack_res;


#define CONST_1 ((u64)1<<46)

#define max(a,b)  ((a) > (b) ? (a) : (b))

#ifdef USE_MACROS
#define lfabs(a)  ((a) < 0 ? (a)*(d64)(-1) : a)
#define get_ab_mod_2_46(a,b)  ((((a) % CONST_1)*((b) % CONST_1)) % CONST_1)

#else

#ifndef TSYSTEM
d64 max(d64 a, d64 b){
  return a > b ? a : b;
}
#endif

d64 lfabs(d64 a){
  return a < 0 ? -a : a ;
}

u64 get_ab_mod_2_46(u64  a, u64  b ){
  return ((a % CONST_1)*(b % CONST_1)) % CONST_1;
}

#endif

u64 fast_raise(u64 a, u64 x){
  u64 res=1;
  while( x > 0 ){
    if ( x & (u64)1) res=get_ab_mod_2_46(res,a);
    a=get_ab_mod_2_46(a,a);
    x >>= 1;
  }
  return res;
}

void make_gist(u64 k,u64 packet_len, u64 a, u64 s, pack_res *res){
  u64 count=0;
  u64 tmp1,tmp2,x0=fast_raise(a,(u64)2*k);
  tmp2=get_ab_mod_2_46(x0,s); //x_0
  for(; count < (u64)packet_len; count++){
    d64 x,y,t;
    tmp1=get_ab_mod_2_46(a,tmp2);          //x_1
    tmp2=get_ab_mod_2_46(a,tmp1);          //y_1
    x=(d64)((s64)(tmp1<<1) - (s64)CONST_1)/(d64)CONST_1;
    y=(d64)((s64)(tmp2<<1) - (s64)CONST_1)/(d64)CONST_1;
    t=x*x+y*y;
    if(t<=1.0){
      t=sqrt(-2.0*log(t)/t);
      x*=t;
      y*=t;
      res->sum[(int)max(lfabs(x),lfabs(y))]++;
      res->Sx+=x;
      res->Sy+=y;
    }
  }
}

tfun void iterator (pack_prm const parm, tout pack_res *res)
{
  int i;
  if(parm.rec_count == 0 ){
    memset((void*)res,0,sizeof(*res));
    make_gist(parm.start,parm.size,(u64)pow(5,13),271828183,res);
  } else {
    pack_prm pt1,pt2;
    tval pack_res res1,res2;
    pt1.rec_count=pt2.rec_count=parm.rec_count-1;
    pt1.start=parm.start;
    pt2.start=parm.start+parm.size/2;
    pt1.size=pt2.size=parm.size/2;

    iterator(pt1,&res1);
    iterator(pt2,&res2);

    for(i=0;i<10;i++) res->sum[i] = res1.sum[i] + res2.sum[i];
    res->Sx = res1.Sx + res2.Sx;
    res->Sy = res1.Sy + res2.Sy;
  }
}

int tmain (int argc, char* argv[])
{
  int i;
  tval pack_res res;
  pack_prm prm = { (u64)10, (u64)0, (u64)1<<28 };

  terrprint("v20 sizeof(u64)=%d sizeof(d64)=%d\nProblem size 0x%Lx\n",
         sizeof(u64),sizeof(d64),prm.size);

  iterator(prm,&res);

  for(i=0;i<10;i++) {
    terrprint("Q[%d]=%Lu\n",i,res.sum[i]);
  }
  terrprint("Sx=%le\nSy=%le\n",res.Sx,res.Sy);
  return 0;
}