1 #include <2geom/chebyshev.h>
3 #include <2geom/sbasis.h>
4 #include <2geom/sbasis-poly.h>
6 #include <vector>
7 using std::vector;
9 #include <gsl/gsl_math.h>
10 #include <gsl/gsl_chebyshev.h>
12 namespace Geom{
14 SBasis cheb(unsigned n) {
15 static std::vector<SBasis> basis;
16 if(basis.empty()) {
17 basis.push_back(Linear(1,1));
18 basis.push_back(Linear(0,1));
19 }
20 for(unsigned i = basis.size(); i <= n; i++) {
21 basis.push_back(Linear(0,2)*basis[i-1] - basis[i-2]);
22 }
24 return basis[n];
25 }
27 SBasis cheb_series(unsigned n, double* cheb_coeff) {
28 SBasis r;
29 for(unsigned i = 0; i < n; i++) {
30 double cof = cheb_coeff[i];
31 //if(i == 0)
32 //cof /= 2;
33 r += cheb(i)*cof;
34 }
36 return r;
37 }
39 SBasis clenshaw_series(unsigned m, double* cheb_coeff) {
40 /** b_n = a_n
41 b_n-1 = 2*x*b_n + a_n-1
42 b_n-k = 2*x*b_{n-k+1} + a_{n-k} - b_{n - k + 2}
43 b_0 = x*b_1 + a_0 - b_2
44 */
46 double a = -1, b = 1;
47 SBasis d, dd;
48 SBasis y = (Linear(0, 2) - (a+b)) / (b-a);
49 SBasis y2 = 2*y;
50 for(int j = m - 1; j >= 1; j--) {
51 SBasis sv = d;
52 d = y2*d - dd + cheb_coeff[j];
53 dd = sv;
54 }
56 return y*d - dd + 0.5*cheb_coeff[0];
57 }
59 SBasis chebyshev_approximant (double (*f)(double,void*), int order, Interval in, void* p) {
60 gsl_cheb_series *cs = gsl_cheb_alloc (order+2);
62 gsl_function F;
64 F.function = f;
65 F.params = p;
67 gsl_cheb_init (cs, &F, in[0], in[1]);
69 SBasis r = compose(clenshaw_series(order, cs->c), Linear(-1,1));
71 gsl_cheb_free (cs);
72 return r;
73 }
75 struct wrap {
76 double (*f)(double,void*);
77 void* pp;
78 double fa, fb;
79 Interval in;
80 };
82 double f_interp(double x, void* p) {
83 struct wrap *wr = (struct wrap *)p;
84 double z = (x - wr->in[0]) / (wr->in[1] - wr->in[0]);
85 return (wr->f)(x, wr->pp) - ((1 - z)*wr->fa + z*wr->fb);
86 }
88 SBasis chebyshev_approximant_interpolating (double (*f)(double,void*),
89 int order, Interval in, void* p) {
90 double fa = f(in[0], p);
91 double fb = f(in[1], p);
92 struct wrap wr;
93 wr.fa = fa;
94 wr.fb = fb;
95 wr.in = in;
96 printf("%f %f\n", fa, fb);
97 wr.f = f;
98 wr.pp = p;
99 return compose(Linear(in[0], in[1]), Linear(fa, fb)) + chebyshev_approximant(f_interp, order, in, &wr) + Linear(fa, fb);
100 }
102 SBasis chebyshev(unsigned n) {
103 static std::vector<SBasis> basis;
104 if(basis.empty()) {
105 basis.push_back(Linear(1,1));
106 basis.push_back(Linear(0,1));
107 }
108 for(unsigned i = basis.size(); i <= n; i++) {
109 basis.push_back(Linear(0,2)*basis[i-1] - basis[i-2]);
110 }
112 return basis[n];
113 }
115 };
117 /*
118 Local Variables:
119 mode:c++
120 c-file-style:"stroustrup"
121 c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
122 indent-tabs-mode:nil
123 fill-column:99
124 End:
125 */
126 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :