1、數值分析程序設計Part IIFortran程序設計1 線性代數1.1 矩陣的加法和乘法矩陣作為數值分析的重要工具，在數值計算中具有非常重要的作用，在Fortran語言中利用數組來表示矩陣。在Fortran 90中可以直接對數組進行計算，一個命令就可以完成矩陣的加法、減法。Real a(m,n), b(m,n),c(m,n)C=a+b ！矩陣加法C=a-b ！矩陣減法在Fortran 77中必須利用循環完成矩陣的加減法：do r=1,m do c=1,n c(r,c)=a(r,c)+b(r,c) end doend do矩陣的乘法，在Fortran 90中不能直接用乘號來做，必須調用庫函數。C

2、=a*b ！這個命令是做c(I,j)=a(I,j)*b(I,j)C=matmul(a,b) ！庫函數matmul可以做矩陣相乘在Fortran 77中則要自己寫矩陣乘法程序代碼進行計算，下面是Fortran 77固定格式編寫的矩陣相乘程序代碼：CC 矩陣乘法范例C By Perng 1997/9/17 PROGRAM MATMUL_DEMO IMPLICIT NONE INTEGER N PARAMETER(N=3) INTEGER A(N,N) ! MATRIX A INTEGER B(N,N) ! MATRIX B INTEGER C(N,N) ! MATRIX C DATA B /1,2

3、,3,4,5,6,7,8,9/ DATA C /9,8,7,6,5,4,3,2,1/ CALL MATMUL(A,B,N,N,C,N,N) WRITE(*,*) 'Matrix A:' CALL OUTPUT(A,N) STOP ENDCC 輸出矩陣的子程序C SUBROUTINE OUTPUT(A,N) IMPLICIT NONE INTEGER N,A(N,N) INTEGER I,J CHARACTER FOR*20 DATA FOR /'(?(1X,I3)'/C 用字符串來設定輸出格式 WRITE( FOR(2:3), '(I2)' )

4、N DO I=1,N WRITE( *, FMT=FOR ) (A(I,J),J=1,N) END DO RETURN ENDCC 矩陣乘法的子程序C SUBROUTINE MATMUL(A,B,BR,BC,C,CR,CC) IMPLICIT NONE INTEGER BR ! Row of Matrix B INTEGER BC ! Column of Matrix B INTEGER B(BR,BC) ! Matrix B INTEGER CR ! Row of Matrix C INTEGER CC ! Column of Matrix C INTEGER C(CR,CC) ! Matr

5、ix C INTEGER A(BR,CC) ! Matrix A INTEGER I,J,K ! 循環的計數器 ! BC若不等于CR, 這兩個矩陣無法相乘 IF ( BC .NE. CR ) THEN WRITE(*,*) 'Matrix size error!' STOP END IF DO I=1,BR DO J=1,CC A(I,J)=0 DO K=1,BC A(I,J)=A(I,J)+B(I,K)*C(K,J) END DO END DO END DO RETURN END需要注意的是：假如矩陣一開始聲明了10*10的大小，而所需要使用的矩陣大小只有2*2，那么一定要避

6、免下面的情況。使用下面的方法，在子程序中會讀到錯誤的矩陣內容。Program mianImplicit noneInteger A(10,10)A(1,1)=1.0A(2,1)=2.0A(1,2)=1.0A(2,2)=2.0Call sub(A,2,2).end program mainsubroutine sub(matrix, row, col)implicit noneinteger row, colinteger matrix(row, col)end subroutine sub實際聲明數組大小為10*10，傳遞到子程序中卻把它聲明為2*2，會有下面的對應情況發生：matrix(1,

7、1)=a(1,1)=1.0matrix(2,1)=a(2,1)=2.0matrix(1,2)=a(3,1)=0.0/=a(1,2)matrix(2,2)=a(4,1)=0.0/=a(2,2)這是由于多維數組的存儲方式決定的，在數據讀取時應當注意。1.2 三角矩陣通過矩陣的初等變換，可以將矩陣化成上三角矩陣、下三角矩陣或者對角矩陣。下面程序將矩陣分別化成上三角矩陣和下三角矩陣：=>=>UPPERLOWER.F90module LinearAlgebra implicit nonecontains! 輸出矩陣的子程序subroutine output(matrix) implicit

8、none integer : m,n real : matrix(:,:) integer : i character(len=20) : for='(?(1x,f6.3)' m = size(matrix,1) n = size(matrix,2) ! 用字符串來設定輸出格式 write( FOR(2:3), '(I2)' ) N do i=1,Nwrite( *, FMT=FOR ) matrix(i,:) end do returnend subroutine output! 求上三角矩陣的子程序subroutine Upper(matrix) impli

9、cit none real : matrix(:,:) integer : M,N integer : I,J real : E M=size(matrix,1) N=size(matrix,2) do I=1,N-1do J=I+1,M E=matrix(J,I)/matrix(I,I) ! 用90的功能可以少一層循環 matrix(J,I:M)=matrix(J,I:M)-matrix(I,I:M)*Eend do end do returnend subroutine Upper! 求下三角矩陣的子程序subroutine Lower(matrix) implicit none real

10、 : matrix(:,:) integer : M,N real : I,J,E M = size(matrix,1) N = size(matrix,2) do I=N,2,-1 do J=I-1,1,-1 E=matrix(J,I)/matrix(I,I) ! 用90的功能可以少一層循環 matrix(J,1:I)=matrix(J,1:I)-matrix(I,1:I)*E end do end do returnend subroutine Lowerend moduleprogram main use LinearAlgebra implicit none integer, para

11、meter : N = 3! Size of Matrix real : A(N,N) = reshape( (/1,2,1,3,2,3,2,3,4/),(/N,N/) ) real : B(N,N) write(*,*) "Matrix A:" call output(A) B=A write(*,*) "Upper:" call Upper(B) call output(B) B=A write(*,*) "Lower:" call Lower(B) call output(B) stopend program程序執行結果：矩陣的

12、三角化可以應用于行列式計算、求解逆矩陣、線性方程組求解等方面。1.3 矩陣的行列式計算方陣的行列式計算，采用矩陣三角化方法，矩陣的行列式等于其三角化矩陣對角線元素的乘積。下面程序先把矩陣化成上三角矩陣，再求行列式的值。DETERMINANT.F90module LinearAlgebra implicit nonecontains! 求矩陣的Determinant值real function Determinant(matrix) real : matrix(:,:) real, allocatable : ma(:,:) integer : i,N N = size(matrix,1) al

13、locate(ma(N,N) ma = matrix call Upper(ma) Determinant = 1.0 do i=1,N Determinant = Determinant*ma(i,i) end doend function! 求上三角矩陣的子程序subroutine Upper(matrix) real : matrix(:,:) integer : M,N integer : I,J real : E M=size(matrix,1) N=size(matrix,2) do I=1,N-1do J=I+1,M E=matrix(J,I)/matrix(I,I) ! 用90

14、的功能可以少一層循環 matrix(J,I:M)=matrix(J,I:M)-matrix(I,I:M)*Eend do end do returnend subroutine Upperend moduleprogram main use LinearAlgebra implicit none integer, parameter : N = 3! Size of Matrix real : A(N,N) = reshape( (/1,2,1,3,2,3,2,3,4/),(/N,N/) ) write(*,"('det(A)=',F6.2)") Deter

15、minant(A) stopend program1.4 Gauss-Jordan消去法求解線性方程組線性方程組的求解可以分為直接法和迭代法。直接法有Gauss消去法和直接三角分解法。迭代法有Jacobi迭代法和Gauss-Seidel迭代法等。Gauss-Jordan消去法是直接求解線性方程組的數值方法。通過矩陣的初等變換，將矩陣化成對角矩陣求解。GAUSS-JORDAN.F90module LinearAlgebra implicit nonecontains! Gauss_Jordan法subroutine Gauss_Jordan(A,S,ANS) implicit none real

16、 : A(:,:) real : S(:) real : ANS(:) real, allocatable : B(:,:) integer : i, N N = size(A,1) allocate(B(N,N) ! 保存原先的矩陣A,及數組S B=A ANS=S ! 把B化成對角線矩陣(除了對角線外,都為0) call Upper(B,ANS,N) ! 先把B化成上三角矩陣 call Lower(B,ANS,N) ! 再把B化成下三角矩陣 ! 求解 forall(i=1:N) ANS(i)=ANS(i)/B(i,i) end forall returnend subroutine Gaus

17、s_Jordan! 輸出等式subroutine output(M,S) implicit none real : M(:,:), S(:) integer : N,i,j N = size(M,1) ! write中加上advance="no",可以中止斷行發生,使下一次的 ! write接續在同一行當中. do i=1,N write(*,"(1x,f5.2,a1)", advance="NO") M(i,1),'A' do j=2,N if ( M(i,j) < 0 ) then write(*,"

18、;('-',f5.2,a1)",advance="NO") -M(i,j),char(64+j) else write(*,"('+',f5.2,a1)",advance="NO") M(i,j),char(64+j) end if end do write(*,"('=',f8.4)") S(i) end do returnend subroutine output! 求上三角矩陣的子程序subroutine Upper(M,S,N) implicit n

19、one integer : N real : M(N,N) real : S(N) integer : I,J real : E do I=1,N-1 do J=I+1,N E=M(J,I)/M(I,I) M(J,I:N)=M(J,I:N)-M(I,I:N)*E S(J)=S(J)-S(I)*E end do end do returnend subroutine Upper! 求下三角矩陣的子程序subroutine Lower(M,S,N) implicit none integer : N real : M(N,N) real : S(N) integer : I,J real : E

20、do I=N,2,-1 do J=I-1,1,-1 E=M(J,I)/M(I,I) M(J,1:N)=M(J,1:N)-M(I,1:N)*E S(J)=S(J)-S(I)*E end do end do returnend subroutine Lowerend module! 求解聯立式program main use LinearAlgebra implicit none integer, parameter : N=3 ! Size of Matrix real : A(N,N)=reshape( (/1,2,3,4,5,6,7,8,8/),(/N,N/) ) real : S(N)=(

21、/12,15,17/) real : ans(N) integer : i write(*,*) 'Equation:' call output(A,S) call Gauss_Jordan(A,S,ANS) write(*,*) 'Ans:' do i=1,N write(*,"(1x,a1,'=',F8.4)") char(64+i),ANS(i) end do stopend program程序執行結果：1.5 逆矩陣求解逆矩陣的算法與Gauss-Jordan消去法求解方程類似，在矩陣的右邊添加一個同階的單位矩陣，將左邊

22、矩陣變換為單位矩陣，右邊的單位矩陣變換為矩陣的逆矩陣。程序如下：INVERSE.F90module LinearAlgebra implicit nonecontains! Gauss_Jordan法subroutine Gauss_Jordan(A,S,ANS) implicit none real : A(:,:) real : S(:) real : ANS(:) real, allocatable : B(:,:) integer : i, N N = size(A,1) allocate(B(N,N) ! 保存原先的矩陣A,及數組S B=A ANS=S ! 把B化成對角線矩陣(除了對

23、角線外,都為0) call Upper(B,ANS,N) ! 先把B化成上三角矩陣 call Lower(B,ANS,N) ! 再把B化成下三角矩陣 ! 求解 forall(i=1:N) ANS(i)=ANS(i)/B(i,i) end forall returnend subroutine Gauss_Jordan! 輸出等式subroutine output(M,S) implicit none real : M(:,:), S(:) integer : N,i,j N = size(M,1) ! write中加上advance="no",可以中止斷行發生,使下一次的

24、! write接續在同一行當中. do i=1,N write(*,"(1x,f5.2,a1)", advance="NO") M(i,1),'A' do j=2,N if ( M(i,j) < 0 ) then write(*,"('-',f5.2,a1)",advance="NO") -M(i,j),char(64+j) else write(*,"('+',f5.2,a1)",advance="NO") M(i,j)

25、,char(64+j) end if end do write(*,"('=',f8.4)") S(i) end do returnend subroutine output! 求上三角矩陣的子程序subroutine Upper(M,S,N) implicit none integer : N real : M(N,N) real : S(N) integer : I,J real : E do I=1,N-1 do J=I+1,N E=M(J,I)/M(I,I) M(J,I:N)=M(J,I:N)-M(I,I:N)*E S(J)=S(J)-S(I)*E

26、end do end do returnend subroutine Upper! 求下三角矩陣的子程序subroutine Lower(M,S,N) implicit none integer : N real : M(N,N) real : S(N) integer : I,J real : E do I=N,2,-1 do J=I-1,1,-1 E=M(J,I)/M(I,I) M(J,1:N)=M(J,1:N)-M(I,1:N)*E S(J)=S(J)-S(I)*E end do end do returnend subroutine Lowerend module! 求解聯立式prog

27、ram main use LinearAlgebra implicit none integer, parameter : N=3 ! Size of Matrix real : A(N,N)=reshape( (/1,2,3,4,5,6,7,8,8/),(/N,N/) ) real : S(N)=(/12,15,17/) real : ans(N) integer : i write(*,*) 'Equation:' call output(A,S) call Gauss_Jordan(A,S,ANS) write(*,*) 'Ans:' do i=1,N w

28、rite(*,"(1x,a1,'=',F8.4)") char(64+i),ANS(i) end do stopend program2 數值積分對于積分進行數值積分，就是利用被積函數在區間上的一些節點處的函數值的線性組合去近似。其中稱為積分節點，稱為積分系數，為近似求積公式。積分節點和積分系數的不同決定了不同的積分公式。常用的有梯形積分公式、Simpson積分公式和Gauss積分公式等。2.1 梯形公式和Simpson公式梯形公式和Simpson公式是插值型積分公式。梯形公式是用一次函數逼近被積函數的積分公式，而Simpson公式是利用二次函數逼近被積函數

29、的積分公式。梯形積分公式：。計算TRAPE.F90module INTEGRAL implicit none real, parameter : PI=3.14159contains! 產生數列 subroutine GenerateData(datas, width, func) real datas(:), widthreal, external : funcreal rinteger i,nn = size(datas,1)width = PI/(n-1)r = 0do i=1,n datas(i) = func(r) r = r+widthend do end subroutine!

30、梯形法積分 real function Trape_Integral(datas, width) implicit none real datas(:) real width ! 每條數據的間隔 real SUM ! 計算所有上底加下底除以二的和 integer i,n n = size(datas,1) SUM = (datas(1)+datas(n)/2.0 do i=2,n-1 SUM=SUM+datas(i) ! 累加邊長 end do Trape_Integral=SUM*width ! 計算面積和 return end function end module! 梯形法積分范例pro

31、gram main use INTEGRAL implicit none integer, parameter : N = 10 real DATAS(N), width real ANS ! 答案 real, intrinsic : sin ! 模擬用來產生數據的函數 call GenerateData(DATAS, width, sin) ANS = Trape_Integral(DATAS, width) ! 計算積分 write(*,"('ans=',F5.2)") ANS ! 寫出答案 stopend programSimpson 積分公式：Sim

32、pson.F90module INTEGRAL implicit none real, parameter : PI=3.14159contains! 產生數列 subroutine GenerateData(datas, width, func) real datas(:), widthreal, external : funcreal rinteger i,nn = size(datas,1)width = PI/(n-1)r = 0do i=1,n datas(i) = func(r) r = r+widthend do end subroutine! 梯形法積分real functio

33、n Simpson_Integral(datas, width) IMPLICIT NONE real datas(:), width real sum integer i,n n = size(datas,1) if ( mod(n,2)=0 ) then write(*,*) "要有奇數條數據"stop end if sum = datas(1) + datas(n) ! 先算出頭尾的和 do i=2,n-1 if ( mod(i,2)=0 ) then sum = sum + 4*datas(i) ! 把4*f(x)的部分累加起來else sum = sum + 2*

34、datas(i) ! 把2*f(x)的部分累加起來end if end do Simpson_Integral = sum * width/3.0 ! SUM再乘上H/3 就好了 returnend functionend module! SIMPSON法積分范例program main use INTEGRAL implicit none integer, parameter : N = 9 real, intrinsic : sin real datas(N), width call GenerateData(datas, width, sin) write(*,"('a

35、ns=',F6.2)") simpson_integral(datas, width) stopend program3 插值與擬合3.1 Lagrange Interpolation所謂插值，就是給定一組點和函數值，求一個函數，使得。插值一般表示成基函數的形式：，其中成作插值基函數，其滿足，。Lagrange插值的基函數為：給定正弦函數的10個值，采用拉格朗日插值法求100個插值節點的值。LAGRANGE.F90module INTERPOLATE_UTILITY implicit none type point real x,y end type real, parame

36、ter : PI=3.14159 real, parameter : xmin = 0.0, xmax = PI*3.0 integer, parameter : N = 10, NP = 100 type(point) : datas(N) type(point) : interpolate(NP)contains! 產生數列 subroutine GenerateData(func)real, external : funcreal r, widthinteger iwidth = (xmax-xmin)/(N-1)r = 0do i=1,N datas(i)%x = r datas(i)

37、%y = func(r) r = r+widthend do end subroutine real function lagrange(x)real x real coeffinteger i,jlagrange = 0 do i=1,n coeff = 1 do j=1,n if ( i/=j ) coeff = coeff * (x-datas(j)%x)/(datas(i)%x-datas(j)%x) end do lagrange = lagrange + coeff*datas(i)%yend do end functionend moduleprogram main use IN

38、TERPOLATE_UTILITY implicit none real, intrinsic : sin real xinc,x integer i call GenerateData(sin) ! 產生數據點 x=0 xinc = (xmax-xmin)/(NP-1) do i=1,NP interpolate(i)%x = x interpolate(i)%y = lagrange(x) ! 插值出f(x)的值x = x+xinc end dowrite(*,"(4X'X',5X)") do I=1,Np write(*,"(F5.3,F10

39、.5)") X(I),Y(I) end do stopend program3.2 最小二乘法（Least-Square Method）最小二乘法用于數據處理，利用給定的數據擬合出函數曲線。與插值法不同的是，擬合曲線一般不通過給定的數據點，數據點分布在曲線的附近。曲線擬合根據事先假定的函數關系不同可以分為：線性和非線性擬合。按照變量的多少可以分為一元和多元擬合。假設一個銅條在0度時長1米，由于溫度的不同熱脹冷縮的長度變化如下：溫度T（C）長度Y（米）51.047101.112151.152201.191251.252假設銅的長度關于溫度是線性膨脹的，求膨脹系數。根據已知條件，可以得到

40、長度和溫度的關系式為。利用最小二乘法的思想是：使得數據點到直線的距離最小。為了便于數學處理，采用距離平方和最小。N個點到直線的距離平方和為：要使距離平方和最小，將上式分別對求偏導數，并令其等于零，求出。上式等價于寫成矩陣的形式為求解方程組，得到系數。程序實現如下：LEAST_SQUARE.F90module datas implicit none integer, parameter : N=5 real : temperature(N) = (/5.0,10.0,15.0,20.0,25.0/) ! 記錄溫度 real : length(N) = (/1.047,1.112,1.1152,1

41、.191,1.252/)! 記錄不同溫度下的長度 real, save : A,B ! 函數L=At+B的系數end module! 最小方差法范例program main use datas implicit none ! 調用最小方差法的子程序來計算系數A,B call least_square(temperature, length, N, A, B) write(*,"('A=',F6.4,' B=',F6.4)") A,B stopENDsubroutine least_square(x, y, n, s, t) implicit

42、none integer n real x(n) ! x上的數據點 real y(n) ! y上的數據點 real s,t ! 所要計算的系數 real A,B,C,D,E,F ! 聯立方程式中的系數 integer I ! 循環計數器 ! 解 As+Bt=E 中的未知數s,t ! Cs+Dt=F ! 先設定好A,B,C,D,E,F的系數值 A = 0; B = 0; E = 0; F = 0 do I=1,N A=A+X(I)*X(I) B=B+X(I) E=E+X(I)*Y(I) F=F+Y(I) END DO C=B D=N ! 二元一次方程式有公式可解 S=(B*F-E*D)/(B*C

43、-A*D) T=(E*C-A*F)/(B*C-A*D) returnend subroutine程序執行結果：4 數據結構與算法4.1 排序4.1.1 冒泡排序法冒泡排序法是最簡單的排序方法，他的步驟如下：（1）從第一個數字開始，依次把兩個相鄰的數據互相比較大小。如果前一個數字比后一個數字大，就交換他們的位置；（2）一直做到每一對相鄰數字都比較過后，結束這一輪工作；（3）回到第一步，再作下一循環的比較。如果有n個數字要排序，就需要重復n-1次的掃描工作。冒泡排序法可以想象成讓重的東西向下沉，輕的東西向上浮。等到狀態穩定，就得到排序的結果。下面是冒泡排序法的程序實現：! 冒泡排序法范例! By

44、Perng 1997/8/29program BUBBLE_SORT_DEMO implicit none integer, parameter : N=10 integer : A(N)=(/6,2,8,4,0,9,3,5,1,7/) ! 待排序的數據 write(*,"('Source=>',10I3)") A call BUBBLE_SORT(A,N) ! 調用排序的子程序 write(*,"('Sort=>',10I3)") A stopend programsubroutine BUBBLE_SORT

45、(A,N) implicit none integer : N,A(N) integer I,J, TEMP do I=N-1,1,-1 ! 開始做N-1次的掃瞄 do J=1,I ! 一對一對的來比較，I之后的數字不用比較 ! 如果A(J) > A(J+1) 就把這兩個數值交換 if ( A(J) > A(J+1) ) then TEMP=A(J) A(J)=A(J+1) A(J+1)=TEMP end if end do end do returnend subroutine4.1.2 選擇排序法選擇排序法的基本原理為：（1）找出全部n個數據中最小的一個，把它和數列的第一個數字

46、交換位置；（2）找出剩余n-1個數據中最小的一個，把它和數列的第二個數字交換位置；（3）找出剩余n-2個數據中最小的一個，把它和數列的第三個數字交換位置；（4）（5）一直做到只剩一個數據為止。計算程序如下：! 選擇排序法范例! By Perng 1997/8/29program SELECTION_SORT_DEMO implicit none integer, parameter : N=10 integer : A(N)=(/6,2,8,4,0,9,3,5,1,7/) ! 排序的數據 write(*,"('Source=>',10I3)") A c

47、all SELECTION_SORT(A,N) ! 調用排序的子程序 write(*,"('Sort=>',10I3)") A stopend program! 選擇排序法的子程序subroutine SELECTION_SORT(A,N) implicit none integer : N,A(N) integer I,J ! 循環計數器 integer MIN ! 找出每一輪中的最小值 integer TEMP ! 交換數據時使用 do I=1,N MIN=A(I) ! 暫時令A(I)是最小值 do J=I+1,N if ( MIN > A(

48、J) ) then ! 發現A(I)不是最小 TEMP=A(J) ! 把A(I)、A(J)交換 A(J)=A(I) A(I)=TEMP MIN=A(I) end ifend do end do returnend subroutine 4.2 搜索4.2.1 順序搜索法順序搜索是一種最簡單的搜索方法，其原理是將數據一個一個拿出來，看他是不是我們所需要的數據。順序搜索法的程序如下：! 順序查找法范例! By Perng 1997/8/31program SEQUENTIAL_SEARCH_DEMO implicit none integer, parameter : N=10 integer :

49、 A(N)=(/6,2,8,4,0,9,3,5,1,7/) ! 存放數據組的類型 integer KEY ! 記錄所要找的值 integer LOC integer, external : SEQUENTIAL_SEARCH write(*,"('Source=>',10I3)") A write(*,*) 'Input KEY:' read (*,*) KEY ! 鍵入待尋數據 ! 調用順序查找的函數 LOC = SEQUENTIAL_SEARCH(A,N,KEY) if ( LOC/=0 ) then write(*,"(

50、'A(',I2,' )='I3)") LOC,KEY else write(*,*) "Not found" end if stopend program! 順序查找法的子程序integer function SEQUENTIAL_SEARCH(A,N,KEY) implicit none integer N, A(N) integer KEY ! 所要尋找的值 integer I ! 循環的計數器 do I=1,N ! 開始做掃瞄, 最多做N次. if ( KEY=A(I) ) then ! 找到了, 返回數字在類型中的位置 SE

51、QUENTIAL_SEARCH=I returnend if end do ! 沒找到時返回-1 SEQUENTIAL_SEARCH=0 returnend function4.2.2 二元搜索二元搜索法必須配合排序好的數據才能使用，假設所要尋找的數據值為K，數據存放在數組中，搜索的步驟為：（1）取出數組的中間值M和K來比較，如果M=K，結束。如果K>M，因為數組早已作好排序，所以K值一定是在數組的后半段。如果K<M，那么K值一定在數組的前半段；（2）根據K值在數組的前半段或后半段來重新分組，再回到第一步進行搜索。分組一直細分到不能再細分為止，而此時若還沒有找到K，代表K值不存在。二元搜索法充分利用數組排序的特性，相比順序搜索法可以減少搜索次數。二元搜索法的程序如下：! 折半查找法范例! By Perng 1997/8/31program BINARY_SE