summaryrefslogtreecommitdiff
path: root/test/c
diff options
context:
space:
mode:
Diffstat (limited to 'test/c')
-rw-r--r--test/c/Makefile4
-rw-r--r--test/c/Results/fftsp1
-rw-r--r--test/c/fftsp.c196
3 files changed, 199 insertions, 2 deletions
diff --git a/test/c/Makefile b/test/c/Makefile
index d14cb23..1444073 100644
--- a/test/c/Makefile
+++ b/test/c/Makefile
@@ -10,8 +10,8 @@ LIBS=$(LIBMATH)
TIME=xtime -o /dev/null -mintime 2.0 # Xavier's hack
#TIME=time >/dev/null # Otherwise
-PROGS=fib integr qsort fft fftw sha1 sha3 aes almabench lists \
- binarytrees fannkuch knucleotide mandelbrot nbody \
+PROGS=fib integr qsort fft fftsp fftw sha1 sha3 aes almabench \
+ lists binarytrees fannkuch knucleotide mandelbrot nbody \
nsieve nsievebits spectral vmach \
bisect chomp perlin siphash24
diff --git a/test/c/Results/fftsp b/test/c/Results/fftsp
new file mode 100644
index 0000000..cbeb099
--- /dev/null
+++ b/test/c/Results/fftsp
@@ -0,0 +1 @@
+4096 points, result OK
diff --git a/test/c/fftsp.c b/test/c/fftsp.c
new file mode 100644
index 0000000..42ae905
--- /dev/null
+++ b/test/c/fftsp.c
@@ -0,0 +1,196 @@
+/*
+ C Program: Test_Fast_Fourier_Transform_in_Double_Precision (TFFTDP.c)
+ by: Dave Edelblute, edelblut@cod.nosc.mil, 05 Jan 1993
+*/
+
+/* Like fft.c, but using single-precision floats */
+
+#include <math.h>
+#include <stdlib.h>
+#include <stdio.h>
+
+#ifndef PI
+#define PI 3.14159265358979323846
+#endif
+
+/********************************************************/
+/* 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(float x[], float y[], int np)
+{
+ float *px,*py;
+ int i,j,k,m,n,i0,i1,i2,i3,is,id,n1,n2,n4;
+ float 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 = cosf(a);
+ ss1 = sinf(a);
+
+ cc3 = cosf(a3);
+ ss3 = sinf(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);
+}
+
+/* Test harness */
+
+float * xr, * xi;
+
+int main(int argc, char ** argv)
+{
+ int n, np, npm, n2, i, j;
+ float enp, t, y, z, zr, zi, zm, a;
+ float * pxr, * pxi;
+
+ if (argc >= 2) n = atoi(argv[1]); else n = 12;
+ np = 1 << n;
+ enp = np;
+ npm = np / 2 - 1;
+ t = PI / enp;
+ xr = calloc(np, sizeof(float));
+ xi = calloc(np, sizeof(float));
+ pxr = xr;
+ pxi = xi;
+ *pxr = (enp - 1.0) * 0.5;
+ *pxi = 0.0;
+ n2 = np / 2;
+ *(pxr+n2) = -0.5;
+ *(pxi+n2) = 0.0;
+ for (i = 1; i <= npm; i++) {
+ j = np - i;
+ *(pxr+i) = -0.5;
+ *(pxr+j) = -0.5;
+ z = t * (double)i;
+ y = -0.5*(cosf(z)/sinf(z));
+ *(pxi+i) = y;
+ *(pxi+j) = -y;
+ }
+ dfft(xr,xi,np);
+ zr = 0.0;
+ zi = 0.0;
+ npm = np-1;
+ for (i = 0; i <= npm; i++ ) {
+ a = fabsf(pxr[i] - i);
+ if (zr < a) zr = a;
+ a = fabsf(pxi[i]);
+ if (zi < a) zi = a;
+ }
+ zm = zr;
+ if (zr < zi) zm = zi;
+ if (zm < 1e-3)
+ printf("%d points, result OK\n", np);
+ else
+ printf("%d points, WRONG result %e\n", np, zm);
+ return 0;
+}