Drobnicky.

diff --git a/Makefile b/Makefile
index c806aa287969e6c7dd6ed10106db84acacadd71a..781c602afb2a7f849d72d0717cb4da082d64815f 100644 (file)
--- a/Makefile
+++ b/Makefile
@@ -11,6 +11,9 @@ xclipcat: xclipcat.c
 xclipsend: xclipsend.c
 prefork: prefork.c
 
+fft: fft.c
+fft: LDFLAGS+=-lm
+
 batt: CFLAGS+=$(shell xosd-config --cflags)
 batt: LDFLAGS+=$(shell xosd-config --libs)
 

diff --git a/fft.c b/fft.c
new file mode 100644 (file)
index 0000000..973c8a8
--- /dev/null
+++ b/fft.c
@@ -0,0 +1,53 @@
+#define _GNU_SOURCE
+#include <stdio.h>
+#include <math.h>
+#include <complex.h>
+
+static int rev(int i, int n)
+{
+  int j = 0;
+  n--;
+  while (n)
+    {
+      j = (j << 1) | (i & 1);
+      i >>= 1;
+      n >>= 1;
+    }
+  return j;
+}
+
+static complex omega(int k, int n)
+{
+  return cos(2*M_PI*k/n) + I*sin(2*M_PI*k/n);
+}
+
+static void fft(double *x, double *y, int n)
+{
+  complex omegas[n];
+  for (int i=0; i<n; i++)
+    omegas[i] = omega(i, n);
+
+  for (int i=0; i<n; i++)
+    y[rev(i,n)] = x[i];
+
+  for (int i=1; i<n; i*=2)
+    for (int j=0; j<n; j+=2*i)
+      for (int k=0; k<i; k++)
+       {
+         complex a = y[j+k];
+         complex b = y[j+k+i];
+         complex o = omegas[(n/(2*i))*k % n];
+         y[j+k] = a + b*o;
+         y[j+k+i] = a - b*o;
+       }
+}
+
+int main(void)
+{
+  double x[8] = { 0, 1, 0, 1, 0, 1, 0, 1 };
+  double y[8];
+  fft(x, y, 8);
+  for (int i=0; i<8; i++)
+    printf("%f + %fi\n", creal(y[i]), cimag(y[i]));
+  return 0;
+}

diff --git a/pipprox.c b/pipprox.c
new file mode 100644 (file)
index 0000000..cf6d266
--- /dev/null
+++ b/pipprox.c
@@ -0,0 +1,21 @@
+#include <stdio.h>
+#include <math.h>
+
+int main(void)
+{
+  double delta = 1e30;
+
+  for (unsigned y=1; y < ~0U/M_PI; y++)
+    {
+      unsigned x = y*M_PI + 0.5;
+      double a = (double)x / y;
+      double d = fabs(a-M_PI);
+      if (d < delta)
+       {
+         delta = d;
+         printf("%u/%u = %.20f (%.20f)\n", x, y, a, d);
+       }
+    }
+
+  return 0;
+}