#include <stdio.h>
#include <math.h>
// Global Value 
extern double BASE ;
extern int NORID;
// Function
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) ;
//==================================================================C
void ADD(double *A, double *B, double *C, int N, int ID)
//------------------------------------------------------------------C
//   Fixed Add or Subtract (C=A+B or A-B) with Normalize by BASE    C
//       A(N)     R*8, In,  First Input for A                       C
//       B(N)     R*8, In,  Second Input for B                      C
//       C(N)     R*8, Out, Ouput for C=A+B or A-B                  C
//       N        I*4, In,  Number of Input(A,B) Elements           C
//       ID       I*4, In,  ID= 1 : C=A+B                           C
//                            =-1 : C=A-B                           C
//------------------------------------------------------------------C
//  Written by Y. Ushiro ( Waseda Uni. and Hitachi )  2003.01.12    C
//==================================================================C
{  int i ; 
   if(ID >= 0) {
      for (i=0; i<N; i++) 
	    { C[i] = A[i] + B[i] ; }  }      //  C=A+B
   else {
      for (i=0; i<N; i++)
		{ C[i] = A[i] - B[i] ; } }       //  C=A-B
   NORM(C,N) ;  }                        //  Fixed Normalize
//==================================================================C
void FADD(double *A, double *B, double *C, int N, int ID)
//------------------------------------------------------------------C
//   Floating Add or Subtract (C=A+B or A-B) with Normalize         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 or A-B                  C
//       N        I*4, In,  Number of Input(A,B) Elements           C
//       ID       I*4, In,  ID= 1 : C=A+B                           C
//                            =-1 : C=A-B                           C
//------------------------------------------------------------------C
//  Written by Y. Ushiro ( Waseda Uni. and Hitachi )  2003.01.12    C
//==================================================================C
{ int  i, K ;
//   Check Exponent
  if (A[0] >= B[0]) {
//   Exponent A >= B
     K = A[0] - B[0] ;
     COPY(A, C, K+1) ;
//    Fixed ADD or SUB
     ADD(&A[K+1], &B[1], &C[K+1], N-K, ID) ;
      }
  else {
     C[0] = B[0] ;
//   Exponent A <  B
     K = B[0] - A[0] ;
     if (ID >= 0) {
        COPY(&B[1], &C[1], K) ; }
     else {
        for (i=1; i<=K; i++)
		  { C[i] = -B[i] ; } }
//    Fixed ADD or SUB
     ADD(&A[1], &B[K+1], &C[K+1], N-K, ID) ;
    }
//   Normalize
  NORM(&C[1], N) ;
  FNORM(C, N) ;  }
//==================================================================C
void FCADD(int IA, double *B, double *C, int N, int ID)
//------------------------------------------------------------------C
//   Floating Constant Add or Subtract (C=IA-B)                     C
//       IA       I*4, In,  Constant Value                          C
//       B(0:N)   R*8, In,  Second Input for B                      C
//       C(0:N)   R*8, Out, Ouput for C=IA+B or IA-B                C
//       N        I*4, In,  Number of B,C Elements                  C
//       ID       I*4, In,  ID= 1 : C=IA+B                          C
//                            =-1 : C=IA-B                          C
//------------------------------------------------------------------C
//  Written by Y. Ushiro ( Waseda Uni. and Hitachi )  2003.01.12    C
//==================================================================C
{  int    i, K ;
   double S ;
//  Check Exponent of B
   S = 1.0 ;
   if (ID < 0) S = -1.0 ;
   if (B[0] >= 1.0) {
//  B(0) >= 1
      C[0] = B[0] ;
      K    = B[0] ;
      for (i=1; i<=N; i++)
		{ C[i] = S*B[i] ; }
      C[K] = IA + C[K] ;   }
   else {
//  B(0) <= 0
      C[0] = 1 ;
      K    = 1 - B[0] ;
      C[1] = IA ;
      SETV(0, &C[2], K-1) ;
      for (i=K+1; i<=N; i++)
		{ C[i] = S*B[i-K] ; }  }      
//   Fixed and Floating Normalize
   NORM(&C[1], N) ;
   FNORM(C, N) ;   }
//==================================================================C
void FCMT(double *A, int IB, double *C, int N)
//------------------------------------------------------------------C
//   Floating Constant Multiplication (C=A*IB)                      C
//       A(0:N)   R*8, In,  First Input for A                       C
//       IB       I*4, In,  Constant value  (IB <= BASE)            C
//       C(0:N)   R*8, Out, Ouput for C=A*B                         C
//       N        I*4, In,  Number of A,C Elements                  C
//------------------------------------------------------------------C
//  Written by Y. Ushiro ( Waseda Uni. and Hitachi )  2003.01.12    C
//==================================================================C
{  int    i ;
   double B, TOP ;
//   Top Multiplication
   B    = IB ;
   TOP  = fabs(A[1]*B) ;
   if (TOP >= BASE) {
//   Fixed Multiplication ( TOP >= BASE )
      C[0] = A[0] + 1 ;
      C[1] = 0 ;
      for (i=1; i<N ; i++) 
		{ C[i+1] = A[i]*B ;}
      C[N] = C[N] + (int)(A[N]*B/BASE) ;  }
   else {
//   Fixed Multiplication ( TOP < BASE )
      C[0] = A[0] ;
      for (i=1; i<=N ; i++) 
		{ C[i] = A[i]*B ; }  }
//   Fixed and Floating Normalize
   NORM(&C[1], N) ;
   FNORM(C, N) ;   }
//==================================================================C
void FCDIV(double *A, int IB, double *C, int N)
//------------------------------------------------------------------C
//   Floating Constant Division (C=A/IB)                            C
//       A(0:N)    R*8, In,  First Input for A                      C
//       IB        I*4, In,  Constant Value  (IB <= BASE)           C
//       C(0:N)    R*8, Out, Ouput for C=A/B                        C
//       N         I*4, In,  Number of A,C Elements                 C
//------------------------------------------------------------------C
//  Written by Y. Ushiro ( Waseda Uni. and Hitachi )  2003.01.12    C
//==================================================================C
{  int    i ;
   double B, BV, DV ;
//   Top Check
   B  = IB ;
   if (fabs(A[1]) >= fabs(B)) {
//   Fixed Division ( |A(1)| >= |B| )
      C[0] = A[0] ;
      BV   = 0 ;
      for (i=1; i<=N; i++) {
         DV   = BV*BASE + A[i] ;
         C[i] = (int)(DV/B) ;
         BV   = DV - C[i]*B ; }  }
   else {
//   Fixed Division ( |A(1)| < |B| )
      C[0] = A[0] - 1 ;
      BV   = A[1] ;
      for (i=1; i<N ; i++) {
         DV   = BV*BASE + A[i+1] ;
         C[i] = (int)(DV/B) ;
         BV   = DV - C[i]*B ; }
      C[N] = (int)(BV*BASE/B) ;  }
//   Fixed and Floating Normalize
   NORM(&C[1], N) ;
   FNORM(C, N) ;   }
//==================================================================C
void CMLT(double *A, double X, int N)
//==================================================================C
//   Fixed Multiplication (A(i)=A(i)*X) with Fixed Normalize        C
//      A(N)      R*8, I/O,  Input and Output Data                  C
//      X         I*4, In,   Constant Value                         C
//      N         I*4, In,   Number of Elements in A                C
//------------------------------------------------------------------C
//  Written by Y. Ushiro ( Waseda Uni. and Hitachi )  2003.01.12    C
//==================================================================C
{  int  i;
   for (i=0; i<N; i++) 
    { A[i] = A[i]*X ;  }            // A=A*X 
   NORM(A, N) ;  }                  // Fixed Normalize
//==================================================================C
void FNORM(double *A, int N)
//------------------------------------------------------------------C
//   A = Floating Normalize(A)                                      C
//     A(0:N)   R*8,Both, Input and Output Data                     C
//     N        I*4, In,  Number of Elements                        C
//------------------------------------------------------------------C
//  Written by Y. Ushiro ( Waseda Uni. and Hitachi )  2003.01.12    C
//------------------------------------------------------------------C
{  int    i, K ;
   double EXP ;
//   Save Exponent
   EXP = A[0] ;
//   Find Number of Zero (k) from Top
   for (i=1; i<=N; i++) 
	  { if(A[i] != 0.0) break ; }
   K = i - 1 ;
//   Shift K Word Left
   for (i=1; i<=N-K; i++) 
	  { A[i] = A[i+K] ; }
   SETV(0, &A[N-K+1], K) ;
//   Set Exponent
   A[0] = EXP - K ;   }
//==================================================================C
void NORM(double *A, int N)
//------------------------------------------------------------------C
//   A = Fixed Normalize(A) by Base-Value                           C
//     A(N)     R*8,Both, Input and Output Data                     C
//     N        I*4, In,  Number of Elements                        C
//------------------------------------------------------------------C
//  Written by Y. Ushiro ( Waseda Uni. and Hitachi )  2003.01.12    C
//------------------------------------------------------------------C
{ int i;
  double AI, B2;
//  Normalize by Base-Value (BASE)
  for (i=N-1; i>=1; i=i-1) {
     AI     = (int)( A[i]/BASE ) ;
	 A[i]   = A[i] - AI*BASE ;
     A[i-1] = A[i-1] + AI ; }
//  Check NORID (>=0; Positive Normalize, < 0 : Nearest Normalize )
  if(NORID >= 0) {
//   Set to A(i) >= 0 ( Positive Normalize )
     for (i=N-1; i>=1; i=i-1) {
        if(A[i] < 0.0) {
           A[i-1] = A[i-1] - 1.0 ;
           A[i]   = A[i] + BASE; } }
         }
  else {  
//   Set to |A(i)| <= BASE/2 ( Nearest Normalize )
     B2 = BASE/2.0 ;
     for (i=N-1; i>=1; i=i-1) {
        if(A[i] > B2) {
           A[i-1] = A[i-1] + 1 ;
           A[i]   = A[i] - BASE ; }
        else {  if (A[i] < -B2) {
           A[i-1] = A[i-1] - 1 ;
           A[i]   = A[i] + BASE ;  }
	   }   }
  }  }
//==================================================================C
void FLOAT(double *A, int N)
//------------------------------------------------------------------C
//   Change A from fixed to Floating                                C
//     A(0:N)   R*8,Both, Input and Output Data                     C
//     N        I*4, In,  Number of Elements                        C
//------------------------------------------------------------------C
//  Written by Y. Ushiro ( Waseda Uni. and Hitachi )  2003.01.12    C
//------------------------------------------------------------------C
//  Set Exponent and Floating Normalize
{  A[0] = N ;
   FNORM(A, N) ;  }
//==================================================================C
void COPY(double *A, double *C, int N)
//------------------------------------------------------------------C
//   Copy from A to C ( C = A )                                     C
//       A(N)     R*8, In,   Input A                                C
//       C(N)     R*8, Out,  Output C                               C
//       N        I*4, In,   Number of Elements of A,C              C
//------------------------------------------------------------------C
//  Written by Y. Ushiro ( Waseda Uni. and Hitachi )  2003.01.12    C
//==================================================================C
{  int i ;
//   Copy from A to C
   for (i=0; i<N; i++)
     { C[i] = A[i] ; }
   }
//==================================================================C
void SETV(int IV, double *C, int N)
//------------------------------------------------------------------C
//   Set Value V to C ( C(i) = IV )                                 C
//       IV       I*4, In,   Constant IV                            C
//       C(N)     R*8, Out,  Output C                               C
//       N        I*4, In,   Number of Elements of C                C
//------------------------------------------------------------------C
//  Written by Y. Ushiro ( Waseda Uni. and Hitachi )  2003.01.03    C
//==================================================================C
{  int i ; 
//   Set IV to C
   for (i=0; i<N; i++)
     { C[i] = IV ; }
   }
//==================================================================C
void FVSET(int IV, double *C, int N)
//------------------------------------------------------------------C
//   Set Floating Value to C  ( C = Floating IV )                   C
//       IV       I*4, In,   Constant IV                            C
//       C(0:N)   R*8, Out,  Output C                               C
//       N        I*4, In,   Number of Elements of C                C
//------------------------------------------------------------------C
//  Written by Y. Ushiro ( Waseda Uni. and Hitachi )  2003.01.12    C
//==================================================================C
//   Set Exponent
  {  C[0] = 1 ;
//   Set IV Value
     C[1] = IV ;
     SETV(0, &C[2], N-1) ; }
//==================================================================C
// N=N1**N2 Function
int power(int N1, int N2)
{ int k, N;
  N = 1;
  for (k=1; k<=N2; k++)
   { N = N*N1 ; }
  return (N) ; }
// SIGN Function
double SIGN(double a, double b)
{ double c ;
  if (b >= 0.0 ) { c = a ; }
  else { c= -a; }
  return (c) ;  }
// fmax Function
double fmax(double A1, double A2)
{ double C ;
  if (A1 >= A2) { C = A1 ; }
  else { C = A2; }
  return (C) ; }
// fmax3 Function
double fmax3(double A1, double A2, double A3)
{ double C ;
  C = fmax(A1, A2) ;
  C = fmax(C,  A3) ;
  return (C) ; }
//  MIN Function
int MIN(int N1, int N2)
{ int N;
  if(N1 <= N2) { N = N1 ; }
  else { N = N2 ; }
  return (N); }
// NB=log2(N) with int
int LOG2F(int N)
 {  int NB;  double DN ;
	 DN = N + 0.5 ;
     NB = log(DN)/log(2.0) ;
     return (NB) ; }