Fixes analysis program

Author: Jason Boutte

lab1/analysis.c

@@ -1,136 +1,85 @@
 #include <math.h>
-#include <float.h>
 #include <stdio.h>
-#include <stdlib.h>
 
-#define SPREAD 4
+#ifdef DEBUG
+#define debugf(x, ...) fprintf(stdout, x, ##__VA_ARGS__);
+#else
+#define debugf(x, ...)
+#endif
+
 #define ABS(x)(((x)<0)?-(x):(x))
 #define MAX_ERROR 0.5*pow(10, -12)
 
-int search(int traps, int tmin, int tmax, double a, double b);
-double calc_true_error(double true_val, double approx_val);
+double calc_approx(double a, double b, int traps);
 double calc_actual(double a, double b);
-double calc_approx(double a, double b, int traps, double h);
 double func(double x);
 double func_int(double x);
 
-double current_min;
-
 int main(int argc, char **argv) {
-  int traps, mid, optimal_traps;
-  double lend, rend, approx, actual;
-
-  printf("Enter left, right, and trapezoids\n");
-
-  scanf("%lf%lf%d", &lend, &rend, &traps);
-
-  mid = (traps/2);
-
-  optimal_traps = search(mid, 0, traps, lend, rend);
-
-  printf("t = %d\n", optimal_traps);
-
-  approx = calc_approx(lend, rend, optimal_traps, (rend-lend)/optimal_traps); 
-
-  actual = calc_actual(lend, rend);
-
-  printf("absolute relative true error %.16lf\n", calc_true_error(actual, approx));
-
-	return 0;
-}
-
-int search(int traps, int tmin, int tmax, double a, double b) {
-  int x, lvalue, rvalue, lmid, rmid, min_index;
-  double h, approx, actual, true_error, min_value;
+	int stride, index;
+	double a, b, actual, approx, true_error;
 
-  if (traps <= tmin || traps > tmax) return traps;
+	printf("Enter a and b\n");
+	
+	scanf("%lf%lf", &a, &b);
 
-  min_value = DBL_MAX;
+	actual = calc_actual(a, b);
 
-  for (x = traps-SPREAD; x < traps+SPREAD+1; ++x) {
-    h = (b-a)/x;
+	stride = 4096;
 
-    approx = calc_approx(a, b, x, h); 
+	while (1) {
+		if (stride == 1) break;
 
-    actual = calc_actual(a, b);
+		printf("Index at %d\n", index);
+		printf("Stride at %d\n", stride);
 
-    true_error = calc_true_error(actual, approx);
+		while (1) {
+			index += stride;
 
-    if (true_error <= MAX_ERROR) {
-      printf("traps %d\n", x);
-      printf("h %.16lf\n", h);
-      printf("approx %.16lf\n", approx);
-      printf("actual %.16lf\n", actual);
-      printf("true error %.13lf\n", true_error);
-      printf("true error %.13lf\n", MAX_ERROR);
+			approx = calc_approx(a, b, index);
 
-      return x;
-    }
+			true_error = ABS(actual-approx)*100/actual;
 
-    if (true_error < min_value) {
-      min_index = x;
+			printf("traps %d\n", index);
+			printf("actual %.10lf\n", actual);
+			printf("approx %.10lf\n", approx);
+			printf("true error %.16e\n", true_error);
 
-      min_value = true_error;
-    }
+			if (true_error <= MAX_ERROR) {
+				printf("Found true error less than equal %.16e\n", MAX_ERROR);
+		
+				break;
+			}
+		}
 
-    printf("traps %d\n", x);
-    printf("h %.16lf\n", h);
-    printf("approx %.16lf\n", approx);
-    printf("actual %.16lf\n", actual);
-    printf("true error %.13lf\n", true_error);
-    printf("true error %.13lf\n", MAX_ERROR);
-  }
+		index -= stride;
 
-  printf("Done looking at values found min index %d\n", min_index);
+		stride /= 2;
+	}
 
-  if (min_value < current_min) {
-    current_min = min_value;
-  } else {
-    min_index = traps+1;
-  }
-  
-  if (min_index <= traps) {
-    lvalue = floor((traps-tmin)/2);
-    lmid = tmin+lvalue;
+	return 0;
+}
 
-    printf("We're going left\n");
-    return search(lmid, tmin, traps, a, b);
-  } else {
-    rvalue = floor((tmax-traps)/2);
-    rmid = traps+rvalue;
+double calc_approx(double a, double b, int n) {
+	int i;
+	double h = (b-a)/n;
+	double approx = (func(a)+func(b))/2;
 
-    printf("We're going right\n");
-    return search(rmid, traps, tmax, a, b);
-  }
-}
+	for (i = 1; i <= n-1; ++i) {
+		approx += func(a+i*h);		
+	}	
 
-double calc_true_error(double true_val, double approx_val) {
-  return ABS(true_val-approx_val)/true_val;
+	return h * approx;
 }
 
 double calc_actual(double a, double b) {
-  return func_int(b)-func_int(a);
-}
-
-double calc_approx(double a, double b, int traps, double h) {
-  int x;
-  double approx = 0.0;
-
-  for (x = 1; x < traps-1; ++x) {
-    approx += func(a+x*h);
-  }
-
-  approx *= 2;
-  approx += func(a)+func(b);
-  approx *= h/2;
-
-  return approx;
+	return func_int(b)-func_int(a);
 }
 
 double func(double x) {
-  return cos(x/3)-2*cos(x/5)+5*sin(x/4)+8;
+	return 8+5*sin(x/4)-2*cos(x/5)+cos(x/3);
 }
 
 double func_int(double x) {
-  return 8*x-10*sin(x/5)+3*sin(x/3)-20*cos(x/4);
+	return 8*x-10*sin(x/5)+3*sin(x/3)-20*cos(x/4);
 }