Нормальное Распределение чисел - C++ - Ответ 1752557

#include <math.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 {
  d64 Sx;
  d64 Sy;
} pack_res;
 
#define CONST_1 ((u64)1<<46)
 
#ifndef max
#define max(a,b)  ((a) > (b) ? (a) : (b))
#endif
 
#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
 
#if 0
#ifndef __tplusplus
d64 max(d64 a, d64 b){
  return a > b ? a : b;
}
#endif
#endif
 
inline d64 lfabs(d64 a) {
    return a < 0 ? -a : a ;
}
 
inline u64 get_ab_mod_2_46(u64 a, u64 b) {
    return ((a % CONST_1)*(b % CONST_1)) % CONST_1;
}
 
#endif
 
inline 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->Sx=x;
      res->Sy=y;
    }
  }
}
 
tfun pack_res iterator (pack_prm parm) {
  int i;
  pack_res res;
  if(parm.rec_count == 0 )
  {
    
    res.Sx = res.Sy = 0;
    make_gist(parm.start,parm.size,(u64)pow(5.0,13),271828183,&res);
  } 
 else 
  {
    pack_prm pt1,pt2;
    tval pack_res res1,res2;
    pt1.rec_count=pt2.rec_count=parm.rec_count-1;
    pt2.start=(pt1.start=parm.start)+(pt1.size=pt2.size=parm.size/2);
 
    res1 = iterator(pt1);
    res2 = iterator(pt2);
 
  
    res.Sx = ((pack_res&)res1).Sx + ((pack_res&)res2).Sx;
    res.Sy = ((pack_res&)res1).Sy + ((pack_res&)res2).Sy;
  }
  return res;
}
 
tfun int main (int argc, char* argv[]) {
  int i=1;
  pack_prm prm = { (u64)10, (u64)0, (u64)1 << 28 };
                      //^^-depth              ^^-size
  while (i < argc) {
    if (argv[i][0] == '+') argv[i][0] = '-';
    if (strcmp(argv[i], "-size") == 0) {
      prm.size = 1ULL << atol(argv[++i]);
    } else if (strcmp(argv[i], "-depth") == 0) {
      prm.rec_count = atol(argv[++i]);
    } else {
      printf(
        "Usage:\n"
        "ep [-size <problem_size>] [-depth <tree depth>]\n"
      );
      return -1;
    }
    i++;
  }
 
  printf("v20 sizeof(u64)=%lu sizeof(d64)=%lu\nProblem size %Lu\n",
         sizeof(u64),sizeof(d64),prm.size);
 
  pack_res res = iterator(prm);
 
  printf("Sx=%le\nSy=%le\n", res.Sx, res.Sy);
  return 0;
}