forked from Github_Repos/cvw
115 lines
2.6 KiB
C
Executable File
115 lines
2.6 KiB
C
Executable File
#include "disp.h"
|
|
|
|
// QSLC is for division by recuerrence for
|
|
// r=2 using a CPA - See 5.109 EL
|
|
int qst (double D, double prem) {
|
|
|
|
int q;
|
|
|
|
// For Debugging
|
|
printf("rw --> %lg\n", prem);
|
|
|
|
if (prem >= 0.5) {
|
|
q = 1;
|
|
} else if (prem >= -0.5) {
|
|
q = 0;
|
|
} else {
|
|
q = -1;
|
|
}
|
|
return q;
|
|
|
|
}
|
|
|
|
/*
|
|
This routine performs a radix-2 SRT division
|
|
algorithm. The user inputs the numerator, the denominator,
|
|
and the number of iterations. It assumes that 0.5 <= D < 1.
|
|
|
|
*/
|
|
|
|
int main(int argc, char* argv[]) {
|
|
|
|
double P, N, D, Q, RQ, RD, RREM, scale;
|
|
int q;
|
|
int num_iter, i;
|
|
int prec;
|
|
int radix = 2;
|
|
|
|
if (argc < 5) {
|
|
fprintf(stderr,
|
|
"Usage: %s numerator denominator num_iterations prec\n",
|
|
argv[0]);
|
|
exit(1);
|
|
}
|
|
sscanf(argv[1],"%lg", &N);
|
|
sscanf(argv[2],"%lg", &D);
|
|
sscanf(argv[3],"%d", &num_iter);
|
|
sscanf(argv[4],"%d", &prec);
|
|
// Round to precision
|
|
N = rne(N, prec);
|
|
D = rne(D, prec);
|
|
printf("N = ");
|
|
disp_bin(N, 3, prec, stdout);
|
|
printf("\n");
|
|
printf("D = ");
|
|
disp_bin(D, 3, prec, stdout);
|
|
printf("\n");
|
|
|
|
Q = 0;
|
|
P = N * pow(2.0, -log2(radix));
|
|
printf("N = %lg, D = %lg, N/D = %lg, num_iter = %d \n\n",
|
|
N, D, N/D, num_iter);
|
|
for (scale = 1, i = 0; i < num_iter; i++) {
|
|
scale = scale * pow(2.0, -log2(radix));
|
|
q = qst(flr(2*D, 1), 2*P);
|
|
printf("2*W[n] = ");
|
|
disp_bin(radix*P, 3, prec, stdout);
|
|
printf("\n");
|
|
printf("q*D = ");
|
|
disp_bin(q*D, 3, prec, stdout);
|
|
printf("\n");
|
|
printf("W[n+1] = ");
|
|
disp_bin(P ,3, prec, stdout);
|
|
printf("\n");
|
|
// Recurrence
|
|
P = radix * P - q * D;
|
|
Q = Q + q*scale;
|
|
printf("i = %d, q = %d, Q = %1.18lf, W = %1.18lf\n", i, q, Q, P);
|
|
printf("i = %d, q = %d", i, q);
|
|
printf(", Q = ");
|
|
disp_bin(Q, 3, prec, stdout);
|
|
printf(", W = ");
|
|
disp_bin(P, 3, prec, stdout);
|
|
printf("\n\n");
|
|
}
|
|
if (P < 0) {
|
|
Q = Q - scale;
|
|
P = P + D;
|
|
printf("\nCorrecting Negative Remainder\n");
|
|
printf("Q = %1.18lf, W = %1.18lf\n", Q, P);
|
|
printf("Q = ");
|
|
disp_bin(Q, 3, prec, stdout);
|
|
printf(", W = ");
|
|
disp_bin(P, 3, prec, stdout);
|
|
printf("\n");
|
|
}
|
|
|
|
// Output Results
|
|
RQ = N/D;
|
|
// Since q_{computed} = q / radix, multiply by radix
|
|
RD = Q * radix;
|
|
printf("true = %1.18lf, computed = %1.18lf, \n", RQ, RD);
|
|
printf("true = ");
|
|
disp_bin(RQ, 3, prec, stdout);
|
|
printf(", computed = ");
|
|
disp_bin(RD, 3, prec, stdout);
|
|
printf("\n\n");
|
|
printf("REM = %1.18lf \n", P);
|
|
printf("REM = ");
|
|
disp_bin(P, 3, prec, stdout);
|
|
printf("\n\n");
|
|
|
|
return 0;
|
|
|
|
}
|