【算法】算法中的趣味数学（三）

求解线性方程

　　用高斯(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 255int 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 10000typedef 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