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