summaryrefslogtreecommitdiffstats
path: root/src/2geom/chebyshev.cpp
diff options
context:
space:
mode:
Diffstat (limited to 'src/2geom/chebyshev.cpp')
-rw-r--r--src/2geom/chebyshev.cpp252
1 files changed, 126 insertions, 126 deletions
diff --git a/src/2geom/chebyshev.cpp b/src/2geom/chebyshev.cpp
index 447c5183f..4ba3e2325 100644
--- a/src/2geom/chebyshev.cpp
+++ b/src/2geom/chebyshev.cpp
@@ -1,126 +1,126 @@
-#include <2geom/chebyshev.h>
-
-#include <2geom/sbasis.h>
-#include <2geom/sbasis-poly.h>
-
-#include <vector>
-using std::vector;
-
-#include <gsl/gsl_math.h>
-#include <gsl/gsl_chebyshev.h>
-
-namespace Geom{
-
-SBasis cheb(unsigned n) {
- static std::vector<SBasis> basis;
- if(basis.empty()) {
- basis.push_back(Linear(1,1));
- basis.push_back(Linear(0,1));
- }
- for(unsigned i = basis.size(); i <= n; i++) {
- basis.push_back(Linear(0,2)*basis[i-1] - basis[i-2]);
- }
-
- return basis[n];
-}
-
-SBasis cheb_series(unsigned n, double* cheb_coeff) {
- SBasis r;
- for(unsigned i = 0; i < n; i++) {
- double cof = cheb_coeff[i];
- //if(i == 0)
- //cof /= 2;
- r += cheb(i)*cof;
- }
-
- return r;
-}
-
-SBasis clenshaw_series(unsigned m, double* cheb_coeff) {
- /** b_n = a_n
- b_n-1 = 2*x*b_n + a_n-1
- b_n-k = 2*x*b_{n-k+1} + a_{n-k} - b_{n - k + 2}
- b_0 = x*b_1 + a_0 - b_2
- */
-
- double a = -1, b = 1;
- SBasis d, dd;
- SBasis y = (Linear(0, 2) - (a+b)) / (b-a);
- SBasis y2 = 2*y;
- for(int j = m - 1; j >= 1; j--) {
- SBasis sv = d;
- d = y2*d - dd + cheb_coeff[j];
- dd = sv;
- }
-
- return y*d - dd + 0.5*cheb_coeff[0];
-}
-
-SBasis chebyshev_approximant (double (*f)(double,void*), int order, Interval in, void* p) {
- gsl_cheb_series *cs = gsl_cheb_alloc (order+2);
-
- gsl_function F;
-
- F.function = f;
- F.params = p;
-
- gsl_cheb_init (cs, &F, in[0], in[1]);
-
- SBasis r = compose(clenshaw_series(order, cs->c), Linear(-1,1));
-
- gsl_cheb_free (cs);
- return r;
-}
-
-struct wrap {
- double (*f)(double,void*);
- void* pp;
- double fa, fb;
- Interval in;
-};
-
-double f_interp(double x, void* p) {
- struct wrap *wr = (struct wrap *)p;
- double z = (x - wr->in[0]) / (wr->in[1] - wr->in[0]);
- return (wr->f)(x, wr->pp) - ((1 - z)*wr->fa + z*wr->fb);
-}
-
-SBasis chebyshev_approximant_interpolating (double (*f)(double,void*),
- int order, Interval in, void* p) {
- double fa = f(in[0], p);
- double fb = f(in[1], p);
- struct wrap wr;
- wr.fa = fa;
- wr.fb = fb;
- wr.in = in;
- printf("%f %f\n", fa, fb);
- wr.f = f;
- wr.pp = p;
- return compose(Linear(in[0], in[1]), Linear(fa, fb)) + chebyshev_approximant(f_interp, order, in, &wr) + Linear(fa, fb);
-}
-
-SBasis chebyshev(unsigned n) {
- static std::vector<SBasis> basis;
- if(basis.empty()) {
- basis.push_back(Linear(1,1));
- basis.push_back(Linear(0,1));
- }
- for(unsigned i = basis.size(); i <= n; i++) {
- basis.push_back(Linear(0,2)*basis[i-1] - basis[i-2]);
- }
-
- return basis[n];
-}
-
-};
-
-/*
- Local Variables:
- mode:c++
- c-file-style:"stroustrup"
- c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
- indent-tabs-mode:nil
- fill-column:99
- End:
-*/
-// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :
+#include <2geom/chebyshev.h>
+
+#include <2geom/sbasis.h>
+#include <2geom/sbasis-poly.h>
+
+#include <vector>
+using std::vector;
+
+#include <gsl/gsl_math.h>
+#include <gsl/gsl_chebyshev.h>
+
+namespace Geom{
+
+SBasis cheb(unsigned n) {
+ static std::vector<SBasis> basis;
+ if(basis.empty()) {
+ basis.push_back(Linear(1,1));
+ basis.push_back(Linear(0,1));
+ }
+ for(unsigned i = basis.size(); i <= n; i++) {
+ basis.push_back(Linear(0,2)*basis[i-1] - basis[i-2]);
+ }
+
+ return basis[n];
+}
+
+SBasis cheb_series(unsigned n, double* cheb_coeff) {
+ SBasis r;
+ for(unsigned i = 0; i < n; i++) {
+ double cof = cheb_coeff[i];
+ //if(i == 0)
+ //cof /= 2;
+ r += cheb(i)*cof;
+ }
+
+ return r;
+}
+
+SBasis clenshaw_series(unsigned m, double* cheb_coeff) {
+ /** b_n = a_n
+ b_n-1 = 2*x*b_n + a_n-1
+ b_n-k = 2*x*b_{n-k+1} + a_{n-k} - b_{n - k + 2}
+ b_0 = x*b_1 + a_0 - b_2
+ */
+
+ double a = -1, b = 1;
+ SBasis d, dd;
+ SBasis y = (Linear(0, 2) - (a+b)) / (b-a);
+ SBasis y2 = 2*y;
+ for(int j = m - 1; j >= 1; j--) {
+ SBasis sv = d;
+ d = y2*d - dd + cheb_coeff[j];
+ dd = sv;
+ }
+
+ return y*d - dd + 0.5*cheb_coeff[0];
+}
+
+SBasis chebyshev_approximant (double (*f)(double,void*), int order, Interval in, void* p) {
+ gsl_cheb_series *cs = gsl_cheb_alloc (order+2);
+
+ gsl_function F;
+
+ F.function = f;
+ F.params = p;
+
+ gsl_cheb_init (cs, &F, in[0], in[1]);
+
+ SBasis r = compose(clenshaw_series(order, cs->c), Linear(-1,1));
+
+ gsl_cheb_free (cs);
+ return r;
+}
+
+struct wrap {
+ double (*f)(double,void*);
+ void* pp;
+ double fa, fb;
+ Interval in;
+};
+
+double f_interp(double x, void* p) {
+ struct wrap *wr = (struct wrap *)p;
+ double z = (x - wr->in[0]) / (wr->in[1] - wr->in[0]);
+ return (wr->f)(x, wr->pp) - ((1 - z)*wr->fa + z*wr->fb);
+}
+
+SBasis chebyshev_approximant_interpolating (double (*f)(double,void*),
+ int order, Interval in, void* p) {
+ double fa = f(in[0], p);
+ double fb = f(in[1], p);
+ struct wrap wr;
+ wr.fa = fa;
+ wr.fb = fb;
+ wr.in = in;
+ printf("%f %f\n", fa, fb);
+ wr.f = f;
+ wr.pp = p;
+ return compose(Linear(in[0], in[1]), Linear(fa, fb)) + chebyshev_approximant(f_interp, order, in, &wr) + Linear(fa, fb);
+}
+
+SBasis chebyshev(unsigned n) {
+ static std::vector<SBasis> basis;
+ if(basis.empty()) {
+ basis.push_back(Linear(1,1));
+ basis.push_back(Linear(0,1));
+ }
+ for(unsigned i = basis.size(); i <= n; i++) {
+ basis.push_back(Linear(0,2)*basis[i-1] - basis[i-2]);
+ }
+
+ return basis[n];
+}
+
+};
+
+/*
+ Local Variables:
+ mode:c++
+ c-file-style:"stroustrup"
+ c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
+ indent-tabs-mode:nil
+ fill-column:99
+ End:
+*/
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :
generated by cgit v1.2.3 (git 2.25.1) at 2026-05-19 09:26:45 +0000