/************************************************
 Expect bugs!
 Please use and enjoy, and let me know of any bugs/mods/improvements 
 that you have found/implemented and I will fix/incorporate them into 
 this file. Thank Mr. Xushiliang once again.

					hujinshan@2002.11.3
				Airforce Engineering University
************************************************/

/*****  #include "CH13.h"  隨機數的產生*****/
#ifndef CH13_H_
#define CH13_H_

#include "stdlib.h"
#include "math.h"
#include "stdio.h"
//*******************************************************************

double mrnd1(double* r);//0-1之間均勻分布的一個隨機數
void mrnds(double r,double p[],int n);//0-1之間均勻分布的隨機數序列
int mrab1(int a,int b,int* r);//任意區間內均勻分布的一個隨機整數
void mrabs(int a,int b,int* r,int p[],int n);//任意區間內均勻分布的隨機整數序列
double mgrn1(double u,double g,double* r);//任意均值與方差的一個正態分布隨機數
void mgrns(double u,double g,double* r,int n,double a[]);//任意均值與方差的正態分布隨機數序列

//*******************************************************************
double mrnd1(double* r)
{ 
	int m;
    double s,u,v,p;
    s=65536.0; u=2053.0; v=13849.0;
    m=(int)(*r/s); *r=*r-m*s;
    *r=u*(*r)+v; m=(int)(*r/s);
    *r=*r-m*s; p=*r/s;
    return(p);
}
/////////////////////////////////////////////////////////////
void mrnds(double* r,double p[],int n)
{
	int i,m;
    double s,u,v;
    s=65536.0; u=2053.0; v=13849.0;
    for (i=0; i<=n-1; i++)
      { *r=u*(*r)+v; m=(int)(*r/s);
        *r=*r-m*s; p[i]=*r/s;
      }
    return;
}
/////////////////////////////////////////////////////////////
int mrab1(int a,int b,int* r)
{ 
	int k,l,m,i,p;
    k=b-a+1; l=2;
    while (l<k) l=l+l;
    m=4*l; k=*r; i=1;
    while (i<=1)
      { k=k+k+k+k+k;
	k=k%m; l=k/4+a;
        if (l<=b) { p=l; i=i+1;}
      }
    *r=k;
    return(p);
}
/////////////////////////////////////////////////////////////
void mrabs(int a,int b,int* r,int p[],int n)
{ 
	int k,l,m,i;
    k=b-a+1; l=2;
    while (l<k) l=l+l;
    m=4*l; k=*r; i=0;
    while (i<=n-1)
      { k=k+k+k+k+k;
	k=k%m; l=k/4+a;
        if (l<=b) { p[i]=l; i=i+1;}
      }
    *r=k;
    return;
}
/////////////////////////////////////////////////////////////
double mgrn1(double u,double g,double* r)
{ 
	int i,m;
    double s,w,v,t;
    s=65536.0; w=2053.0; v=13849.0;
    t=0.0;
    for (i=1; i<=12; i++)
      { *r=(*r)*w+v; m=(int)(*r/s);
        *r=*r-m*s; t=t+(*r)/s;
      }
    t=u+g*(t-6.0);
    return(t);
}
/////////////////////////////////////////////////////////////
void mgrns(double u,double g,double* r,int n,double a[])
{ 
	int i,k,m;
    double s,w,v,t;
    s=65536.0; w=2053.0; v=13849.0;
    for (k=0; k<=n-1; k++)
      { t=0.0;
        for (i=1; i<=12; i++)
          { *r=(*r)*w+v; m=(int)(*r/s);
            *r=*r-m*s; t=t+(*r)/s;
          }
        a[k]=u+g*(t-6.0);
      }
    return;
}
/////////////////////////////////////////////////////////////
#endif