Don't like ads? PRO users don't see any ads ;-)
Guest

Untitled

By: a guest on Jun 13th, 2012  |  syntax: None  |  size: 2.07 KB  |  hits: 26  |  expires: Never
download  |  raw  |  embed  |  report abuse  |  print
Text below is selected. Please press Ctrl+C to copy to your clipboard. (⌘+C on Mac)
  1. #include <stdio.h>
  2. #include <math.h>
  3.  
  4. unsigned int evaluations; /* This EVIL global variable stores the number of times curve_length has been called for the current curve. */
  5.  
  6. typedef struct {double x, y;} point;
  7.  
  8. typedef point (*curve_func)(double t);
  9.  
  10. point line(double t) {
  11.         point p;
  12.         p.x = t;
  13.         p.y = t * 2;
  14.         return p;
  15. }
  16.  
  17. point parabola(double t) {
  18.         point p;
  19.         p.x = t;
  20.         p.y = t * t;
  21.         return p;
  22. }
  23.  
  24. point sqrt_sin(double t) {
  25.         point p;
  26.         p.x = t;
  27.         p.y = sqrt(sin(t));
  28.         return p;
  29. }
  30.  
  31. double distance(point a, point b) {
  32.         double dx = (a.x - b.x);
  33.         double dy = (a.y - b.y);
  34.         return sqrt(dx * dx + dy * dy);
  35. }
  36.  
  37. double curve_length(double a, double b, curve_func curve, double tolerance, unsigned int depth) {
  38.         if (depth == 1) {
  39.                 evaluations = 1;
  40.         } else {
  41.                 evaluations++;
  42.         }
  43.         double avg_ab = (a + b) / 2;
  44.         point curve_at_a = (*curve)(a);
  45.         point curve_at_b = (*curve)(b);
  46.         double coarse = distance(curve_at_a, curve_at_b);
  47.         point curve_at_avg_ab = (*curve)(avg_ab);
  48.         double fine = (
  49.                 distance(
  50.                         curve_at_a,
  51.                         curve_at_avg_ab
  52.                 ) +
  53.                 distance(
  54.                         curve_at_avg_ab,
  55.                         curve_at_b
  56.                 )
  57.         );
  58.         if (depth == 30) {
  59.                 return fine;
  60.         } else {
  61.                 if ((depth >= 3) && (fabs(fine - coarse) <= tolerance)) {
  62.                         return fine;
  63.                 } else {
  64.                         double new_tolerance = tolerance / 2;
  65.                         unsigned int new_depth = depth + 1;
  66.                         return (
  67.                                 curve_length(a,      avg_ab, curve, new_tolerance, new_depth) +
  68.                                 curve_length(avg_ab, b,      curve, new_tolerance, new_depth)
  69.                         );
  70.                 }
  71.         }
  72. }
  73.  
  74. void print_curve_length(double a, double b, curve_func curve, double tolerance) {
  75.         double length = curve_length(a, b, curve, tolerance, 1);
  76.         printf("  tol %4f : %f, %d evaluations\n", tolerance, length, evaluations);
  77. }
  78.  
  79. int main() {
  80.         printf("Line\n");
  81.         print_curve_length(0, 2, &line, .01);
  82.         print_curve_length(0, 2, &line, .0001);
  83.        
  84.         printf("Parabola\n");
  85.         print_curve_length(0, 1, &parabola, .01);
  86.         print_curve_length(0, 1, &parabola, .0001);
  87.        
  88.         printf("Sqrt sin\n");
  89.         print_curve_length(0, 1, &sqrt_sin, .01);
  90.         print_curve_length(0, 1, &sqrt_sin, .0001);
  91.        
  92.         return 0;
  93. }