多数桁乗除算Cプログラム
2003/01/10 日立製作所 & 早稲田大学 後 保範 ( Ushiro Yasunori )
--------------------------------------------------------------
#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] ; }
}