From 4b0dfce7e0526346df2daa268d787bb7e80f4036 Mon Sep 17 00:00:00 2001
From: Martin Mares <mj@ucw.cz>
Date: Sun, 26 Oct 2008 00:43:51 +0200
Subject: [PATCH] Drobnicky.

---
 Makefile  |  3 +++
 fft.c     | 53 +++++++++++++++++++++++++++++++++++++++++++++++++++++
 pipprox.c | 21 +++++++++++++++++++++
 3 files changed, 77 insertions(+)
 create mode 100644 fft.c
 create mode 100644 pipprox.c

diff --git a/Makefile b/Makefile
index c806aa2..781c602 100644
--- 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
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
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;
+}
-- 
2.39.2