src/2geom/chebyshev.cpp

   1 #include <2geom/chebyshev.h>
   2
   3 #include <2geom/sbasis.h>
   4 #include <2geom/sbasis-poly.h>
   5
   6 #include <vector>
   7 using std::vector;
   8
   9 #include <gsl/gsl_math.h>
  10 #include <gsl/gsl_chebyshev.h>
  11
  12 namespace Geom{
  13
  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     }
  23
  24     return basis[n];
  25 }
  26
  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     }
  35
  36     return r;
  37 }
  38
  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     */
  45
  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     }
  55
  56     return y*d - dd + 0.5*cheb_coeff[0];
  57 }
  58
  59 SBasis chebyshev_approximant (double (*f)(double,void*), int order, Interval in, void* p) {
  60     gsl_cheb_series *cs = gsl_cheb_alloc (order+2);
  61
  62     gsl_function F;
  63
  64     F.function = f;
  65     F.params = p;
  66
  67     gsl_cheb_init (cs, &F, in[0], in[1]);
  68
  69     SBasis r = compose(clenshaw_series(order, cs->c), Linear(-1,1));
  70
  71     gsl_cheb_free (cs);
  72     return r;
  73 }
  74
  75 struct wrap {
  76     double (*f)(double,void*);
  77     void* pp;
  78     double fa, fb;
  79     Interval in;
  80 };
  81
  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 }
  87
  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 }
 101
 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     }
 111
 112     return basis[n];
 113 }
 114
 115 };
 116
 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:fileencoding=utf-8:textwidth=99 :