diff --git a/pipelined/srt/stine/Makefile b/pipelined/srt/stine/Makefile new file mode 100644 index 000000000..a019b3880 --- /dev/null +++ b/pipelined/srt/stine/Makefile @@ -0,0 +1,18 @@ + +CC = gcc +CFLAGS = -lm +LIBS = +OBJS = disp.o srt4div.o + +srt4div: $(OBJS) + $(CC) -g -O3 -o srt4div $(OBJS) $(CFLAGS) + +disp.o: disp.h disp.c + $(CC) -g -c -o disp.o disp.c $(CFLAGS) + +srt4div.o: srt4div.c + $(CC) -g -c -o srt4div.o srt4div.c $(CFLAGS) + +clean: + rm -f *.o *~ + rm -f core diff --git a/pipelined/srt/stine/disp.c b/pipelined/srt/stine/disp.c new file mode 100755 index 000000000..66652cf01 --- /dev/null +++ b/pipelined/srt/stine/disp.c @@ -0,0 +1,60 @@ +#include "disp.h" + +double rnd_zero(double x, double bits) { + if (x < 0) + return ceiling(x, bits); + else + return flr(x, bits); +} + +double rne(double x, double precision) { + double scale, x_round; + scale = pow(2.0, precision); + x_round = rint(x * scale) / scale; + return x_round; +} + +double flr(double x, double precision) { + double scale, x_round; + scale = pow(2.0, precision); + x_round = floor(x * scale) / scale; + return x_round; +} + +double ceiling(double x, double precision) { + double scale, x_round; + scale = pow(2.0, precision); + x_round = ceil(x * scale) / scale; + return x_round; +} + +void disp_bin(double x, int bits_to_left, int bits_to_right, FILE *out_file) { + + double diff; + int i; + if (fabs(x) < pow(2.0, -bits_to_right)) { + for (i = -bits_to_left + 1; i <= bits_to_right; i++) { + fprintf(out_file,"0"); + } + return; + } + if (x < 0.0) { + // fprintf(out_file, "-"); + // x = - x; + x = pow(2.0, ((double) bits_to_left)) + x; + } + for (i = -bits_to_left + 1; i <= bits_to_right; i++) { + diff = pow(2.0, -i); + if (x < diff) { + fprintf(out_file, "0"); + } + else { + fprintf(out_file, "1"); + x -= diff; + } + if (i == 0) { + fprintf(out_file, "."); + } + } +} + diff --git a/pipelined/srt/stine/disp.h b/pipelined/srt/stine/disp.h new file mode 100755 index 000000000..bf2c78872 --- /dev/null +++ b/pipelined/srt/stine/disp.h @@ -0,0 +1,18 @@ +#include +#include +#include + +#ifndef DISP +#define DISP + +double rnd_zero(double x, double bits); + +double rne(double x, double precision); + +double flr(double x, double precision); + +double ceiling(double x, double precision); + +void disp_bin(double x, int bits_to_left, int bits_to_right, FILE *out_file); + +#endif diff --git a/pipelined/srt/stine/srt4div b/pipelined/srt/stine/srt4div new file mode 100755 index 000000000..ff76638fa Binary files /dev/null and b/pipelined/srt/stine/srt4div differ diff --git a/pipelined/srt/stine/srt4div.c b/pipelined/srt/stine/srt4div.c new file mode 100755 index 000000000..cbf0d8224 --- /dev/null +++ b/pipelined/srt/stine/srt4div.c @@ -0,0 +1,325 @@ +#include "disp.h" +#include + +int qslc (double prem, double d) { + + int q; + + printf("d --> %lg\n", d); + printf("rw --> %lg\n", prem); + if ((d>=0.0)&&(d<1.0)) { + if (prem>=1.0) + q = 2; + else if (prem>=0.25) + q = 1; + else if (prem>=-0.25) + q = 0; + else if (prem >= -1) + q = -1; + else + q = -2; + return q; + } + + if ((d>=1.0)&&(d<2.0)) { + if (prem>=2.0) + q = 2; + else if (prem>=0.66667) + q = 1; + else if (prem>=-0.6667) + q = 0; + else if (prem >= -2) + q = -1; + else + q = -2; + return q; + } + + if ((d>=2.0)&&(d<3.0)) { + if (prem>=4.0) + q = 2; + else if (prem>=1.25) + q = 1; + else if (prem>=-1.25) + q = 0; + else if (prem >= -4) + q = -1; + else + q = -2; + return q; + } + + if ((d>=3.0)&&(d<4.0)) { + if (prem>=5.0) + q = 2; + else if (prem>=2.0) + q = 1; + else if (prem>=-2.0) + q = 0; + else if (prem >= -5) + q = -1; + else + q = -2; + return q; + } + + if ((d>=4.0)&&(d<5.0)) { + if (prem>=6.66667) + q = 2; + else if (prem>=2.0) + q = 1; + else if (prem>=-2.0) + q = 0; + else if (prem >= -6.66667) + q = -1; + else + q = -2; + return q; + } + + if ((d>=5.0)&&(d<6.0)) { + if (prem>=8.0) + q = 2; + else if (prem>=2.0) + q = 1; + else if (prem>=-2.0) + q = 0; + else if (prem >= -8.0) + q = -1; + else + q = -2; + return q; + } + + if ((d>=6.0)&&(d<7.0)) { + if (prem>=10.0) + q = 2; + else if (prem>=4.0) + q = 1; + else if (prem>=-4.0) + q = 0; + else if (prem >= -10.0) + q = -1; + else + q = -2; + return q; + } + + if ((d>=7.0)&&(d<8.0)) { + if (prem>=11.0) + q = 2; + else if (prem>=4.0) + q = 1; + else if (prem>=-4.0) + q = 0; + else if (prem >= -11.0) + q = -1; + else + q = -2; + return q; + } + + if ((d>=8.0)&&(d<9.0)) { + if (prem>=12.0) + q = 2; + else if (prem>=4.0) + q = 1; + else if (prem>=-4.0) + q = 0; + else if (prem >= -12.0) + q = -1; + else + q = -2; + return q; + } + + if ((d>=9.0)&&(d<10.0)) { + if (prem>=15.0) + q = 2; + else if (prem>=4.0) + q = 1; + else if (prem>=-4.0) + q = 0; + else if (prem >= -15.0) + q = -1; + else + q = -2; + return q; + } + + if ((d>=10.0)&&(d<11.0)) { + if (prem>=15.0) + q = 2; + else if (prem>=4.0) + q = 1; + else if (prem>=-4.0) + q = 0; + else if (prem >= -15.0) + q = -1; + else + q = -2; + return q; + } + + if ((d>=11.0)&&(d<12.0)) { + if (prem>=16.0) + q = 2; + else if (prem>=4.0) + q = 1; + else if (prem>=-4.0) + q = 0; + else if (prem >= -16.0) + q = -1; + else + q = -2; + return q; + } + + if ((d>=12.0)&&(d<13.0)) { + if (prem>=20.0) + q = 2; + else if (prem>=8.0) + q = 1; + else if (prem>=-8.0) + q = 0; + else if (prem >= -20.0) + q = -1; + else + q = -2; + return q; + } + + if ((d>=13.0)&&(d<14.0)) { + if (prem>=20.0) + q = 2; + else if (prem>=8.0) + q = 1; + else if (prem>=-8.0) + q = 0; + else if (prem >= -20.0) + q = -1; + else + q = -2; + return q; + } + + if ((d>=14.0)&&(d<15.0)) { + if (prem>=20.0) + q = 2; + else if (prem>=8.0) + q = 1; + else if (prem>=-8.0) + q = 0; + else if (prem >= -20.0) + q = -1; + else + q = -2; + return q; + } + + if ((d>=15.0)&&(d<16.0)) { + if (prem>=24.0) + q = 2; + else if (prem>=8.0) + q = 1; + else if (prem>=-8.0) + q = 0; + else if (prem >= -24.0) + q = -1; + else + q = -2; + return q; + } + +} + +/* + This routine performs a radix-4 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; + 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*0.25; + 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++) { + // Shift by r + scale = scale*0.25; + q = qslc(flr((4*P)*16,3), D*16); + //q = -q; + printf("4*W[n] = "); + disp_bin(4*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 = 4*P - q*D; + // OTFC + 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"); + } + RQ = flr(N/D, (double) prec); + RD = Q*4; + 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; + +}