001package cnslab.cnsmath;
002public class Cnsran
003{
004        public  static int idum2=123456789,iy=0;
005        public  static int[] iv=new int[32];
006        private static int iset=0;
007        private static double gset;
008
009        public static double ran2(Seed idum)
010        {       
011                int IM1=2147483563,IM2=2147483399;
012                int IA1=40014,IA2=40692,IQ1=53668,IQ2=52774;
013                int IR1=12211,IR2=3791,NTAB=32,IMM1=IM1-1;
014                int NDIV=1+IMM1/NTAB;
015                double EPS=3.0e-16,RNMX=1.0-EPS,AM=1.0/(double)(IM1);
016                int j,k;
017                double temp;
018
019                if (idum.seed <= 0) {
020                        idum.seed=(idum.seed==0 ? 1 : -idum.seed);
021                        idum2=idum.seed;
022                        for (j=NTAB+7;j>=0;j--) {
023                                k=idum.seed/IQ1;
024                                idum.seed=IA1*(idum.seed-k*IQ1)-k*IR1;
025                                if (idum.seed < 0) idum.seed += IM1;
026                                if (j < NTAB) iv[j] = idum.seed;
027                        }
028                        iy=iv[0];
029                }
030                k=idum.seed/IQ1;
031                idum.seed=IA1*(idum.seed-k*IQ1)-k*IR1;
032                if (idum.seed < 0) idum.seed += IM1;
033                k=idum2/IQ2;
034                idum2=IA2*(idum2-k*IQ2)-k*IR2;
035                if (idum2 < 0) idum2 += IM2;
036                j=iy/NDIV;
037                iy=iv[j]-idum2;
038                iv[j] = idum.seed;
039                if (iy < 1) iy += IMM1;
040                if ((temp=AM*iy) > RNMX) return RNMX;
041                else return temp;
042        }
043
044
045        public static double ran1(Seed idum)
046        {       
047                int IA=16807,IM=2147483647,IQ=127773,IR=2836,NTAB=32;
048                int NDIV=(1+(IM-1)/NTAB);
049                double EPS=3.0e-16,AM=1.0/IM,RNMX=(1.0-EPS);
050                int j,k;
051                double temp;
052
053                if (idum.seed <= 0 || iy==0) {
054                        if (-idum.seed < 1) idum.seed=1;
055                        else idum.seed = -idum.seed;
056                        for (j=NTAB+7;j>=0;j--) {
057                                k=idum.seed/IQ;
058                                idum.seed=IA*(idum.seed-k*IQ)-IR*k;
059                                if (idum.seed < 0) idum.seed += IM;
060                                if (j < NTAB) iv[j] = idum.seed;
061                        }
062                        iy=iv[0];
063                }
064                k=idum.seed/IQ;
065                idum.seed=IA*(idum.seed-k*IQ)-IR*k;
066                if (idum.seed < 0) idum.seed += IM;
067                j=iy/NDIV;
068                iy=iv[j];
069                iv[j] = idum.seed;
070                if ((temp=AM*iy) > RNMX) return RNMX;
071                else return temp;
072        }
073
074        public static double gasdev(Seed idum)
075        {       
076                double fac,rsq,v1,v2;
077
078                if (idum.seed < 0) iset=0;
079                if (iset == 0) {
080                        do {    
081                                v1=2.0*ran1(idum)-1.0;
082                                v2=2.0*ran1(idum)-1.0;
083                                rsq=v1*v1+v2*v2;
084                        } while (rsq >= 1.0 || rsq == 0.0);
085                        fac=Math.sqrt(-2.0*Math.log(rsq)/rsq);
086                        gset=v1*fac;
087                        iset=1;
088                        return v2*fac;
089                } else {
090                        iset=0;
091                        return gset;
092                }
093        }
094}