Version C des tests Cminor

git-svn-id: https://yquem.inria.fr/compcert/svn/compcert/trunk@40 fca1b0fc-160b-0410-b1d3-a4f43f01ea2e

1 files changed, 158 insertions, 0 deletions

diff --git a/test/c/fft.c b/test/c/fft.c
new file mode 100644
index 0000000..6f24af5
--- /dev/null
+++ b/test/c/fft.c
@@ -0,0 +1,158 @@
+/*
+ C Program: Test_Fast_Fourier_Transform_in_Double_Precision (TFFTDP.c)
+ by: Dave Edelblute, edelblut@cod.nosc.mil, 05 Jan 1993
+*/
+
+#include <math.h>
+
+#ifndef PI
+#define PI  3.14159265358979323846
+#endif
+
+extern double cos_static(double), sin_static(double);
+
+#define cos cos_static
+#define sin sin_static
+
+/********************************************************/
+/*     A Duhamel-Hollman split-radix dif fft            */
+/*     Ref: Electronics Letters, Jan. 5, 1984           */
+/*     Complex input and output data in arrays x and y  */
+/*     Length is n.                                     */
+/********************************************************/
+
+int dfft(double x[], double y[], int np)
+{
+double *px,*py;
+int i,j,k,m,n,i0,i1,i2,i3,is,id,n1,n2,n4;
+double  a,e,a3,cc1,ss1,cc3,ss3,r1,r2,s1,s2,s3,xt,tpi;
+
+  px = x - 1; 
+  py = y - 1;
+  i = 2; 
+  m = 1; 
+  
+  while (i < np) 
+  { 
+  i = i+i; 
+  m = m+1; 
+  }
+  
+  n = i; 
+  
+  if (n != np) 
+  {
+     for (i = np+1; i <= n; i++)  
+     { 
+     *(px + i) = 0.0; 
+     *(py + i) = 0.0; 
+     }
+     /*printf("nuse %d point fft",n); */
+  }
+
+  n2 = n+n;
+  tpi = 2.0 * PI;
+  for (k = 1;  k <= m-1; k++ ) 
+  {
+    n2 = n2 / 2; 
+    n4 = n2 / 4; 
+    e  = tpi / (double)n2; 
+    a  = 0.0;
+
+    for (j = 1; j<= n4 ; j++) 
+    {
+      a3 = 3.0 * a; 
+      cc1 = cos(a); 
+      ss1 = sin(a);
+      
+      cc3 = cos(a3); 
+      ss3 = sin(a3); 
+      a = e * (double)j; 
+      is = j; 
+      id = 2 * n2;
+	  
+	  while ( is < n ) 
+	  {
+	     for (i0 = is; i0 <= n-1; i0 = i0 + id) 
+	     {
+	     i1 = i0 + n4; 
+	     i2 = i1 + n4; 
+	     i3 = i2 + n4;
+	     r1 = *(px+i0) - *(px+i2);
+	     *(px+i0) = *(px+i0) + *(px+i2);
+	     r2 = *(px+i1) - *(px+i3);
+	     *(px+i1) = *(px+i1) + *(px+i3);
+	     s1 = *(py+i0) - *(py+i2);
+	     *(py+i0) = *(py+i0) + *(py+i2);
+	     s2 = *(py+i1) - *(py+i3);
+	     *(py+i1) = *(py+i1) + *(py+i3);
+	     s3 = r1 - s2; r1 = r1 + s2; 
+	     s2 = r2 - s1; r2 = r2 + s1;
+	     *(px+i2) = r1*cc1 - s2*ss1; 
+	     *(py+i2) = -s2*cc1 - r1*ss1;
+	     *(px+i3) = s3*cc3 + r2*ss3;
+	     *(py+i3) = r2*cc3 - s3*ss3;
+	     }
+	  is = 2 * id - n2 + j; 
+	  id = 4 * id;
+	  }
+    }
+  }
+
+/************************************/
+/*  Last stage, length=2 butterfly  */
+/************************************/
+  is = 1; 
+  id = 4;
+  
+  while ( is < n) 
+  {
+     for (i0 = is; i0 <= n; i0 = i0 + id) 
+     {
+      i1 = i0 + 1; 
+      r1 = *(px+i0);
+      *(px+i0) = r1 + *(px+i1);
+      *(px+i1) = r1 - *(px+i1);
+      r1 = *(py+i0);
+      *(py+i0) = r1 + *(py+i1);
+      *(py+i1) = r1 - *(py+i1);
+     }
+  is = 2*id - 1; 
+  id = 4 * id; 
+  }
+
+/*************************/
+/*  Bit reverse counter  */
+/*************************/
+  j = 1; 
+  n1 = n - 1;
+  
+  for (i = 1; i <= n1; i++) 
+  {
+    if (i < j) 
+    {
+    xt = *(px+j);
+    *(px+j) = *(px+i); 
+    *(px+i) = xt;
+    xt = *(py+j); 
+    *(py+j) = *(py+i);
+    *(py+i) = xt;
+    }
+    
+    k = n / 2; 
+    while (k < j) 
+    { 
+    j = j - k; 
+    k = k / 2; 
+    }
+    j = j + k;
+  }
+
+/*
+  for (i = 1; i<=16; i++) printf("%d  %g   %gn",i,*(px+i),(py+i));
+*/
+  
+  return(n);
+}
+
+