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