参阅我的文章：
成都创新互联公司长期为成百上千家客户提供的网站建设服务，团队从业经验10年，关注不同地域、不同群体，并针对不同对象提供差异化的产品和服务；打造开放共赢平台，与合作伙伴共同营造健康的互联网生态环境。为富源企业提供专业的网站设计、网站建设，富源网站改版等技术服务。拥有10多年丰富建站经验和众多成功案例,为您定制开发。
#include "stdafx.h" //VS 预编译头文件，其他系统请删除
#includestdio.h
#includestdlib.h
#includememory.h
#includemath.h
#includetime.h
//VS 2013 否决了 scanf 等函数，为了使用，加上下句。
//其他系统请删除
#pragma warning(disable:4996)
int GaussJordanElimination(int n, const double *pCoef, double *pOut);
//VS 主函数签名格式。其他系统请改变签名，如：
//int main()
int _tmain(int argc, _TCHAR* argv[])
{
double cf[3][4] = { {-0.02, 2.0, 2.0, 0.4}, {1.0, 0.78125, 0.0, 1.3816}, {3.996, 5.526, 4.0, 7.4178} };
double rs[3];
int i;
i = GaussJordanElimination(3, (double*)cf, rs);
printf("x1 = %lf, x2 = %lf, x3 = %lf\n", rs[0], rs[1], rs[2]);
system("pause"); //避免窗口一闪而退
return 0;
}
//绝对值函数
__inline double _abs(double v)
{
return v 0 ? -v : v;
}
//线性方程组列主元高斯消元法
//n 方程元数；pCoef 系数，必须以行主序方式存放的二维数组；
//pOut 长度为 n 的一维数组(调用者负责维护)，用于输出数据
//返回值：0 成功，-1 无解，1 申请内存失败, 2 不定解。
int GaussJordanElimination(int n, const double *pCoef, double *pOut)
{
double *pcf;
int rows = n, columns = n + 1;
//pcf = new double[rows * columns];
pcf = (double*)malloc(rows * columns * sizeof(double));
if (pcf == 0) return 1; //巧妇难为无米之炊，内存都申请不到，还能干嘛！
memcpy(pcf, pCoef, (rows * columns) * sizeof(double)); //据说这个运行效率很高
int r, c, i; //循环变量
int a, b;
double x, y;
//开始消元，将 pcf 方阵区处理成到直角三角形（直角在右上角）矩阵
for (r = 0; r rows - 1; r++)
{
//选取主元
a = r; x = _abs(pcf[r * columns + r]);
for (i = r + 1; i rows; i++)
{ //查找主元在哪行
if (x _abs(pcf[i * columns + r])) a = i;
}
if (a r)
{ //主元不是当前行（r）,比较麻烦，需要将第 a 行与第 r 行兑换
//第 r 列前面的就不要对换了，因为这些项已经被消元，变成 0 了
for (c = r; c columns; c++)
{
x = pcf[r * columns + c];
pcf[r * columns + c] = pcf[a * columns + c];
pcf[a * columns + c] = x;
}
}
//开始消元
a = r * columns; //记住将主元的行地址偏移量，以提高程序运行效率
x = -pcf[a + r]; //要多次使用，记下她，以提高程序运行效率
if (x == 0) //主元居然为 0，纯粹是想坑爹，岂能上当！
continue; //继续后面的消元，以便最终判断是无解还是任意解
for (i = r + 1; i rows; i++)
{ //正在消元
b = i * columns;//记住将要消元的行地址偏移量，以提高程序运行效率
y = pcf[b + r]; //要多次使用，记下她，以提高程序运行效率
if (y != 0)
{ //y == 0，本行不需要消元
y /= x; //要多次使用，记下她，以提高程序运行效率
pcf[b + r] = 0; //肯定为 0，不用计算。
for (c = r + 1; c columns; c++)
pcf[b + c] += pcf[a + c] * y;
}
}
}//至此，pcf 方阵区已经处理成到直角三角形（直角在右上角）矩阵
//回代，将 pcf 方阵区处理成主对角线为 1，其他为 0 的矩阵
int columns_1 = c = columns - 1; //多次用到，提高效率
for (r = rows - 1; r = 1; r--)
{
b = r * columns;
if (pcf[b + r] == 0)
{ //经过前面的消元，除主元外，其他元应该都为 0
if (pcf[b + columns_1] == 0)
{ //常数项为 0，方程有不定解
free(pcf);
return 2;
}
else
{ //常数项为 0，方程有无解
free(pcf); //释放内存
return -1;
}
}
pcf[b + columns_1] /= pcf[b + r];
pcf[b + r] = 1; //肯定为 1，不用计算。
y = -pcf[b + columns_1];
//回代
for (i = r - 1; i = 0; i--)
{
pcf[i * columns + columns_1] += pcf[i * columns + r] * y;
pcf[i * columns + r] = 0; //已经回代，此项已消，置为 0。
}
}
//处理第一行数据
pcf[columns_1] /= pcf[0];
pcf[0] = 1;
//至此，回代过程结束，pcf 矩阵的最后一列就是结果
//返回结果
for (r = 0; r rows; r++)
pOut[r] = pcf[r * columns + columns_1];
free(pcf);
return 0;
}
下面共有两个程序，程序2 加入了图形显示
程序1
这个程序就是你要的。
# include "stdio.h"
# include "math.h"
# include "stdlib.h"
# include "math.h"
# include "dos.h"
# define MAX_N 3000 /*这个值为N可以定义的最大长度*/
# define N 100 /*产生随机序列的点数，注意不要大于MAX_N*/
/*产生均匀分布的随机变量*/
void randa(float *x,int num);
/*产生瑞利分布的随机变量*/
void randr(float *x,int num);
/*产生标准高斯分布的随机变量*/
void randn(float *x,int num);
/*产生莱斯分布的随机变量*/
void randl(float *x, float a, float b, int num);
void fshow(char *name,float *x,int num);
main()
{
float x[N];
int i;
/*
randa(x,N);
randr(x,N);
randl(x,10,10,N);
*/
randn(x,N);
/*此时x[N]就是所需要的高斯分布的序列*/
/*显示该序列*/
fshow("x",x,N);
getch();
}
void randa(float *x,int num)
{
int i;
struct time stime;
unsigned seed;
gettime(stime);
seed=stime.ti_hund*stime.ti_min*stime.ti_hour;
srand(seed);
for(i=0;inum;i++)
{
x[i]=rand();
x[i]=x[i]/32768;
}
}
void randr(float *x,int num)
{
float x1[MAX_N];
int i;
struct time stime;
unsigned seed;
gettime(stime);
seed=stime.ti_hund*stime.ti_min*stime.ti_hour;
srand(seed);
for(i=0;inum;i++)
{
x1[i]=rand();
x[i]=x1[i]/32768;
x[i]=sqrt(-2*log(x[i]));
}
}
void randn(float *x,int num)
{
float x1[MAX_N],x2[MAX_N];
int i;
struct time stime;
unsigned seed;
gettime(stime);
seed=stime.ti_hund*stime.ti_min*stime.ti_hour;
srand(seed);
for(i=0;inum;i++)
{
x1[i]=rand();
x2[i]=rand();
x1[i]=x1[i]/32768;
x2[i]=x2[i]/32768;
x[i]=sqrt(-2*log(x1[i]))*cos(x2[i]*M_PI);
}
}
void randl(float *x, float a, float b, int num)
{
float x1[MAX_N],x2[MAX_N];
float temp[MAX_N];
int i;
struct time stime;
unsigned seed;
gettime(stime);
seed=stime.ti_hund*stime.ti_min*stime.ti_hour;
srand(seed);
for(i=0;inum;i++)
{
x1[i]=rand();
x2[i]=rand();
x1[i]=x1[i]/32768;
x2[i]=x2[i]/32768;
temp[i]=sqrt(-2*log(x1[i]))*cos(x2[i]*M_PI);
x2[i]=sqrt(-2*log(x1[i]))*sin(x2[i]*M_PI);
x1[i]=temp[i];
x[i]=sqrt((a+x1[i])*(a+x1[i])+(b+x2[i])*(b+x2[i]));
}
}
void fshow(char *name,float *x,int num)
{
int i,sign,L;
float temp;
printf("\n");
printf(name);
printf(" = ");
L=6;
/*按照每行6个数据的格式显示*/
for(i=0;inum;i++)
{
temp=i/L;
sign=temp;
if((i-sign*L)==0) printf("\n");
if(x[i]0) printf(" %f ",x[i]);
else printf("%f ",x[i]);
}
}
程序 2
以下程序加入了图形显示的效果，因此更加直观，你可以参考一下。
/* 作者 Leo_nanjing
时间 2008.5.10
功能 生成各种分布的随机变量,并显示
*/
# include "stdio.h"
# include "math.h"
# include "graphics.h"
# include "math.h"
# include "dos.h"
# define MAX_N 3000
# define N 1000
void randa(float *x,int num);
void randr(float *x,int num);
void randn(float *x,int num);
void randl(float *x, float a, float b, int num);
void fshow(char *name,float *x,int num);
/*用于图形显示的部分*/
void init_graphic(unsigned color);
void plotxy(float *x, float *y, int num,int mode);
void plot(float *y,int num, int mode);
float max(float *x, int num);
float min(float *x, int num);
/*画出该随机序列的分布函数曲线*/
void plotpdf(float *x,int num,int part,int mode);
main()
{
float x[N];
int i;
randn(x,N);
fshow("x",x,N);
getch();
/*以下为图形显示部分*/
init_graphic(0);
/*显示随机序列*/
plot(x,N,1);
getch();
/*显示其分布函数*/
plotpdf(x,N,20,0);
getch();
}
void randa(float *x,int num)
{
int i;
struct time stime;
unsigned seed;
gettime(stime);
seed=stime.ti_hund*stime.ti_min*stime.ti_hour;
srand(seed);
for(i=0;inum;i++)
{
x[i]=rand();
x[i]=x[i]/32768;
}
}
void randr(float *x,int num)
{
float x1[MAX_N];
int i;
struct time stime;
unsigned seed;
gettime(stime);
seed=stime.ti_hund*stime.ti_min*stime.ti_hour;
srand(seed);
for(i=0;inum;i++)
{
x1[i]=rand();
x[i]=x1[i]/32768;
x[i]=sqrt(-2*log(x[i]));
}
}
void randn(float *x,int num)
{
float x1[MAX_N],x2[MAX_N];
int i;
struct time stime;
unsigned seed;
gettime(stime);
seed=stime.ti_hund*stime.ti_min*stime.ti_hour;
srand(seed);
for(i=0;inum;i++)
{
x1[i]=rand();
x2[i]=rand();
x1[i]=x1[i]/32768;
x2[i]=x2[i]/32768;
x[i]=sqrt(-2*log(x1[i]))*cos(x2[i]*M_PI);
}
}
void randl(float *x, float a, float b, int num)
{
float x1[MAX_N],x2[MAX_N];
float temp[MAX_N];
int i;
struct time stime;
unsigned seed;
gettime(stime);
seed=stime.ti_hund*stime.ti_min*stime.ti_hour;
srand(seed);
for(i=0;inum;i++)
{
x1[i]=rand();
x2[i]=rand();
x1[i]=x1[i]/32768;
x2[i]=x2[i]/32768;
temp[i]=sqrt(-2*log(x1[i]))*cos(x2[i]*M_PI);
x2[i]=sqrt(-2*log(x1[i]))*sin(x2[i]*M_PI);
x1[i]=temp[i];
x[i]=sqrt((a+x1[i])*(a+x1[i])+(b+x2[i])*(b+x2[i]));
}
}
void fshow(char *name,float *x,int num)
{
int i,sign,L;
float temp;
printf("\n");
printf(name);
printf(" = ");
L=6;
for(i=0;inum;i++)
{
temp=i/L;
sign=temp;
if((i-sign*L)==0) printf("\n");
if(x[i]0) printf(" %f ",x[i]);
else printf("%f ",x[i]);
}
}
/*以下为图形显示的函数*/
void init_graphic(unsigned color)
{
int graphicdriver,graphicmode;
graphicdriver=DETECT;
graphicmode=1;
initgraph(graphicdriver,graphicmode,"E:\\turboc2\\");
setbkcolor(color);
}
void plotxy(float *x, float*y, int num,int mode)
{
int i;
float max_x,max_y,min_x,min_y;
float x0,y0,x1,y1;
clrscr(0);
cleardevice();
setbkcolor(0);
max_x=max(x,num);
max_y=max(y,num);
min_x=min(x,num);
min_y=min(y,num);
setlinestyle(0,2,1);
line(65,35,65,445);
line(65,445,575,445);
setlinestyle(3,0,1);
line(65,35,575,35);
line(575,35,575,445);
setlinestyle(0,2,1);
if(max_x==min_x)
x0=320;
else
x0=(x[0]-min_x)*500/(max_x-min_x)+70;
if(max_y==min_y)
y0=240;
else
y0=480-((y[0]-min_y)*400/(max_y-min_y)+40);
if(mode==0) circle(x0,y0,2);
for(i=1;inum;i++)
{
if(max_x==min_x)
x1=320;
else
x1=(x[i]-min_x)*500/(max_x-min_x)+70;
if(max_y==min_y)
y1=240;
else
y1=480-((y[i]-min_y)*400/(max_y-min_y)+40);
if(mode==0) circle(x1,y1,2);
line(x0,y0,x1,y1);
x0=x1;y0=y1;
}
printf("\n\n");
printf("%f",max_y);
printf("\n\n\n\n\n\n\n\n\n\n");
printf("\n\n\n");
printf("%f",(max_y+min_y)/2);
printf("\n\n\n\n\n\n\n\n\n\n");
printf("\n\n");
printf("%f",min_y);
printf("\n %f",min_x);
printf(" ");
printf("%f",(max_x+min_x)/2);
printf(" ");
printf("%f",max_x);
}
void plot(float*y, int num,int mode)
{
int i;
float max_x,max_y,min_x,min_y;
float x0,y0,x1,y1;
float x[MAX_N];
clrscr(0);
cleardevice();
setbkcolor(0);
for(i=0;inum;i++) x[i]=i+1;
max_x=max(x,num);
max_y=max(y,num);
min_x=min(x,num);
min_y=min(y,num);
setlinestyle(0,2,1);
line(65,35,65,445);
line(65,445,575,445);
setlinestyle(3,0,1);
line(65,35,575,35);
line(575,35,575,445);
setlinestyle(0,2,1);
if(max_x==min_x)
x0=320;
else
x0=(x[0]-min_x)*500/(max_x-min_x)+70;
if(max_y==min_y)
y0=240;
else
y0=480-((y[0]-min_y)*400/(max_y-min_y)+40);
if(mode==0) circle(x0,y0,2);
for(i=1;inum;i++)
{
if(max_x==min_x)
x1=320;
else
x1=(x[i]-min_x)*500/(max_x-min_x)+70;
if(max_y==min_y)
y1=240;
else
y1=480-((y[i]-min_y)*400/(max_y-min_y)+40);
if(mode==0) circle(x1,y1,2);
line(x0,y0,x1,y1);
x0=x1;y0=y1;
}
printf("\n\n");
printf("%f",max_y);
printf("\n\n\n\n\n\n\n\n\n\n");
printf("\n\n\n");
printf("%f",(max_y+min_y)/2);
printf("\n\n\n\n\n\n\n\n\n\n");
printf("\n\n");
printf("%f",min_y);
printf("\n %f",min_x);
printf(" ");
printf("%f",(max_x+min_x)/2);
printf(" ");
printf("%f",max_x);
}
void plotpdf(float *x,int num,int part,int mode)
{
int i,j;
float max_x,min_x,round,deltax,up,down,sum;
float xl[MAX_N],yl[MAX_N];
sum=0;
max_x=max(x,num);
min_x=min(x,num);
round=max_x-min_x;
deltax=round/part;
xl[0]=min_x;
for(i=1;i=part;i++)
{
xl[i]=min_x+deltax*i;
yl[i-1]=0;
up=xl[i];
down=xl[i-1];
for(j=0;jnum;j++)
{
if((x[j]up) (x[j]=down)) yl[i-1]=yl[i-1]+1;
}
yl[i-1]=yl[i-1]/num/deltax;
}
for(i=0;ipart;i++) sum=sum+yl[i];
plotxy(xl,yl,part,mode);
}
float max(float *x, int num)
{
int i;
float max;
max=x[0];
for(i=1;inum;i++)
{
if(x[i]max) max=x[i];
}
return max;
}
float min(float *x, int num)
{
int i;
float min;
min=x[0];
for(i=1;inum;i++)
{
if(x[i]min) min=x[i];
}
return min;
}
数学函数库，一些数学计算的公式的具体实现是放在math.h里，具体有：
1、 三角函数
double sin(double);正弦
double cos(double);余弦
double tan(double);正切
2 、反三角函数
double asin (double); 结果介于[-PI/2,PI/2]
double acos (double); 结果介于[0,PI]
double atan (double); 反正切（主值），结果介于[-PI/2,PI/2]
double atan2 (double,double); 反正切（整圆值），结果介于[-PI,PI]
3 、双曲三角函数
double sinh (double);
double cosh (double);
double tanh (double);
4 、指数与对数
double frexp(double value,int *exp);这是一个将value值拆分成小数部分f和（以2为底的）指数部分exp，并返回小数部分f，即f*2^exp。其中f取值在0.5~1.0范围或者0。
double ldexp(double x,int exp);这个函数刚好跟上面那个frexp函数功能相反，它的返回值是x*2^exp
double modf(double value,double *iptr);拆分value值，返回它的小数部分，iptr指向整数部分。
double log (double); 以e为底的对数
double log10 (double);以10为底的对数
double pow(double x,double y);计算以x为底数的y次幂
float powf(float x,float y); 功能与pow一致，只是输入与输出皆为浮点数
double exp (double);求取自然数e的幂
double sqrt (double);开平方
5 、取整
double ceil (double); 取上整，返回不比x小的最小整数
double floor (double); 取下整，返回不比x大的最大整数，即高斯函数[x]
6 、绝对值
int abs(int i); 求整型的绝对值
double fabs (double);求实型的绝对值
double cabs(struct complex znum);求复数的绝对值
7 、标准化浮点数
double frexp (double f,int *p); 标准化浮点数，f = x * 2^p，已知f求x,p (x介于[0.5,1])
double ldexp (double x,int p); 与frexp相反，已知x,p求f
8 、取整与取余
double modf (double,double*); 将参数的整数部分通过指针回传，返回小数部分
double fmod (double,double); 返回两参数相除的余数
