using System; 
public class RRan { // even distribution from 0 to sizetf-1 !!! 
    const int NTAB = 32;  
    const ulong IA1=2521, IM1=4503599627370449, capa1=IM1-1,
                    IA2=2503, IM2=4503599627370353, capa2=IM2-1, 
                    NDIV=(1+(capa1/NTAB))/123456; // === nw 
    ulong[] idums = new ulong[NTAB];
    ulong _idum, _idum2; // iy is usually 0 .. capa1-1 
    ulong iy = 0, viewlimit, size; 
        // ============= construction ==================  
    public RRan(ulong sizetf, ulong initial) { // for size: 15 zeros possible !!! 
        ulong ulm = ulong.MaxValue; 
        if (size > capa1) {
            System.Windows.Forms.MessageBox.Show("random size too large !!! "); }
        size = sizetf; 
        ulong folds = IM1 / size;
        viewlimit = folds * size; 
        // 
        if (initial == 0) {_idum=IA1; } else {_idum=initial; }  
        _idum2=_idum; 
        for (var j=0; j < 10; j++) { 
            kstep(ref _idum, IA1, IM1);   
            kstep(ref _idum2, IA2, IM2); } 
        for (var j=0; j < NTAB; j++) { 
            kstep(ref _idum, IA1, IM1);   
            idums[j]=_idum; } 
    } // ======================================================= 
    void kstep(ref ulong _idum, ulong IA, ulong IM) { // 1 .. IM-1 will be the output !!! 
        _idum *= IA; 
        _idum %= IM; 
    } // ======================================================= 
    public ulong step() { // === 0 .. size-1 will be answered !!! 
        kstep(ref _idum, IA1, IM1); // idum: 1..capa1 
        kstep(ref _idum2, IA2, IM2); // idum2: 1..capa2 
        ulong j1 = iy / NDIV; // j is 0..32*123456 - 1
        ulong j = j1%32; // j is 0..31 
        if (idums[j] >= _idum2) {
            iy = idums[j] - _idum2; 
        } else {
            ulong aaa = capa1 - _idum2; // === step is required here !!! 
            iy = aaa + idums[j]; } 
        idums[j]=_idum; 
        if (iy<viewlimit) {
            ulong toret = iy % size;
            return toret; 
        } else {
            return step(); } 
    } // ======================================================= 
} // ############################################################# 
public class DRan01 { // === flat distribution 0 to 1 !!!
    RRan rran = new RRan(1000000000000000, 253); // === 10e15 !!! 
    public double step() {
        double toret=1e-15*rran.step();
        return toret; 
    } // ========================================================= 
} // ############################################################# 
public class DRanExp { // === exponential distribution of sigma=1 !!! 
    DRan01 ran=new DRan01(); 
    public double step() { 
        double rn=ran.step();   
        double R = -Math.Log(rn); 
        return R;  
    } // ========================================================= 
} // ############################################################# 
public class DRanNormal { // === normal distribution of sigma=1 !!! 
    DRan01 ran=new DRan01(); 
    public double step() { 
        double rn=ran.step();  
        double R = -Math.Log(rn); 
            // 半径２乗に関する累積分布関数である指数関数（半径２乗　→　頻度）の逆関数である
            // 対数関数（頻度　→　半径２乗）を利用して、半径２乗の分布を生成！！！
        double r = Math.Sqrt(2 * R);   
            // 半径を復元！！！
        double tt = Math.PI * 2 * ran.step();     
            // 角度方向のランダムな分布を生成、極座標的にプロットし、
        return r*Math.Sin(tt); 
            // 一軸の成分を抽出！！！
    } // ========================================================= 
} // #############################################################