Check out the Pastebin Gadgets Shop. We have thousands of fun, geeky & affordable gadgets on sale :-)Want more features on Pastebin? Sign Up, it's FREE!
tweet

# Untitled

By: a guest on Feb 26th, 2012  |  syntax: C  |  size: 4.57 KB  |  views: 24  |  expires: Never
Text below is selected. Please press Ctrl+C to copy to your clipboard. (⌘+C on Mac)
1. #include <stdio.h>
2. #include <stdlib.h>
3. #include <math.h>
4.
5. const double EPS1 = 1e-5;
6. const double EPS2 = 1e-5;
7.
8. typedef struct Slice {
9.     int n;
10.     double *keys, *values;
11. } Slice;
12.
13. double f(int k, double x) {
14.     switch(k) {
15.         case 0: return exp(x) + 2;
16.         case 1: return -1/x;
17.         case 2: return -2*(x + 1)/3;
18.         default: return 0;
19.     }
20. }
21.
22. double fp(int k, double x) {
23.      switch(k) {
24.         case 0: return exp(x);
25.         case 1: return 1/(x * x);
26.         case 2: return -2*x/3;
27.         default: return 0;
28.     }
29. }
30.
31. double rootChord(int f1, int f2, double l, double r, double lv, double rv) {
32.     double m = (l * rv - r * lv) / (rv - lv);
33.     if (fabs(m - l) < EPS1 || fabs(r - m) < EPS1) {
34.         return m;
35.     }
36.     double mv = f(f1, m) - f(f2, m);
37.     if ((lv < 0) != (mv < 0)) {
38.         return rootChord(f1, f2, l, m, lv, mv);
39.     } else {
40.         return rootChord(f1, f2, m, r, mv, rv);
41.     }
42. }
43.
44. double rootTangent(int f1, int f2, double l, double r, double lv, double rv, double lp, double rp) {
45.     //fprintf(stderr, "%.4lf, %.4lf\n", l, r);
46.     //system("pause");
47.     double m1 = l - lv / lp;
48.     double m2 = r - rv / rp;
49.     double m = (l < m1 && m1 < r) ? m1 : m2;
50.     if (fabs(m - l) < EPS1 || fabs(r - m) < EPS1) {
51.         return m;
52.     }
53.     double mv = f(f1, m) - f(f2, m);
54.     double mp = fp(f1, m) - fp(f2, m);
55.     if ((lv < 0) != (mv < 0)) {
56.         return rootTangent(f1, f2, l, m, lv, mv, lp, mp);
57.     } else {
58.         return rootTangent(f1, f2, m, r, mv, rv, mp, rp);
59.     }
60. }
61.
62. double root(int f1, int f2, double l, double r) {
63.     //fprintf(stderr, "Root: %.4lf, %.4lf\n", l, r);
64.     double lv = f(f1, l) - f(f2, l);
65.     double rv = f(f1, r) - f(f2, r);
66.     double lp = fp(f1, l) - fp(f2, l);
67.     double rp = fp(f1, r) - fp(f2, r);
68.     //fprintf(stderr, "C\n");
69.     double root1 = rootChord(f1, f2, l, r, lv, rv);
70.     //fprintf(stderr, "T\n");
71.     double root2 = rootTangent(f1, f2, l, r, lv, rv, lp, rp);
72.     //fprintf(stderr, "OK\n");
73.     return root1;
74. }
75.
76. void expandArray(int fn, Slice *s) {
77.     s->n *= 2;
78.     s->keys = realloc(s->keys, (s->n + 1) * sizeof(double));
79.     s->values = realloc(s->values, (s->n + 1) * sizeof(double));
80.     for (int i = s->n; i >= 2; i -= 2) {
81.         s->keys[i] = s->keys[i / 2];
82.         s->values[i] = s->values[i / 2];
83.     }
84.     for (int i = 1; i < s->n; i += 2) {
85.         s->keys[i] = (s->keys[i - 1] + s->keys[i + 1]) / 2.0;
86.         s->values[i] = f(fn, s->keys[i]);
87.     }
88. }
89.
90. Slice newSlice(int fn, double l, double r) {
91.     Slice res;
92.     res.n = 1;
93.     res.keys = malloc(2 * sizeof(double));
94.     res.values = malloc(2 * sizeof(double));
95.     res.keys[0] = l;
96.     res.values[0] = f(fn, l);
97.     res.keys[1] = r;
98.     res.values[1] = f(fn, r);
99.     expandArray(fn, &res);
100.     return res;
101. }
102.
103. double integrateArray(Slice s) {
104.     double res = s.values[0] + s.values[s.n];
105.     for (int i = 1; i < s.n; ++i) {
106.         if (i % 2) {
107.             res += 2 * s.values[i];
108.         } else {
109.             res += 4 * s.values[i];
110.         }
111.     }
112.     res *= fabs(s.keys[0] - s.keys[s.n]) / (s.n * 3);
113.     return res;
114. }
115.
116. double integrate(int fn, double l, double r) {
117.     Slice sl = newSlice(fn, l, r);
118.     double prev, cur = integrateArray(sl);
119.     do {
120.         prev = cur;
121.         expandArray(fn, &sl);
122.         cur = integrateArray(sl);
123.     } while (fabs(cur - prev) > EPS2);
124.     return prev;
125. }
126.
127. int main(void) {
128.     double RANGES[3][2];
129.     RANGES[0][0] = -3.0;
130.     RANGES[0][1] = -1.0;
131.     RANGES[1][0] = -5.0;
132.     RANGES[1][1] = -3.0;
133.    // fprintf(stderr, "To send: %.4lf, %.4lf\n", RANGES[1][0], RANGES[1][1]);
134.     RANGES[2][0] = -1.0;
135.     RANGES[2][1] = -0.1;
136.     double cuts[3];
137.     //fprintf(stderr, "---\n");
138.     cuts[0] = root(1, 2, RANGES[0][0], RANGES[0][1]);
139.     cuts[1] = root(0, 2, RANGES[1][0], RANGES[1][1]);
140.     cuts[2] = root(0, 1, RANGES[2][0], RANGES[2][1]);
141.     //fprintf(stderr, "---\n");
142.     double ints[3];
143.     ints[0] = integrate(0, cuts[1], cuts[2]);
144.     ints[1] = integrate(1, cuts[0], cuts[2]);
145.     ints[2] = integrate(2, cuts[0], cuts[1]);
146.     //fprintf(stderr, "---\n");
147.     double result;
148.     if ((cuts[0] < cuts[1]) == (cuts[1] < cuts[2])) {
149.         result = ints[1] - ints[0] - ints[2];
150.     } else if ((cuts[0] < cuts[2]) == (cuts[2] < cuts[1])) {
151.         result = ints[2] - ints[0] - ints[1];
152.     } else {
153.         result = ints[0] - ints[1] - ints[2];
154.     }
155.     //fprintf(stderr, "---\n");
156.     printf("Result: %.4lf\n", fabs(result));
157.     printf("Points: %.4lf %.4lf %.4lf\n", cuts[0], cuts[1], cuts[2]);
158.     return 0;
159. }
clone this paste RAW Paste Data
Top