本文給出一個使用Fortran語言的計算大數階乘的程式，該程式可以計算出1-21萬之間的數的階乘。Fortran表示公式翻譯語言，是最古老的高階語言，主要用途是科學計算，曾經很流行，但現在用的少了，他的大部分市場被C語言取代。在編寫本程式之前，我從未寫過Fortran程式，這個程式是現學現編。為了方便大家學習，我給出這個程式的一些註釋。

　1. Fortran語言的一些語法特點，主要和C語言比。

       1.1 運算子主要有+,-,*,/和**,分別表示加，減，乘，除和乘方，他沒有求餘數的運算子，求餘數使用函式MOD來做。

       1.2 陣列的訪問形式是array(i)，他使用圓括號，而c語言使用方括號。

       1.3 陣列可以定義下標的範圍，下標可以為負數，預設的下標從1開始。而c語言的下標總是從0開始。

       1.4 變數定義的特點，INTEGER類似於c中的int, 變數定義要寫成“變數型別::變數”的形式,不可改變的變數的定義形式為“變數型別,PARAMETER::變數名=”，動態陣列

           申明的格式為“變數型別,DIMENSION(:),ALLOCATABLE::變數名”

       1.5 語句後面不需要任何符號，c語言中，語句後面必須使用分號。

       1.6 動態陣列的分配使用ALLOCATE，釋放使用DEALLOCATE,不需要使用指標。

       1.7 輸入使用READ,輸出使用WRITE。WRITE支援格式符，在語言中，格式符中的特殊格式開始於“%”，可以在格式符中使用字串，而在Fortran中，格式符需要用（）括起來，格式符中的字串必須使用“”，如 '(I6"!="I4.1\)')。在C語言中，輸出預設是不換行,Fortran恰好相反。如果要換行，需在格式符最後加上"\". ,另外，也可以用advance="no"的語法來說明，不換行。

      1.8  c是強型別語言，變數必須定義，然後才可以使用。預設的，Fortran中的變數可以不定義而使用，語句“IMPLICIT NONE”表示變數必須定義。

      1.9  c的行註釋是//,Fortran是 ！

      1.10 Fortran程式中，程式，迴圈塊必須配對，程式以PROGRAM開始，以END PROGRAM結束。

      1.11 Fortran中，關鍵字和變數是大小不敏感的，Fortran90推薦使用小寫格式。這裡則保持了大小格式。

  2. 主要演算法描述

     2.1  階乘的結果是一個很大的數，這裡的大數用陣列來表示，大數使用10000進位制儲存，低位在前，高位再後，和書寫數序相反。buff[1]表示大數的最低4位，buff[2]表示次低4位，以此類推。

    2.2  n!的有多少為10進位制數。

    然後計算N！萬進製表示的長度。因為n！小於(n+1)/2的n次方，故首先使用函式log10,估算(n+1)/2)^n的10進位制的位數，然後除以4+1得到萬進製表示的位數的一個上限。n！小於(n+1)/2的n次方的證明見下：

Caes 1:n是奇數
    n!可表示為連續n個數的乘積，中間的那個數mid=(n+1)/2, 除了這個數外，我們可以將1到n之間的數分成n/2組，每組的兩個數為 mid-i和mid+i (i=1到mid-1),如1,2,3,4,5,6,7，中間的那個數是4，其餘的3對數是(3,5),(2,6)和(1,7),容易知道，每對數的積都於小mid*mid，故n!小於mid的n次方，即 n! <((n+1)/2)^n。
 
Case 2: n 是個偶數，
    n!可表示為連續n個數的乘積，則中間的兩個數(n－1)/2和(n+1)/2, 我們將(n+1)/2記做mid,則中間的一對數是(mid-1)和mid,其它的幾對數是(mid-2,mid+1),(mid-3)(mid+2)等等，故n!小於mid 的n次方,n! <((n+1)/2)^n。

   2.3 具體計算過程就不說了，請看《大數階乘之計算從入門到精通》系列中的入門篇。

3. Fortran開發環境的搭建

   3.1 整合開發環境

     在Windows, 我使用微軟Development Studio的早期版本4.0,這個版本包含了VC4.0 和一個Fortran語言的整合開發環境Fortran PowerStation4.0，是一個32位的編譯器。我下載的是http://download.csdn.net/download/nethouyi/1897519，而我主要參考的文件來自這本數，主要的參考文件則來源於http://download.csdn.net/download/goudanbaoma/2134239。在編寫這個程式的過程中，遇到一些困難，最終是用百度來搜尋得到的。 這個整合開發環境和VC6類似，關於在整合開發環境下編譯和連結的方法，這裡就不多說了。

   3.2 命令列編譯。編譯使用 C:\MSDEV\BIN\ML32.exe,語法和微軟的編譯器cl.exe 差不多。倒是link.exe幾乎會使人瘋掉。這個連結其的版本3.00.5270,不支援在命令列輸入libpath來指定庫檔案的路徑，用這個程式編譯的時候，你需要顯式指定一大堆lib檔案，最終，我不的不放棄這個連結器，而使用了一個更高的版本。下面是我的命令列編譯過程。

        3.2.1 首先執行下面的行，將編譯器和連結器的路徑加入到搜尋路徑。也可以建立批處理檔案，方便以後使用。c:\masm32\bin下包括了一個版本更高的連結器，它是Microsoft (R) Incremental Linker Version 5.12.8078. C:\MSDEV\BIN目錄下包含了fortran語言的編譯器和連結器。因為路徑 C:\masm32\bin 優先於C:\MSDEV\BIN，故如果你這行了這個批出裡，當執行link的時候，他總是使用 C:\masm32\bin下的link.

        set PATH=C:\masm32\bin;C:\MSDEV\BIN;%PATH%;

    3.2.2 編譯和連結程式命令，我的源程式為fac.f90

        fl32 /c fac.f90 

        link fac.obj /machine:i386 /subsystem:console /libpath="c:\msdev\lib" /out:"fac.exe"

 下面是完整的原始碼

PROGRAM PROG1
  IMPLICIT NONE   
  !Calc fraction
  INTEGER::I,J,N,C,PROD  !C:進位,PROD:乘積
  INTEGER::LEN,BUFF_LEN
  INTEGER,PARAMETER::RAD=10000  !使用萬進位制
  INTEGER,DIMENSION(:),ALLOCATABLE::BUFF
  
  WRITE(*,'(" n!="\)')     !不知道"後面的第一個字元為什麼不顯示，故"後跟一個空格
  READ *,N  BUFF_LEN=INT(REAL(N) * LOG10( (REAL(N)+1.0)/2.0)) ! n!的10進位制位數的上限
  BUFF_LEN=BUFF_LEN/4+1     ! n!的萬進位制位數的上限
  
  ALLOCATE(BUFF(BUFF_LEN))

  LEN=1
  BUFF(1)=1

  DO I=1,N
    C=0
    DO J=1,LEN
      PROD=BUFF(J)*I+C
      BUFF(J)=MOD(PROD,RAD)
      C=PROD/RAD 
    END DO
     
    DO
      IF (C==0) EXIT
      BUFF(LEN+1)=MOD(C,RAD)
      C=C/RAD
      LEN=LEN+1
    END DO
      
  END DO
  
  Write(*,'(I6"!="I4.1\)') N,BUFF(LEN)  !\表示不換行
  DO I=LEN-1,1,-1
    WRITE(*,"(I4.4)",advance="no") BUFF(I)
  END DO
  PRINT *
  DEALLOCATE(BUFF)
END PROGRAM PROG1