From cf3a29feb5887344b6633ead1b4b6d5657a15a4b Mon Sep 17 00:00:00 2001
From: Patrick Simianer <p@simianer.de>
Date: Sun, 15 Jun 2014 03:24:33 +0200
Subject: old stuff: algorithms

---
 algorithms/polynomial.c | 68 +++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 68 insertions(+)
 create mode 100644 algorithms/polynomial.c

(limited to 'algorithms/polynomial.c')

diff --git a/algorithms/polynomial.c b/algorithms/polynomial.c
new file mode 100644
index 0000000..ecd9027
--- /dev/null
+++ b/algorithms/polynomial.c
@@ -0,0 +1,68 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+
+#define EPSILON 1E-9
+
+
+double
+eval_poly(double x, double coeffs[], int order)
+{
+  double res;
+  int i;
+
+  res = coeffs[0];
+
+  for ( i = 1; i <= order; i++ ) {
+    res = res * x + coeffs[i];
+  }
+
+  return res;
+}
+
+
+double
+find_zero(double left, double right, double coeffs[], int order)
+{
+  double c, lp, rp, cp, eps=EPSILON;
+
+  if(eval_poly(left, coeffs, order) * eval_poly(right, coeffs, order) >= 0) {
+    exit(1);
+  }
+
+  while(1) {
+    c = (left+right)/2;
+
+    lp = eval_poly(left, coeffs, order);
+    rp = eval_poly(right, coeffs, order);
+    cp = eval_poly(c, coeffs, order);
+
+    if(fabs(cp) < eps) {
+      break;
+    }
+
+    if(cp * lp < 0) {
+      right = c;
+    }
+    else if(cp * rp < 0) {
+      left = c;
+    }
+    else {
+      exit(1);
+    }
+  }
+
+  return c;
+}
+
+int
+main(void)
+{
+  double coeffs[] = { 7, 5, 3, 2, 1, 1 };
+
+  printf("eval_poly: %f\n", eval_poly(3, coeffs, 5));
+  printf("%f\n", find_zero(-100, 100, coeffs, 5));
+
+  return 0;
+}
+
-- 
cgit v1.2.3