Add r=4 C code

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

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, ".");
+    }
+  }
+}
+

pipelined/srt/stine/disp.h

@@ -0,0 +1,18 @@
+#include <stdlib.h>
+#include <math.h>
+#include <stdio.h>
+
+#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 

pipelined/srt/stine/srt4div.c

@@ -0,0 +1,325 @@
+#include "disp.h"
+#include <math.h>
+
+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;
+
+}