#include <stdio.h>
#include <math.h>
// Global Value 
extern double BASE ;
extern int NORID ;
extern FILE *FT1 ;
// Function definition
double MTCFM(double *A, double *B, double *C, int N, double *WK) ;
void CFMT4(double *A, int N, int ID, double *WK) ;
void CFMT(double *A, int N, int ID, double *WK) ;
void CFMTTB(int N2, int ID, double *TBL) ;
void CFMTSB(double *A, int M, int L, double *TB, double *WK) ;
void FADD(double *A, double *B, double *C, int N, int ID) ;
void FCADD(int IA, double *B, double *C, int N, int ID) ;
void FCMT(double *A, int IB, double *C, int N) ;
void FCDIV(double *A, int IB, double *C, int N) ;
void NORM(double *A, int N) ;
void FNORM(double *A, int N) ;
void COPY(double *A, double *C, int N) ;
void SETV(int IV, double *C, int N) ;
int LOG2F(int N) ;
double SIGN(double a, double b) ;
double fmax(double A1, double A2) ;
double fmax3(double A1, double A2, double A3) ;
int MIN(int N1, int N2) ;
//==================================================================C
double FMTCFM(double *A, double *B, double *C, int N, double *WK)
//------------------------------------------------------------------C
//   Floating Multiplication (C=A*B) by N-Length Complex FMT        C
//       A(0:N)    R*8, In,  First Input for A                      C
//       B(0:N)    R*8, In,  Second Input for B                     C
//       C(0:2*N)  R*8, Out, Ouput for C=A*B                        C
//       N         I*4, In,  Number of Input(A,B) Elements          C
//       ERR       R*8, Out, Maximum Rounding Error to Convert Int. C
//       WK(2,N,3) R*8, Work, Work for FMT                          C
//------------------------------------------------------------------C
//  Written by Y. Ushiro ( Waseda Uni. and Hitachi )  2003.01.12    C
//==================================================================C
{  double ERR ;
//   Exponent Part
   C[0] = A[0] + B[0] ;
//   Fixed Multiplication
   ERR = MTCFM(&A[1], &B[1], &C[1], N, WK) ; 
//   fprintf(FT1,"C=A*B; %f %f %f %f \n",C[0],C[1],C[2],C[3]);
//   Floating Normalize
   FNORM(C, 2*N) ;
   return (ERR) ;  }
//==================================================================C
double MTCFM(double *A, double *B, double *C, int N, double *WK)
//------------------------------------------------------------------C
//   Fixed Multiplication (C=A*B) by N-Length Complex FMT           C
//       A(N)      R*8, In,  First Input for A                      C
//       B(N)      R*8, In,  Second Input for B                     C
//       C(2*N)    R*8, Out, Ouput for C=A*B                        C
//       N         I*4, In,  Number of Input(A,B) Elements          C
//       ERR       R*8, Out, Maximum Rounding Error to Convert Int. C
//       WK(2,N,3) R*8, Work, Work for FMT                          C
//------------------------------------------------------------------C
//  Written by Y. Ushiro ( Waseda Uni. and Hitachi )  2003.01.12    C
//==================================================================C
{  int    k ;
   double ERR, WR, WI, CR, CI, D5 ;
//  C=A*B using N-Length Complex FMT
//   Move from A,B to WK(1,k), C(2k-1) and Set 0 to WK(2,k), C(2k)
   for ( k=0; k<N; k++)
    { WK[2*k]   = A[k] ;
      WK[2*k+1] = 0.0 ;
      C[2*k]    = B[k] ;
      C[2*k+1]  = 0.0 ; }
//   WK=FMT(A), C=FMT(B)
   CFMT4(WK, N, 0, &WK[2*N]) ;
   CFMT4(C,  N, 0, &WK[2*N]) ;
//   W(i)=W(i)*C(i) with Complex
   for ( k=0; k<N; k++)
    { WR        = WK[2*k] ;
      WI        = WK[2*k+1] ;
      WK[2*k]   = C[2*k]*WR - C[2*k+1]*WI ;
      WK[2*k+1] = C[2*k]*WI + C[2*k+1]*WR ; }
//   WK=FMT**(-1)(WK)
   CFMT4(WK, N, 1, &WK[2*N]) ;
//   fprintf(FT1,"WK; %f %f %f %f\n",WK[0],WK[1],WK[2],WK[3]) ;
//   C=INT(WK) and Change Order of (Real,Imaginary)
   ERR = 0.0 ;
   D5  = 0.5 ;
   for ( k=0; k<N; k++)
    { CR     = floor( WK[2*k]   + D5 ) ;
      CI     = floor( WK[2*k+1] + D5 ) ;
      ERR    = fmax3( ERR, fabs(CR-WK[2*k]), fabs(CI-WK[2*k+1]) ) ;
      C[k]   = CI ;
      C[k+N] = CR ; }
//   C=Normalize(C)
//   fprintf(FT1,"C; %f %f %f %f\n",C[0],C[1],C[2],C[3]) ;
   NORM(C, N*2) ;
   return (ERR) ;  }
//==================================================================C
double FINV(double *A, double *C, int N, double *FW, double *WK)
//------------------------------------------------------------------C
//   Floating Inverse (C=1/A) by N-Length Complex FMT and Mult.     C
//     C <-- X = X + X*(1-A*X) : Iteration for Conversion           C
//       A(0:N)      R*8, In,  First Input for A                    C
//       C(0:N)      R*8, Out, Ouput for C=1/A                      C
//       N           I*4, In,  Number of Elements (A,C)             C
//       ERR         R*8, Out, Maximum Rounding Error to Int.       C
//       FW(0:2*N,2) R*8, Work, Work for Inverse                    C
//       WK(2,N,3)   R*8, Work, Work for FMT                        C
//------------------------------------------------------------------C
//  Written by Y. Ushiro ( Waseda Uni. and Hitachi )  2003.01.12    C
//==================================================================C
{  int    k, NL, NN, NLP, NFP ; 
   double EXP, X, XI, ERR, ERR1, ERR2 ;
//  Set Initial Value
//    Exponent part
    EXP  = A[0] ;
    C[0] = 1.0 - EXP ;
//    Approximation
    X    = BASE/(A[1] + A[2]/BASE) ;
    XI   = (int)(X) ;
    C[1] = XI ;
    C[2] = (int)((X - XI)*BASE) ;
    SETV(0, &C[3], N-2) ;
//  Iteration Loop
    NLP = LOG2F(N) ;
    NL  = 2 ;
    NFP = 2*N + 1 ;
    ERR = 0.0 ;
    for (k=1; k<=NLP+1; k++)
//   FW(1) = A*C,  FW(2) = 1 - FW(1)
     { NN = MIN(N,2*NL) ;
       ERR1 = FMTCFM(A, C, FW, NL, WK) ;
       FCADD(1, FW, &FW[NFP], 2*NL, -1) ;
//   FW(1) = C*FW(2),  C = C + FW(1)
       ERR2 = FMTCFM(C, &FW[NFP], FW, NL, WK) ;
       ERR  = fmax3(ERR ,ERR1, ERR2) ;
       FADD(C, FW, C, NN, 1) ;
       if(NL < N) { NL = 2*NL ; }   }
    return (ERR) ;  }
//==================================================================C
double FDIV(double *A, double *B, double *C, int N, double *FW,
		  double *WK)
//------------------------------------------------------------------C
//   Floating Division (C=A/B) by N-Length Complex FMT and Mult.    C
//       A(0:N)      R*8, In,  First Input for A                    C
//       B(0:N)      R*8, In,  Second Input for B                   C
//       C(0:N)      R*8, Out, Ouput for C=A/B                      C
//       N           I*4, In,  Number of Elements  (A,B,C)          C
//       ERR         R*8, Out, Maximum Rounding Error to  Int.      C
//       FW(0:2*N,2) R*8, Work, Work for Inverse                    C
//       WK(2,N,3)   R*8, Work, Work for FMT                        C
//------------------------------------------------------------------C
//  Written by Y. Ushiro ( Waseda Uni. and Hitachi )  2003.01.12    C
//==================================================================C
{  double ERR,ERR1, ERR2 ;
//   C=1/B
    ERR1 = FINV(B, C, N, FW, WK) ;
//   C=FW=A/B 
    ERR2 = FMTCFM(A, C, FW, N, WK) ;
    COPY(FW, C, N) ;
    ERR = fmax(ERR1,ERR2) ;
    return (ERR) ; }
//==================================================================C
double FCSQRT(int IA, double *C, int N, int ID,	double *FW,
			  double *WK)
//------------------------------------------------------------------C
//   Floating Constant SQRT (C=SQRT(IA)) by Complex FMT and Mult.   C
//     C = 1/SQRT(IA) <-- X = X + X*(1-IA*X*X)/2 : Iteration        C
//       IA          I*4, In,  Constant Input                       C
//       C(0:N)      R*8, Out, Ouput for C=SQRT(IA) or 1/SQRT(IA)   C
//       N           I*4, In,  Number of Output Elements            C
//       ID          I*4, In,  ID= 1 ; C=SQRT(IA)                   C
//                                -1 ; C=1/SQRT(IA)                 C
//       ERR         R*8, Out, Maximum Rounding Error to Int.       C
//       FW(0:2*N,2) R*8, Work, Work for SQRT                       C
//       WK(2,N,3)   R*8, Work, Work for FMT                        C
//------------------------------------------------------------------C
//  Written by Y. Ushiro ( Waseda Uni. and Hitachi )  2003.01.12    C
//==================================================================C
{  int k, NL, NN, NLP, NFP ;
   double X, XI, ERR, ERR1, ERR2 ;
//  Set Initial Value
   C[0] = 0.0 ;
   X    = BASE/sqrt(IA*1.0) ;
   XI   = (int)(X) ;
   C[1] = XI ;
   C[2] = (int)((X - XI)*BASE) ;
   SETV(0, &C[3], N-2) ;
//  Iteration Number
   NLP = LOG2F(N) ;
   NFP = 2*N + 1 ;
   NL  = 2 ;
   ERR = 0.0 ;
   for (k=1; k<=NLP+1; k++)
//   FW(1) = C*C,  FW(2) = IA*FW(1),  FW(1) = 1 - FW(2)
    { NN = MIN(N, 2*NL) ;
      ERR1 = FMTCFM(C, C, FW, NL, WK) ;
      FCMT(FW, IA, &FW[NFP], NN) ;
      FCADD(1, &FW[NFP], FW, NL*2, -1) ;
//   FW(2) = C*FW(1),  FW(1) = FW(2)/2,  C = C + FW(1)
      ERR2 = FMTCFM(C, FW, &FW[NFP], NL, WK) ;
      FCDIV(&FW[NFP], 2, FW, NN) ;
      ERR = fmax3(ERR, ERR1, ERR2) ;
      FADD(C, FW, C, NN, 1) ;
      if (NL < N) { NL = 2*NL ; }    }
   if(ID >= 0) { 
//   C=SQRT(IA)=IA*(1/SQRT(IA)),  C=FW=IA*C
      FCMT(C, IA, FW, N) ;
      COPY(FW, C, N+1) ;   }
    return (ERR) ; }
//==================================================================C
double FSQRT(double *A, double *C, int N,int ID, double *FW,double *WK)
//------------------------------------------------------------------C
//   Floating SQRT (C=SQRT(A)) by Complex FMT and Mult.             C
//     C = 1/SQRT(A) <-- X = X + X*(1-IA*X*X)/2 : Iteration         C
//       A(0:N)      R*8, In,  Input Data A                         C
//       C(0:N)      R*8, Out, Ouput for C=SQRT(A) or 1/SQRT(A)     C
//       N           I*4, In,  Number of Output Elements            C
//       ID          I*4, In,  ID= 1 ; C=SQRT(A)                    C
//                                -1 ; C=1/SQRT(A)                  C
//       ERR         R*8, Out, Maximum Rounding Error to Int.       C
//       FW(0:2*N,2) R*8, Work, Work for SQRT                       C
//       WK(2,N,3)   R*8, Work, Work for FMT                        C
//------------------------------------------------------------------C
//  Written by Y. Ushiro ( Waseda Uni. and Hitachi )  2003.01.12    C
//==================================================================C
{  int    k, NL, NLP, NN, NFP, IA0 ;
   double OEXP, BS1, X, XI, ERR, ERR1, ERR2, ERR3 ;
//  Set Initial Value
   IA0  = A[0] ;
   OEXP = IA0%2 ;
   if(fabs(OEXP) == 1.0) {
      C[0] = (-A[0] + 1.0)/2.0 ;
      BS1  = BASE ;  }
   else {
      C[0] = -A[0]/2.0 + 1.0 ;
      BS1  = sqrt(BASE) ;   }
   X    = BS1/sqrt(A[1]+A[2]/BASE) ;
   XI   = (int)(X) ;
   C[1] = XI ;
   C[2] = (int)((X - XI)*BASE) ;
   SETV(0, &C[3], N-2) ;
//  Iteration Loop
   NLP = LOG2F(N) ;
   NFP = 2*N + 1 ;
   NL  = 2 ;
   ERR = 0.0 ;
   for (k=1; k<=NLP+1; k++)
//   FW(1) = C*C,  FW(2) = A*FW(1),  FW(1) = 1 - FW(2)
    { NN = MIN(N, 2*NL) ;
      ERR1 = FMTCFM(C, C, FW, NL, WK) ;
      ERR2 = FMTCFM(FW, A, &FW[NFP], NL, WK) ;
      ERR = fmax3(ERR, ERR1, ERR2) ;
      FCADD(1, &FW[NFP], FW, 2*NL, -1) ;
//   FW(2) = C*FW(1),  FW(1) = FW(2)/2,  C = C + FW(1)
      ERR3 = FMTCFM(C, FW, &FW[NFP], NL, WK) ;
      FCDIV(&FW[NFP], 2, FW, NN) ;
      ERR = fmax(ERR, ERR3) ;
      FADD(C, FW, C, NN, 1) ;
      if(NL < N) { NL = 2*NL ; }   }
   if(ID >= 0) { 
//   C=SQRT(A)=A*(1/SQRT(A)),  C=FW=A*C
      ERR1 = FMTCFM(C, A, FW, N, WK) ;
      ERR = fmax(ERR, ERR1) ;
      COPY(FW, C, N+1) ;  }
   return (ERR) ;  }
//==================================================================C
void CFMT4(double *A, int N, int ID, double *WK)
//------------------------------------------------------------------C
//   A = Complex FMT for Real Convolution                           C
//       A(N)     C*16, I/O,  Input and Output Data for FMT         C
//       N        I*4,  In,   Number of FMT Elements, N=2**NB       C
//       ID       I*4,  In,   ID=0 : Normal FMT                     C
//                            ID=1 : Reverse FMT with 1/N           C
//       WK(N,2)  C*16, Work, Work for FMT                          C
//------------------------------------------------------------------C
//  Written by Y. Ushiro ( Waseda Uni. and Hitachi )  2003.01.12    C
//==================================================================C
{  int    k ;
   double PIN, ARG, AR, AI, OR, OI ; 
//   Set Constant
   PIN = 2.0*atan(1.0)/N ;
//   Test ID
   if(ID != 1) { 
//  Normal FMT : A=FMT(A*OMG(k))
//    A(k)=CEXP(Pi*(k)*i/(2N))
      for (k=0; k<N; k++)
       { ARG = PIN*(N-1 - k) ;
         OR  = cos(ARG) ;   OI = sin(ARG) ;
         AR  = A[2*k] ;     AI = A[2*k+1] ;
	     A[2*k]   = OR*AR - OI*AI ;
         A[2*k+1] = OR*AI + OI*AR ;	  }	 
//    A=FMT(A)
      CFMT(A, N, ID, WK) ;  }
    else { 
//  Reverse FMT : A=FMT(A)*OMG(-k)
//    A=FMT(A)
       CFMT(A, N, ID, WK) ;
//    A=A*OMG(-k)
       PIN = - PIN ;
       for (k=0; k<N; k++)
        { ARG    = PIN*(N-1 - k) ;
         OR  = cos(ARG) ;   OI = sin(ARG) ;
         AR  = A[2*k] ;     AI = A[2*k+1] ;
	     A[2*k]   = OR*AR - OI*AI ;
         A[2*k+1] = OR*AI + OI*AR ;	  }  }	 
  }
//==================================================================C
void CFMT(double *A, int N, int ID, double *WK)
//------------------------------------------------------------------C
//   A = FMT(A) or Reverse FMT(A) , Complex FMT                     C
//       A(N)     C*16, I/O,  Input and Output Data for FMT         C
//       N        I*4,  In,   Number of FMT Elements, N=2**NB       C
//       ID       I*4,  In,   ID=0 : Normal FMT                     C
//                            ID=1 : Reverse FMT with 1/N           C
//       WK(N,2)  C*16, Work, Work for FMT                          C
//------------------------------------------------------------------C
//  Written by Y. Ushiro ( Waseda Uni. and Hitachi )  2003.01.12    C
//==================================================================C
{  int    k, NB, L, M ;
//   Set NB ; N=2**NB 
   NB = LOG2F(N) ;
//   Set FMT Table 
   CFMTTB(N/2, ID, WK) ; 
//   Main loop for FMT
   L = 1 ;
   M = N/2 ;
   for (k=0; k<NB; k++)
    { CFMTSB(A, M, L, WK, &WK[2*N]);
      L = L*2 ;
      M = M/2 ;  }
//   Reverse FMT, A=A/N
   if (ID == 1) {
      for (k=0; k<N; k++)
        { A[2*k]   = A[2*k]/N ;
          A[2*k+1] = A[2*k+1]/N ; }  }
   }
//==================================================================C
void CFMTTB(int N2, int ID, double *TBL)
//------------------------------------------------------------------C
//   Set Complex FMT TBL,  TBL(k)=CEXP(Pi*(k-1)*i/N2)               C
//       N2        I*4, In,  Half Number of FMT Elements            C
//       ID        I*4, In,  ID=0 : Normal FMT                      C
//                           ID=1 : Reverse FMT with 1/N            C
//       TBL(2,N2) R*8, In,  Translation Table                      C
//------------------------------------------------------------------C
//  Written by Y. Ushiro ( Waseda Uni. and Hitachi )  2003.01.12    C
//==================================================================C
{  int    k ;
   double PIN, ARG ;
//   TBL(k)=CEXP(Pi*(k-1)*i/N2)
   PIN = 4.0*atan(1.0)/N2 ;
   if (ID == 1) PIN = -PIN ;
   for (k=0; k<N2; k++)
    { ARG        = PIN*(N2-1 - k) ;
      TBL[2*k]   = cos(ARG) ;
      TBL[2*k+1] = sin(ARG) ;  }
  }
//==================================================================C
void CFMTSB(double *A, int M, int L, double *TB, double *WK)
//------------------------------------------------------------------C
//   FMT Complex Translation Subroutine                             C
//       A(M,2,L)  C*16, I/O,  Input and Output Data                C
//       M         I*4,  In,   Inner Size of A                      C
//       L         I*4,  In,   Outer Size of A                      C
//       TB(M,L)   C*16, In,   Translation Table                    C
//       WK(M,L,2) C*16, Work, Work                                 C
//------------------------------------------------------------------C
//  Written by Y. Ushiro ( Waseda Uni. and Hitachi )  2003.01.13    C
//==================================================================C
{  int    i, j, IA, IW, ML2 ;
   double OR, OI, AR, AI, WR, WI ; 
//   FMT Translation
   ML2 = 2*M*L ;
   for (j=0; j<L; j++)
    { OR = TB[2*M*(j+1)-2] ;
	  OI = TB[2*M*(j+1)-1] ;
      for (i=0; i<M; i++)
       { IA = 2*(2*M*j+i) ;  
         IW = 2*(M*j+i) ;
         AR = A[IA] ;  AI = A[IA+1] ;
		 WR = AR*OR - AI*OI ;
		 WI = AR*OI + AI*OR ;
         WK[IW]       = A[IA+2*M]   - WR ;
         WK[IW+1]     = A[IA+2*M+1] - WI ;
         WK[IW+ML2]   = A[IA+2*M]   + WR ;
         WK[IW+ML2+1] = A[IA+2*M+1] + WI ; 
       }  }
//   Copy Back
   for (i=0; i<ML2; i++)
    { A[2*i]   = WK[2*i] ;
      A[2*i+1] = WK[2*i+1] ;  }
  }