求解线性方程

　　用高斯(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
　　　最后对解向量中的元素顺序进行调整。
　　　程序如下：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
#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.
　　程序如下：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
#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”
　　按照此数据结构，可以从两个头结点开始，顺序依次对应相加，求出所需要的进位后，代入下一个结点运算。
　　程序如下：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
#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);
     }
}