# Untitled

By: a guest on Jun 13th, 2012  |  syntax: None  |  size: 2.07 KB  |  hits: 26  |  expires: Never
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. }