Subversion Repositories Scribus

Rev

Rev 15347 | Details | Compare with Previous | Last modification | View Log | RSS feed

Rev Author Line No. Line
15347 fschmid 1
/*
2
 *  sbasis.h - S-power basis function class
12065 fschmid 3
 *
4
 *  Authors:
5
 *   Nathan Hurst <njh@mail.csse.monash.edu.au>
6
 *   Michael Sloan <mgsloan@gmail.com>
7
 *
8
 * Copyright (C) 2006-2007 authors
9
 *
10
 * This library is free software; you can redistribute it and/or
11
 * modify it either under the terms of the GNU Lesser General Public
12
 * License version 2.1 as published by the Free Software Foundation
13
 * (the "LGPL") or, at your option, under the terms of the Mozilla
14
 * Public License Version 1.1 (the "MPL"). If you do not alter this
15
 * notice, a recipient may use your version of this file under either
16
 * the MPL or the LGPL.
17
 *
18
 * You should have received a copy of the LGPL along with this library
19
 * in the file COPYING-LGPL-2.1; if not, write to the Free Software
18122 mrdocs 20
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
12065 fschmid 21
 * You should have received a copy of the MPL along with this library
22
 * in the file COPYING-MPL-1.1
23
 *
24
 * The contents of this file are subject to the Mozilla Public License
25
 * Version 1.1 (the "License"); you may not use this file except in
26
 * compliance with the License. You may obtain a copy of the License at
27
 * http://www.mozilla.org/MPL/
28
 *
29
 * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
30
 * OF ANY KIND, either express or implied. See the LGPL or the MPL for
31
 * the specific language governing rights and limitations.
32
 */
33
 
34
#ifndef SEEN_SBASIS_H
35
#define SEEN_SBASIS_H
36
#include <vector>
37
#include <cassert>
38
#include <iostream>
39
 
40
#include "linear.h"
41
#include "interval.h"
42
#include "utils.h"
43
#include "exception.h"
44
 
15347 fschmid 45
namespace Geom {
12065 fschmid 46
 
47
/*** An empty SBasis is identically 0. */
15347 fschmid 48
class SBasis : public std::vector<Linear>{
12065 fschmid 49
public:
50
    SBasis() {}
51
    explicit SBasis(double a) {
52
        push_back(Linear(a,a));
53
    }
54
    SBasis(SBasis const & a) :
15347 fschmid 55
        std::vector<Linear>(a)
12065 fschmid 56
    {}
57
    SBasis(Linear const & bo) {
58
        push_back(bo);
59
    }
60
 
61
    //IMPL: FragmentConcept
62
    typedef double output_type;
63
    inline bool isZero() const {
64
        if(empty()) return true;
65
        for(unsigned i = 0; i < size(); i++) {
66
            if(!(*this)[i].isZero()) return false;
67
        }
68
        return true;
69
    }
70
    inline bool isConstant() const {
71
        if (empty()) return true;
15347 fschmid 72
        for (unsigned i = 0; i < size(); i++) {
73
            if(!(*this)[i].isConstant()) return false;
12065 fschmid 74
        }
75
        return true;
76
    }
77
 
78
    bool isFinite() const;
79
    inline double at0() const {
80
        if(empty()) return 0; else return (*this)[0][0];
81
    }
82
    inline double at1() const{
83
        if(empty()) return 0; else return (*this)[0][1];
84
    }
85
 
86
    double valueAt(double t) const {
87
        double s = t*(1-t);
88
        double p0 = 0, p1 = 0;
15347 fschmid 89
        double sk = 1;
90
//TODO: rewrite as horner
91
        for(unsigned k = 0; k < size(); k++) {
92
            p0 += sk*(*this)[k][0];
93
            p1 += sk*(*this)[k][1];
94
            sk *= s;
12065 fschmid 95
        }
96
        return (1-t)*p0 + t*p1;
97
    }
15347 fschmid 98
    double valueAndDerivative(double t, double &der) const {
99
        double s = t*(1-t);
100
        double p0 = 0, p1 = 0;
101
        double sk = 1;
102
//TODO: rewrite as horner
103
        for(unsigned k = 0; k < size(); k++) {
104
            p0 += sk*(*this)[k][0];
105
            p1 += sk*(*this)[k][1];
106
            sk *= s;
107
        }
108
        // p0 and p1 at this point form a linear approximation at t
109
        der = p1 - p0;
110
        return (1-t)*p0 + t*p1;
111
    }
12065 fschmid 112
    double operator()(double t) const {
113
        return valueAt(t);
114
    }
115
 
15347 fschmid 116
    std::vector<double> valueAndDerivatives(double /*t*/, unsigned /*n*/) const {
117
        //TODO
118
        throwNotImplemented(0);
119
    }
12065 fschmid 120
 
121
    SBasis toSBasis() const { return SBasis(*this); }
122
 
123
    double tailError(unsigned tail) const;
124
 
125
// compute f(g)
126
    SBasis operator()(SBasis const & g) const;
127
 
15347 fschmid 128
    Linear operator[](unsigned i) const {
129
        assert(i < size());
130
        return std::vector<Linear>::operator[](i);
131
    }
132
 
12065 fschmid 133
//MUTATOR PRISON
15347 fschmid 134
    Linear& operator[](unsigned i) { return this->at(i); }
135
 
12065 fschmid 136
    //remove extra zeros
137
    void normalize() {
138
        while(!empty() && 0 == back()[0] && 0 == back()[1])
139
            pop_back();
140
    }
141
    void truncate(unsigned k) { if(k < size()) resize(k); }
142
};
143
 
144
//TODO: figure out how to stick this in linear, while not adding an sbasis dep
145
inline SBasis Linear::toSBasis() const { return SBasis(*this); }
146
 
147
//implemented in sbasis-roots.cpp
15347 fschmid 148
Interval bounds_exact(SBasis const &a);
149
Interval bounds_fast(SBasis const &a, int order = 0);
150
Interval bounds_local(SBasis const &a, const Interval &t, int order = 0);
12065 fschmid 151
 
152
inline SBasis reverse(SBasis const &a) {
15347 fschmid 153
    SBasis result;
154
    result.reserve(a.size());
12065 fschmid 155
    for(unsigned k = 0; k < a.size(); k++)
15347 fschmid 156
       result.push_back(reverse(a[k]));
12065 fschmid 157
    return result;
158
}
159
 
160
//IMPL: ScalableConcept
161
inline SBasis operator-(const SBasis& p) {
162
    if(p.isZero()) return SBasis();
15347 fschmid 163
    SBasis result;
164
    result.reserve(p.size());
12065 fschmid 165
 
166
    for(unsigned i = 0; i < p.size(); i++) {
15347 fschmid 167
        result.push_back(-p[i]);
12065 fschmid 168
    }
169
    return result;
170
}
171
SBasis operator*(SBasis const &a, double k);
172
inline SBasis operator*(double k, SBasis const &a) { return a*k; }
173
inline SBasis operator/(SBasis const &a, double k) { return a*(1./k); }
174
SBasis& operator*=(SBasis& a, double b);
175
inline SBasis& operator/=(SBasis& a, double b) { return (a*=(1./b)); }
176
 
177
//IMPL: AddableConcept
178
SBasis operator+(const SBasis& a, const SBasis& b);
179
SBasis operator-(const SBasis& a, const SBasis& b);
180
SBasis& operator+=(SBasis& a, const SBasis& b);
181
SBasis& operator-=(SBasis& a, const SBasis& b);
182
 
183
//TODO: remove?
15347 fschmid 184
inline SBasis operator+(const SBasis & a, Linear const & b) {
12065 fschmid 185
    if(b.isZero()) return a;
186
    if(a.isZero()) return b;
187
    SBasis result(a);
188
    result[0] += b;
189
    return result;
190
}
191
inline SBasis operator-(const SBasis & a, Linear const & b) {
192
    if(b.isZero()) return a;
193
    SBasis result(a);
194
    result[0] -= b;
195
    return result;
196
}
197
inline SBasis& operator+=(SBasis& a, const Linear& b) {
198
    if(a.isZero())
199
        a.push_back(b);
200
    else
201
        a[0] += b;
202
    return a;
203
}
204
inline SBasis& operator-=(SBasis& a, const Linear& b) {
205
    if(a.isZero())
206
        a.push_back(-b);
207
    else
208
        a[0] -= b;
209
    return a;
15347 fschmid 210
}
12065 fschmid 211
 
212
//IMPL: OffsetableConcept
213
inline SBasis operator+(const SBasis & a, double b) {
214
    if(a.isZero()) return Linear(b, b);
215
    SBasis result(a);
216
    result[0] += b;
217
    return result;
218
}
219
inline SBasis operator-(const SBasis & a, double b) {
220
    if(a.isZero()) return Linear(-b, -b);
221
    SBasis result(a);
222
    result[0] -= b;
223
    return result;
224
}
225
inline SBasis& operator+=(SBasis& a, double b) {
226
    if(a.isZero())
15347 fschmid 227
        a.push_back(Linear(b,b));
12065 fschmid 228
    else
229
        a[0] += b;
230
    return a;
231
}
232
inline SBasis& operator-=(SBasis& a, double b) {
233
    if(a.isZero())
15347 fschmid 234
        a.push_back(Linear(-b,-b));
12065 fschmid 235
    else
236
        a[0] -= b;
237
    return a;
238
}
239
 
240
SBasis shift(SBasis const &a, int sh);
241
SBasis shift(Linear const &a, int sh);
242
 
243
inline SBasis truncate(SBasis const &a, unsigned terms) {
244
    SBasis c;
245
    c.insert(c.begin(), a.begin(), a.begin() + std::min(terms, (unsigned)a.size()));
246
    return c;
247
}
248
 
249
SBasis multiply(SBasis const &a, SBasis const &b);
250
 
251
SBasis integral(SBasis const &c);
252
SBasis derivative(SBasis const &a);
253
 
254
SBasis sqrt(SBasis const &a, int k);
255
 
256
// return a kth order approx to 1/a)
257
SBasis reciprocal(Linear const &a, int k);
258
SBasis divide(SBasis const &a, SBasis const &b, int k);
259
 
260
inline SBasis operator*(SBasis const & a, SBasis const & b) {
261
    return multiply(a, b);
262
}
263
 
264
inline SBasis& operator*=(SBasis& a, SBasis const & b) {
265
    a = multiply(a, b);
266
    return a;
267
}
268
 
15347 fschmid 269
//valuation: degree of the first non zero coefficient.
12065 fschmid 270
inline unsigned
271
valuation(SBasis const &a, double tol=0){
272
    unsigned val=0;
273
    while( val<a.size() &&
274
           fabs(a[val][0])<tol &&
275
           fabs(a[val][1])<tol )
276
        val++;
277
    return val;
278
}
279
 
280
// a(b(t))
281
SBasis compose(SBasis const &a, SBasis const &b);
282
SBasis compose(SBasis const &a, SBasis const &b, unsigned k);
283
SBasis inverse(SBasis a, int k);
284
//compose_inverse(f,g)=compose(f,inverse(g)), but is numerically more stable in some good cases...
285
//TODO: requires g(0)=0 & g(1)=1 atm. generalization should be obvious.
286
SBasis compose_inverse(SBasis const &f, SBasis const &g, unsigned order=2, double tol=1e-3);
287
 
288
inline SBasis portion(const SBasis &t, double from, double to) { return compose(t, Linear(from, to)); }
289
 
290
// compute f(g)
291
inline SBasis
292
SBasis::operator()(SBasis const & g) const {
293
    return compose(*this, g);
294
}
295
 
296
inline std::ostream &operator<< (std::ostream &out_file, const Linear &bo) {
297
    out_file << "{" << bo[0] << ", " << bo[1] << "}";
298
    return out_file;
299
}
300
 
301
inline std::ostream &operator<< (std::ostream &out_file, const SBasis & p) {
302
    for(unsigned i = 0; i < p.size(); i++) {
303
        out_file << p[i] << "s^" << i << " + ";
304
    }
305
    return out_file;
306
}
307
 
308
SBasis sin(Linear bo, int k);
309
SBasis cos(Linear bo, int k);
310
 
311
std::vector<double> roots(SBasis const & s);
312
std::vector<std::vector<double> > multi_roots(SBasis const &f,
313
                                 std::vector<double> const &levels,
314
                                 double htol=1e-7,
315
                                 double vtol=1e-7,
316
                                 double a=0,
317
                                 double b=1);
15347 fschmid 318
 
12065 fschmid 319
}
320
 
321
/*
322
  Local Variables:
323
  mode:c++
324
  c-file-style:"stroustrup"
325
  c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
326
  indent-tabs-mode:nil
327
  fill-column:99
328
  End:
329
*/
330
// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :
331
#endif