【算法】算法中的趣味数学（三）
求解线性方程
用高斯(Guass)消去法求解N阶线性方程组Ax=B。
实例解析：
高斯消去法解线性代数方程的基本原理如下。
对于线性方程组：
其中系数矩阵为A，未知量为X，值向量为B。计算的方法分为两步进行。
第１步，消去过程，对于k从0到n -2做以下3步。
从系数矩阵A的第k行、第k列开始的右下角子阵中选取绝对值最大的元素，并通过行交换与列交换把它交换到主元素(即对角线元素)的位置上。
归一法： j=k+1,…..,n-1
消去： j=k+1,…..,n-1
第２步，回代过程。
i=n-2,…..,1,0
最后对解向量中的元素顺序进行调整。
程序如下：
#include "stdlib.h"
#include "math.h"
#include "stdio.h"
#define MAX 255
int Guass(double a[], double b[], int n)
{
int *js, m, k, i, j, is, p, q;
double d, t;
js = (int *)malloc(n * sizeof(int));
m = 1;
for(k = 0; k <= n-2; k++)
{
d = 0.0;
/*下面是换主元部分，即从系数矩阵A的第K行，第K列之下的部分选出
绝对值最大的元，交换到对角线上。*/
for(i = k; i <= n-1; i++)
{
for(j = k; j <= n-1; i++)
{
t = fabs(a[i*n+j]);
if(t > d)
{
d = t;
js[k] = j;
is = i;
}
}
if(d +1.0 == 1.0)
m = 0; //主元为0
else
{
if(js[k] != k)
for(i = 0; i <= n-1; i++)
{
p = i*n + k;
q = i*n + js[k];
t = a[p]; a[p] = a[q];
a[q] = t;
}
if(is != k)
{
for(j = k; j <= n-1; j++)
{
p = k*n + j;
q = is*n + j;
t = a[p]; a[p] = a[q];
a[q] = t;
}
t = b[k]; b[k] = b[is]; b[is] = t;
}
}
}
if(m == 0)
{
free(js);
printf("fail\n");
return(0);
}
d = a[k*n + k];
//下面为归一化部分
for(j = k+1; j <= n-1; j++)
{
p = k*n +j;
a[p] = a[p]/d;
}
b[k] = b[k]/d;
//下面为矩阵A,B消元部分
for(i = k+1; i <= n-1; i++)
{
for(j = k+1; j <= n-1; j++)
{
p = i*n + j;
a[p] = a[p] - a[i*n+k]*a[k*n+j];
}
b[i] = b[i] - a[i*n+k]*b[k];
}
}
d = a[(n-1)*n + n-1];
//矩阵无解或有无限多解
if(fabs(d) + 1.0 == 1.0)
{
free(js);
printf("该矩阵为奇异矩阵\n");
return(0);
}
b[n-1] = b[n-1]/d;
//下面为迭代消元
for(i = n-2; i >= 0; i--)
{
t = 0.0;
for(j = i+1; j <= n-1; j++)
t = t + a[i*n+j]*b[j];
b[i] = b[i] - t;
}
js[n-1] = n-1;
for(k = n-1; k >= 0; k--)
if(js[k] != k)
{
t = b[k]; b[k] = b[js[k]]; b[js[k]] = t;
}
free(js);
return 1;
}
int main()
{
int i,n;
double A[MAX];
double B[MAX];
printf(">>Please input the order n (>1): ");
scanf("%d",&n);
printf(">>Please input the %d elements of atrix A(%d*%d) \n", n*n,n,n);
for(i = 0; i < n*n; i++)
scanf("%lf",&A[i]);
printf(">>Please input the %d elements of matrix B(%d*1) one by one:\n",n,n);
for(i = 0; i < n; i++)
scanf("%lf", &B[i]);
if (Guass(A,B,n) != 0) //调用Guass()，1为计算成功
printf(" >> The solution of Ax=B is x(%d*1):\n",n);
for(i = 0; i < n; i++) //打印结果
printf("x(%d)=%f ", i, B[i]);
puts("\n Press any key to quit...");
getch();
return 0;
}
求积分
用变长辛普森法求定积分。
实例解析：：
用变长辛普森(Simpson)法则求定积分值的实现方法如下。
(1) 用梯形公式计算，其中n = 1，h = b - a，且令Sn = Tn
(2) 用变步长梯形法则计算
(3) 用辛普森求积公式计算S2n = (4T2n - Tn) /3
若｜S2n-Sn|≥ε，则令2n→n, h/2→h，转到步骤(2)继续进行计算；否则结束，S2n即为所求积分的近似值。其中ε为事先给定的求积分精度。
在本例中，使用辛普森法则计算定积分：。精度ε＝0.000001.
程序如下：
#include "stdio.h"
#include "math.h"
double fsimpf(double x) //要进行计算的被积函数
{
double y;
y = log(1.0+x)/(1.0+x*x);
return(y);
}
double fsimp(double a,double b,double eps) //辛普森算法
//a为积分下限，b为积分上限，eps是希望达到的精度
{
int n, k;
double h, t1, t2, s1, s2, ep, p, x;
n = 1;
h = b - a;
//用梯形公式求出一个大概的估值
t1 = h*(fsimpf(a) + fsimpf(b))/2.0;
s1 = t1;
ep = eps + 1.0;
while(ep >= eps)
{ //用梯形法则计算
p = 0.0;
for(k = 0; k <= n-1; k++)
{
x = a + (k + 0.5)*h;
p = p + fsimpf(x);
}
t2 = (t1 + h*p)/2.0;
//用辛普森公式求精
s2 = (4.0*t2-t1)/3.0;
ep = fabs(s2-s1);
t1 = t2;
s1 = s2;
n = n + n;
h = h/2.0;
}
return s2;
}
int main()
{
double a,b,eps,t;
a = 0.0;
b = 1.0;
eps = 0.0000001;
t = fsimp(a, b, eps);
puts("\n----------------------------------------");
printf(">>The result of definite integral is : \n");
printf(">>SIGMA(0,1)ln(1+x)/(1+x^2)dx = ");
printf("%e\n",t);
puts("-------------------------------------------");
printf("\n Press any key to quit...");
getch();
return 0;
}
超长整数的加法
实现超长正整数的加法运算。
实例解析：
首先设计一种数据结构来表示一个超长的正整数，然后才能够设计算法。
首先采用一个带头结点的环形链来表示一个非负的超大整数，如果从低位开始为每个数字编号，则第1～4位、第5～8位……的每4位组成的数字，依次放到链表的第1个、第2个……结点中，不足4位的最高位存在链表的最后一个结点中，头结点的值规定为-1。例如：大整数“587890987654321”
按照此数据结构，可以从两个头结点开始，顺序依次对应相加，求出所需要的进位后，代入下一个结点运算。
程序如下：
#include<stdio.h>
#include<stdlib.h>
#define HUNTHOU 10000
typedef struct node
{
int data;
struct node *next;
}NODE; //定义链表结构
struct number
{
int num;
struct number *np;
};
void freelist(NODE* llist);
NODE *insert_after(NODE *u,int num); //在u结点后插入
NODE *AddInt(NODE *p,NODE *q);//完成加法操作返回指向结果的指针
void printint(NODE *s);
NODE *inputint(void);
int main()
{
NODE *s1,*s2,*s;
printf(">> Input S1= ");
s1 = inputint(); //输入被加数
if(s1 == NULL)
{
return 0;
}
printf(">> Input S2= ");
s2 = inputint(); //输入加数
if(s2 == NULL)
{
freelist(s1);
return 0;
}
printf(">> The addition result is as follows.\n\n");
printf(" S1=");
printint(s1); //显示被加数
putchar('\n');
printf(" S2=");
printint(s2); //显示加数
putchar('\n');
s = AddInt(s1,s2); //求和
if(s == NULL)
{
freelist(s1);
freelist(s2);
return 0;
}
printf(" S1+S2=");
printint(s); //输出结果
putchar('\n');
freelist(s1);
freelist(s2);
freelist(s);
printf("\n\n Press any key to quit...");
return 0;
}
NODE *insert_after(NODE *u,int num)
{
NODE *v;
v = (NODE *)malloc(sizeof(NODE)); //申请一个NODE
if( v == NULL)
{
return NULL;
}
v->data = num; //赋值
u->next = v; //在u结点后插入一个NODE
return v;
}
NODE *AddInt(NODE *p,NODE *q) //返回指向*p+*q结果的指针
{
NODE *pp,*qq,*r,*s,*t,*tmp;
int total,number,carry;
pp = p->next;
qq = q->next;
s = (NODE *)malloc(sizeof(NODE)); //建立存放和的链表头
if( s == NULL)
{
return NULL;
}
s->data = -1;
t= tmp = s;
carry = 0; //carry: 进位
while(pp->data != -1 && qq->data != -1)
{ //均不是头结点
total = pp->data + qq->data + carry; //对应位求和
number = total%HUNTHOU; //求出存入链中部分的数值
carry = total/HUNTHOU; //算出进位
tmp = insert_after(t, number); //将部分和存入s指向的链中
if(tmp == NULL)
{
t->next = s;
freelist(s);
return NULL;
}
t = tmp;
pp = pp->next; //分别取后面的加数
qq = qq->next;
}
r = (pp->data != -1) ? pp:qq; //取尚未自理完毕的链指针
while(r->data != -1)
{ //处理加数中较大的数
total = r->data + carry; //与进位相加
number = total%HUNTHOU; //求出存入链中部分的数值
carry = total/HUNTHOU; //算出进位
tmp = insert_after(t,number); //将部分和存入s指向的链中
if(tmp == NULL)
{
t->next = s;
freelist(s);
return NULL;
}
t = tmp;
r = r->next; //取后面的值
}
if(carry)
tmp = insert_after(t, 1); //处理最后一次进位
if(tmp == NULL)
{
t->next = s;
freelist(s);
return NULL;
}
t = tmp;
t->next = s; //完成存和的链表
return s; //返回指向和的结构指针
}
NODE *inputint(void) //输入超长正整数
{
NODE *s,*ps,*qs;
struct number *p,*q,*pre;
int i ,k;
long sum;
char c;
p = NULL;//用来指向输入的整数，链首为整数的最低位，链尾为最高位
while((c=getchar()) != '\n') //输入整数，按字符接收数字
{
if(c >= '0' && c <= '9')
{ //若为数字则存入
q = (struct number*)malloc(sizeof(struct number));
if(q == NULL)
{
freelist2(p);
return NULL;
}
q->num = c-'0'; //存入一位整数
q->np = p; //建立指针
p = q;
}
s = (NODE *)malloc(sizeof(NODE));
if(s == NULL)
{
freelist2(p);
return NULL;
}
s->data = -1; //建立表求超长正整数的链头
ps = s;
while(p != NULL)
{ //转换
sum = 0;
i = 0;
k = 1;
while(i < 4 && p != NULL)
{ //取出低四位
sum = sum + k*(p->num);
i++;
pre = p;
p = p->np;
k = k*10;
free(pre); //释放前一个已经用完的结点
}
qs = (NODE *)malloc(sizeof(NODE));
if(qs == NULL)
{
ps->next = s;
freelist(s);
return NULL;
}
qs->data = sum; //赋值，建立链表
ps->next = qs;
ps = qs;
}
ps->next = s;
return s;
}
}
void printint(NODE *s)
{
int i, k;
if(s->next->data != -1)
{ //若不是头结点，则输出
printint(s->next); //递归输出
if(s->next->next->data == -1)
printf("%d", s->next->data);
else
{
k = HUNTHOU;
for(i = 1; i <= 4; i++,k/=10)
putchar('0' + s->next->data % k/(k/10));
}
}
}
void freelist(NODE* llist)
{
NODE* pnode = llist->next;
NODE* p;
while(pnode != llist)
{
p = pnode;
pnode = pnode->next;
free(p);
}
}
void freelist2(struct number* llist)
{
struct number *p;
while(llist)
{
p = llist;
llist = llist->np;
free(p);
}
}
本文出自 “成鹏致远” 博客，请务必保留此出处http://infohacker.blog.51cto.com/6751239/1171363
