多数桁乗除算FORTRANプログラム

C==================================================================C
      SUBROUTINE FMTCFM(A,B,C,N,ERR,WK)
C------------------------------------------------------------------C
C   Floating Multiplication (C=A*B) by N-Length Complex FMT        C
C       A(0:N)    R*8, In,  First Input for A                      C
C       B(0:N)    R*8, In,  Second Input for B                     C
C       C(0:2*N)  R*8, Out, Ouput for C=A*B                        C
C       N         I*4, In,  Number of Input(A,B) Elements          C
C       ERR       R*8, Out, Maximum Rounding Error to Convert Int. C
C       WK(2,N,3) R*8, Work, Work for FMT                          C
C------------------------------------------------------------------C
C  Written by Y. Ushiro ( Waseda Uni. and Hitachi )  2003.01.05    C
C==================================================================C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION A(0:N), B(0:N), C(0:2*N), WK(2,N,3)
C   Exponent Part
      C(0) = A(0) + B(0)
C   Fixed Multiplication
      CALL MTCFM(A(1),B(1),C(1),N,ERR,WK) 
C   Floating Normalize
      CALL FNORM(C,2*N)
C
      RETURN
      END
C==================================================================C
      SUBROUTINE MTCFM(A,B,C,N,ERR,WK)
C------------------------------------------------------------------C
C   Fixed Multiplication (C=A*B) by N-Length Complex FMT           C
C       A(N)      R*8, In,  First Input for A                      C
C       B(N)      R*8, In,  Second Input for B                     C
C       C(2*N)    R*8, Out, Ouput for C=A*B                        C
C       N         I*4, In,  Number of Input(A,B) Elements          C
C       ERR       R*8, Out, Maximum Rounding Error to Convert Int. C
C       WK(2,N,3) R*8, Work, Work for FMT                          C
C------------------------------------------------------------------C
C  Written by Y. Ushiro ( Waseda Uni. and Hitachi )  2003.01.05    C
C==================================================================C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION A(N), B(N), C(2*N), WK(2,N,3)
C  C=A*B using N-Length Complex FMT
C   Move from A,B to WK(1,k), C(2k-1) and Set 0 to WK(2,k), C(2k)
      DO k=1,N
        WK(1,k,1) = A(k)
        WK(2,k,1) = 0.0
        C(2*k-1)  = B(k)
        C(2*k)    = 0.0
      END DO
C   WK=FMT(A), C=FMT(B)
      CALL CFMT4(WK,N,0,WK(1,1,2))
      CALL CFMT4(C, N,0,WK(1,1,2))
C   W(i)=W(i)*C(i) with Complex
      DO k=1,N
        WR        = WK(1,k,1)
        WC        = WK(2,k,1)
        WK(1,k,1) = C(2*k-1)*WR - C(2*k)*WC
        WK(2,k,1) = C(2*k-1)*WC + C(2*k)*WR
      END DO
C   WK=FMT**(-1)(WK)
      CALL CFMT4(WK,N,1,WK(1,1,2))
C   C=INT(WK) and Change Order of (Real,Imaginary)
      ERR = 0.0
      D5  = 0.5
      DO k=1,N
        CR     = DINT( WK(1,k,1) + DSIGN(D5,WK(1,k,1)) )
        CI     = DINT( WK(2,k,1) + DSIGN(D5,WK(2,k,1)) )
        ERR    = MAX( ERR, ABS(CR-WK(1,k,1)), ABS(CI-WK(2,k,1)) )
        C(k)   = CI
        C(k+N) = CR
      END DO
C   C=Normalize(C)
      CALL NORM(C,N*2)
C
      RETURN
      END
C==================================================================C
      SUBROUTINE FINV(A,C,N,ERR,FW,WK)
C------------------------------------------------------------------C
C   Floating Inverse (C=1/A) by N-Length Complex FMT and Mult.     C
C     C <-- X = X + X*(1-A*X) : Iteration for Conversion           C
C       A(0:N)      R*8, In,  First Input for A                    C
C       C(0:N)      R*8, Out, Ouput for C=1/A                      C
C       N           I*4, In,  Number of Elements (A,C)             C
C       ERR         R*8, Out, Maximum Rounding Error to Int.       C
C       FW(0:2*N,2) R*8, Work, Work for Inverse                    C
C       WK(2,N,3)   R*8, Work, Work for FMT                        C
C------------------------------------------------------------------C
C  Written by Y. Ushiro ( Waseda Uni. and Hitachi )  2003.01.05    C
C==================================================================C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION A(0:N), C(0:N), FW(0:2*N,2), WK(2,N,3)
      COMMON /VBASE/BASE,NORID
C  Set Initial Value
C    Exponent part
      EXP  = A(0)
      C(0) = 1.0 - EXP
C    Approximation
      X    = BASE/(A(1) + A(2)/BASE)
      XI   = DINT(X)
      C(1) = XI
      C(2) = DINT((X - XI)*BASE)
      CALL SETV(0,C(3),N-2)
C  Iteration Loop
      CALL LOG2F(N,NLP)
      NL  = 2
      ERR = 0.0
      DO k=1,NLP+1
C   FW(1) = A*C,  FW(2) = 1 - FW(1)
        NN = MIN(N,2*NL)
        CALL FMTCFM(A,C,FW(0,1),NL,ERR1,WK)
        CALL FCADD(1,FW(0,1),FW(0,2),2*NL,-1)
C   FW(1) = C*FW(2),  C = C + FW(1)
        CALL FMTCFM(C,FW(0,2),FW(0,1),NL,ERR2,WK)
        ERR = MAX(ERR,ERR1,ERR2)
        CALL FADD(C,FW(0,1),C,NN,1)
        IF(NL.LT.N) NL = 2*NL
      END DO 
C  
      RETURN
      END
C==================================================================C
      SUBROUTINE FDIV(A,B,C,N,ERR,FW,WK)
C------------------------------------------------------------------C
C   Floating Division (C=A/B) by N-Length Complex FMT and Mult.    C
C       A(0:N)      R*8, In,  First Input for A                    C
C       B(0:N)      R*8, In,  Second Input for B                   C
C       C(0:N)      R*8, Out, Ouput for C=A/B                      C
C       N           I*4, In,  Number of Elements  (A,B,C)          C
C       ERR         R*8, Out, Maximum Rounding Error to  Int.      C
C       FW(0:2*N,2) R*8, Work, Work for Inverse                    C
C       WK(2,N,3)   R*8, Work, Work for FMT                        C
C------------------------------------------------------------------C
C  Written by Y. Ushiro ( Waseda Uni. and Hitachi )  2003.01.05    C
C==================================================================C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION A(0:N), B(0:N), C(0:N), FW(0:2*N,2), WK(2,N,3)
C   C=1/B
      CALL FINV(B,C,N,ERR1,FW,WK)
C   C=FW=A/B 
      CALL FMTCFM(A,C,FW,N,ERR2,WK)
      CALL COPY(FW,C,N)
      ERR = MAX(ERR1,ERR2)
C
      RETURN
      END
C==================================================================C
      SUBROUTINE FCSQRT(IA,C,N,ID,ERR,FW,WK)
C------------------------------------------------------------------C
C   Floating Constant SQRT (C=SQRT(IA)) by Complex FMT and Mult.   C
C     C = 1/SQRT(IA) <-- X = X + X*(1-IA*X*X)/2 : Iteration        C
C       IA          I*4, In,  Constant Input                       C
C       C(0:N)      R*8, Out, Ouput for C=SQRT(IA) or 1/SQRT(IA)   C
C       N           I*4, In,  Number of Output Elements            C
C       ID          I*4, In,  ID= 1 ; C=SQRT(IA)                   C
C                                -1 ; C=1/SQRT(IA)                 C
C       ERR         R*8, Out, Maximum Rounding Error to Int.       C
C       FW(0:2*N,2) R*8, Work, Work for SQRT                       C
C       WK(2,N,3)   R*8, Work, Work for FMT                        C
C------------------------------------------------------------------C
C  Written by Y. Ushiro ( Waseda Uni. and Hitachi )  2003.01.05    C
C==================================================================C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION C(0:N), FW(0:2*N,2), WK(2,N,3)
      COMMON /VBASE/BASE,NORID
C  Set Initial Value
      C(0) = 0.0
      X    = BASE/DSQRT(DFLOAT(IA))
      XI   = DINT(X)
      C(1) = XI
      C(2) = DINT((X - XI)*BASE)
      CALL SETV(0,C(3),N-2)
C  Iteration Number
      CALL LOG2F(N,NLP)
      NL  = 2
      ERR = 0.0
      DO k=1,NLP+1
C   FW(1) = C*C,  FW(2) = IA*FW(1),  FW(1) = 1 - FW(2)
        NN = MIN(N,2*NL)
        CALL FMTCFM(C,C,FW(0,1),NL,ERR1,WK)
        CALL FCMT(FW(0,1),IA,FW(0,2),NN)
        CALL FCADD(1,FW(0,2),FW(0,1),NL*2,-1)
C   FW(2) = C*FW(1),  FW(1) = FW(2)/2,  C = C + FW(1)
        CALL FMTCFM(C,FW(0,1),FW(0,2),NL,ERR2,WK)
        CALL FCDIV(FW(0,2),2,FW(0,1),NN)
        ERR = MAX(ERR,ERR1,ERR2)
        CALL FADD(C,FW(0,1),C,NN,1)
        IF(NL.LT.N) NL = 2*NL
      END DO 
      IF(ID.GE.0) THEN 
C   C=SQRT(IA)=IA*(1/SQRT(IA)),  C=FW=IA*C
        CALL FCMT(C,IA,FW,N)
        CALL COPY(FW,C,N+1)
      END IF
C
      RETURN
      END
C==================================================================C
      SUBROUTINE FSQRT(A,C,N,ID,ERR,FW,WK)
C------------------------------------------------------------------C
C   Floating SQRT (C=SQRT(A)) by Complex FMT and Mult.             C
C     C = 1/SQRT(A) <-- X = X + X*(1-IA*X*X)/2 : Iteration         C
C       A(0:N)      R*8, In,  Input Data A                         C
C       C(0:N)      R*8, Out, Ouput for C=SQRT(A) or 1/SQRT(A)     C
C       N           I*4, In,  Number of Output Elements            C
C       ID          I*4, In,  ID= 1 ; C=SQRT(A)                    C
C                                -1 ; C=1/SQRT(A)                  C
C       ERR         R*8, Out, Maximum Rounding Error to Int.       C
C       FW(0:2*N,2) R*8, Work, Work for SQRT                       C
C       WK(2,N,3)   R*8, Work, Work for FMT                        C
C------------------------------------------------------------------C
C  Written by Y. Ushiro ( Waseda Uni. and Hitachi )  2003.01.05    C
C==================================================================C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION A(0:N), C(0:N), FW(0:2*N,2), WK(2,N,3)
      COMMON /VBASE/BASE,NORID
C  Set Initial Value
      OEXP = MOD(A(0),2.0D0)
      IF(ABS(OEXP).EQ.1.0D0) THEN
        C(0) = (-A(0)+1.0)/2.0 
        BS1  = BASE
      ELSE
        C(0) = -A(0)/2.0 + 1.0
        BS1  = SQRT(BASE)
      END IF 
      X    = BS1/DSQRT(A(1)+A(2)/BASE)
      XI   = DINT(X)
      C(1) = XI
      C(2) = DINT((X - XI)*BASE)
      CALL SETV(0,C(3),N-2)
C  Iteration Loop
      CALL LOG2F(N,NLP)
      NL  = 2
      ERR = 0.0
      DO k=1,NLP+1
C   FW(1) = C*C,  FW(2) = A*FW(1),  FW(1) = 1 - FW(2)
        NN = MIN(N,2*NL)
        CALL FMTCFM(C,C,FW(0,1),NL,ERR1,WK)
        CALL FMTCFM(FW(0,1),A,FW(0,2),NL,ERR2,WK)
        CALL FCADD(1,FW(0,2),FW(0,1),2*NL,-1)
C   FW(2) = C*FW(1),  FW(1) = FW(2)/2,  C = C + FW(1)
        CALL FMTCFM(C,FW(0,1),FW(0,2),NL,ERR3,WK)
        CALL FCDIV(FW(0,2),2,FW(0,1),NN)
        ERR = MAX(ERR,ERR1,ERR2,ERR3)
        CALL FADD(C,FW(0,1),C,NN,1)
        IF(NL.LT.N) NL = 2*NL
      END DO 
      IF(ID.GE.0) THEN 
C   C=SQRT(A)=A*(1/SQRT(A)),  C=FW=A*C
        CALL FMTCFM(C,A,FW,N,ERR1,WK)
        ERR = MAX(ERR,ERR1)
        CALL COPY(FW,C,N+1)
      END IF
C
      RETURN
      END
C==================================================================C
      SUBROUTINE CFMT4(A,N,ID,WK)
C------------------------------------------------------------------C
C   A = Complex FMT for Real Convolution                           C
C       A(N)     C*16, I/O,  Input and Output Data for FMT         C
C       N        I*4,  In,   Number of FMT Elements, N=2**NB       C
C       ID       I*4,  In,   ID=0 : Normal FMT                     C
C                            ID=1 : Reverse FMT with 1/N           C
C       WK(N,2)  C*16, Work, Work for FMT                          C
C------------------------------------------------------------------C
C  Written by Y. Ushiro ( Waseda Uni. and Hitachi )  2002.10.29    C
C==================================================================C
      IMPLICIT REAL*8 (A-H,O-Z)
      COMPLEX*16 A(N), WK(N,2), OMG
C   Set Constant
      PIN = 2.0*DATAN(1.0D0)/N
C   Test ID
      IF(ID.NE.1) THEN
C  Normal FMT : A=FMT(A*OMG(k))
C    A(k)=CEXP(Pi*(k-1)*i/(2N))
        DO k=1,N
          ARG  = PIN*(N - k)
          OMG  = DCMPLX(DCOS(ARG),DSIN(ARG))
          A(k) = A(k)*OMG
        END DO
C    A=FMT(A)
        CALL CFMT(A,N,ID,WK)
      ELSE
C  Reverse FMT : A=FMT(A)*OMG(-k)
C    A=FMT(A)
        CALL CFMT(A,N,ID,WK) 
C    A=A*OMG(-k)
        PIN = - PIN
        DO k=1,N
          ARG  = PIN*(N - k)
          OMG  = DCMPLX(DCOS(ARG),DSIN(ARG))
          A(k) = A(k)*OMG
        END DO
      END IF
C
      RETURN
      END
C==================================================================C
      SUBROUTINE CFMT(A,N,ID,WK)
C------------------------------------------------------------------C
C   A = FMT(A) or Reverse FMT(A) , Complex FMT                     C
C       A(N)     C*16, I/O,  Input and Output Data for FMT         C
C       N        I*4,  In,   Number of FMT Elements, N=2**NB       C
C       ID       I*4,  In,   ID=0 : Normal FMT                     C
C                            ID=1 : Reverse FMT with 1/N           C
C       WK(N,2)  C*16, Work, Work for FMT                          C
C------------------------------------------------------------------C
C  Written by Y. Ushiro ( Waseda Uni. and Hitachi )  2002.10.29    C
C==================================================================C
      IMPLICIT COMPLEX*16(A-H,O-Z)
      DIMENSION A(N), WK(N,2)
C   Set NB ; N=2**NB 
      CALL LOG2F(N,NB)
C   Set FMT Table 
      CALL CFMTTB(N/2,ID,WK)
C   Main loop for FMT
      L = 1
      M = N/2
      DO k=1,NB
        CALL CFMTSB(A,M,L,WK(1,1),WK(1,2))
        L = L*2
        M = M/2
      END DO
C   Reverse FMT, A=A/N
      IF(ID.EQ.1) THEN
        DO k=1,N
          A(k) = A(k)/N
        END DO
      END IF
C
      RETURN
      END
C==================================================================C
      SUBROUTINE CFMTTB(N2,ID,TBL)
C------------------------------------------------------------------C
C   Set Complex FMT TBL,  TBL(k)=CEXP(Pi*(k-1)*i/N2)               C
C       N2        I*4, In,  Half Number of FMT Elements            C
C       ID        I*4, In,  ID=0 : Normal FMT                      C
C                           ID=1 : Reverse FMT with 1/N            C
C       TBL(2,N2) R*8, In,  Translation Table                      C
C------------------------------------------------------------------C
C  Written by Y. Ushiro ( Waseda Uni. and Hitachi )  2002.10.29    C
C==================================================================C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION TBL(2,N2)
C   TBL(k)=CEXP(Pi*(k-1)*i/N2)
      PIN = 4.0*DATAN(1.0D0)/N2
      IF(ID.EQ.1) PIN = -PIN
      DO k=1,N2
        ARG      = PIN*(N2 - k)
        TBL(1,k) = DCOS(ARG)
        TBL(2,k) = DSIN(ARG)
      END DO
C
      RETURN
      END
C==================================================================C
      SUBROUTINE CFMTSB(A,M,L,TB,WK)
C------------------------------------------------------------------C
C   FMT Complex Translation Subroutine                             C
C       A(M,2,L)  C*16, I/O,  Input and Output Data                C
C       M         I*4,  In,   Inner Size of A                      C
C       L         I*4,  In,   Outer Size of A                      C
C       TB(M,L)   C*16, In,   Translation Table                    C
C       WK(M,L,2) C*16, Work, Work                                 C
C------------------------------------------------------------------C
C  Written by Y. Ushiro ( Waseda Uni. and Hitachi )  2002.10.29    C
C==================================================================C
      IMPLICIT COMPLEX*16(A-H,O-Z)
      DIMENSION A(M,2,L), TB(M,L), WK(M,L,2)
C   FMT Translation
      DO j=1,L
        OMG = TB(M,j)
        DO i=1,M
          A1 = A(i,1,j)*OMG
          WK(i,j,1) = A(i,2,j) - A1
          WK(i,j,2) = A(i,2,j) + A1
        END DO
      END DO
C   Copy Back
      DO i=1,M*L*2
        A(i,1,1) = WK(i,1,1)
      END DO
C
      RETURN
      END