1 /*
2 * sbasis-math.cpp - some std functions to work with (pw)s-basis
3 *
4 * Authors:
5 * Jean-Francois Barraud
6 *
7 * Copyright (C) 2006-2007 authors
8 *
9 * This library is free software; you can redistribute it and/or
10 * modify it either under the terms of the GNU Lesser General Public
11 * License version 2.1 as published by the Free Software Foundation
12 * (the "LGPL") or, at your option, under the terms of the Mozilla
13 * Public License Version 1.1 (the "MPL"). If you do not alter this
14 * notice, a recipient may use your version of this file under either
15 * the MPL or the LGPL.
16 *
17 * You should have received a copy of the LGPL along with this library
18 * in the file COPYING-LGPL-2.1; if not, write to the Free Software
19 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
20 * You should have received a copy of the MPL along with this library
21 * in the file COPYING-MPL-1.1
22 *
23 * The contents of this file are subject to the Mozilla Public License
24 * Version 1.1 (the "License"); you may not use this file except in
25 * compliance with the License. You may obtain a copy of the License at
26 * http://www.mozilla.org/MPL/
27 *
28 * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
29 * OF ANY KIND, either express or implied. See the LGPL or the MPL for
30 * the specific language governing rights and limitations.
31 */
33 //this a first try to define sqrt, cos, sin, etc...
34 //TODO: define a truncated compose(sb,sb, order) and extend it to pw<sb>.
35 //TODO: in all these functions, compute 'order' according to 'tol'.
37 #include <2geom/sbasis-math.h>
39 #include <2geom/d2-sbasis.h>
40 #include <stdio.h>
41 #include <math.h>
42 //#define ZERO 1e-3
45 namespace Geom {
48 //-|x|-----------------------------------------------------------------------
49 /** Return the absolute value of a function pointwise.
50 \param f function
51 */
52 Piecewise<SBasis> abs(SBasis const &f){
53 return abs(Piecewise<SBasis>(f));
54 }
55 /** Return the absolute value of a function pointwise.
56 \param f function
57 */
58 Piecewise<SBasis> abs(Piecewise<SBasis> const &f){
59 Piecewise<SBasis> absf=partition(f,roots(f));
60 for (unsigned i=0; i<absf.size(); i++){
61 if (absf.segs[i](.5)<0) absf.segs[i]*=-1;
62 }
63 return absf;
64 }
66 //-max(x,y), min(x,y)--------------------------------------------------------
67 /** Return the greater of the two functions pointwise.
68 \param f, g two functions
69 */
70 Piecewise<SBasis> max( SBasis const &f, SBasis const &g){
71 return max(Piecewise<SBasis>(f),Piecewise<SBasis>(g));
72 }
73 /** Return the greater of the two functions pointwise.
74 \param f, g two functions
75 */
76 Piecewise<SBasis> max(Piecewise<SBasis> const &f, SBasis const &g){
77 return max(f,Piecewise<SBasis>(g));
78 }
79 /** Return the greater of the two functions pointwise.
80 \param f, g two functions
81 */
82 Piecewise<SBasis> max( SBasis const &f, Piecewise<SBasis> const &g){
83 return max(Piecewise<SBasis>(f),g);
84 }
85 /** Return the greater of the two functions pointwise.
86 \param f, g two functions
87 */
88 Piecewise<SBasis> max(Piecewise<SBasis> const &f, Piecewise<SBasis> const &g){
89 Piecewise<SBasis> max=partition(f,roots(f-g));
90 Piecewise<SBasis> gg =partition(g,max.cuts);
91 max = partition(max,gg.cuts);
92 for (unsigned i=0; i<max.size(); i++){
93 if (max.segs[i](.5)<gg.segs[i](.5)) max.segs[i]=gg.segs[i];
94 }
95 return max;
96 }
98 /** Return the more negative of the two functions pointwise.
99 \param f, g two functions
100 */
101 Piecewise<SBasis>
102 min( SBasis const &f, SBasis const &g){ return -max(-f,-g); }
103 /** Return the more negative of the two functions pointwise.
104 \param f, g two functions
105 */
106 Piecewise<SBasis>
107 min(Piecewise<SBasis> const &f, SBasis const &g){ return -max(-f,-g); }
108 /** Return the more negative of the two functions pointwise.
109 \param f, g two functions
110 */
111 Piecewise<SBasis>
112 min( SBasis const &f, Piecewise<SBasis> const &g){ return -max(-f,-g); }
113 /** Return the more negative of the two functions pointwise.
114 \param f, g two functions
115 */
116 Piecewise<SBasis>
117 min(Piecewise<SBasis> const &f, Piecewise<SBasis> const &g){ return -max(-f,-g); }
120 //-sign(x)---------------------------------------------------------------
121 /** Return the sign of the two functions pointwise.
122 \param f function
123 */
124 Piecewise<SBasis> signSb(SBasis const &f){
125 return signSb(Piecewise<SBasis>(f));
126 }
127 /** Return the sign of the two functions pointwise.
128 \param f function
129 */
130 Piecewise<SBasis> signSb(Piecewise<SBasis> const &f){
131 Piecewise<SBasis> sign=partition(f,roots(f));
132 for (unsigned i=0; i<sign.size(); i++){
133 sign.segs[i] = (sign.segs[i](.5)<0)? Linear(-1.):Linear(1.);
134 }
135 return sign;
136 }
138 //-Sqrt----------------------------------------------------------
139 static Piecewise<SBasis> sqrt_internal(SBasis const &f,
140 double tol,
141 int order){
142 SBasis sqrtf;
143 if(f.isZero() || order == 0){
144 return Piecewise<SBasis>(sqrtf);
145 }
146 if (f.at0()<-tol*tol && f.at1()<-tol*tol){
147 return sqrt_internal(-f,tol,order);
148 }else if (f.at0()>tol*tol && f.at1()>tol*tol){
149 sqrtf.resize(order+1, Linear(0,0));
150 sqrtf[0] = Linear(std::sqrt(f[0][0]), std::sqrt(f[0][1]));
151 SBasis r = f - multiply(sqrtf, sqrtf); // remainder
152 for(unsigned i = 1; int(i) <= order && i<r.size(); i++) {
153 Linear ci(r[i][0]/(2*sqrtf[0][0]), r[i][1]/(2*sqrtf[0][1]));
154 SBasis cisi = shift(ci, i);
155 r -= multiply(shift((sqrtf*2 + cisi), i), SBasis(ci));
156 r.truncate(order+1);
157 sqrtf[i] = ci;
158 if(r.tailError(i) == 0) // if exact
159 break;
160 }
161 }else{
162 sqrtf = Linear(std::sqrt(fabs(f.at0())), std::sqrt(fabs(f.at1())));
163 }
165 double err = (f - multiply(sqrtf, sqrtf)).tailError(0);
166 if (err<tol){
167 return Piecewise<SBasis>(sqrtf);
168 }
170 Piecewise<SBasis> sqrtf0,sqrtf1;
171 sqrtf0 = sqrt_internal(compose(f,Linear(0.,.5)),tol,order);
172 sqrtf1 = sqrt_internal(compose(f,Linear(.5,1.)),tol,order);
173 sqrtf0.setDomain(Interval(0.,.5));
174 sqrtf1.setDomain(Interval(.5,1.));
175 sqrtf0.concat(sqrtf1);
176 return sqrtf0;
177 }
179 /** Compute the sqrt of a function.
180 \param f function
181 */
182 Piecewise<SBasis> sqrt(SBasis const &f, double tol, int order){
183 return sqrt(max(f,Linear(tol*tol)),tol,order);
184 }
186 /** Compute the sqrt of a function.
187 \param f function
188 */
189 Piecewise<SBasis> sqrt(Piecewise<SBasis> const &f, double tol, int order){
190 Piecewise<SBasis> result;
191 Piecewise<SBasis> zero = Piecewise<SBasis>(Linear(tol*tol));
192 zero.setDomain(f.domain());
193 Piecewise<SBasis> ff=max(f,zero);
195 for (unsigned i=0; i<ff.size(); i++){
196 Piecewise<SBasis> sqrtfi = sqrt_internal(ff.segs[i],tol,order);
197 sqrtfi.setDomain(Interval(ff.cuts[i],ff.cuts[i+1]));
198 result.concat(sqrtfi);
199 }
200 return result;
201 }
203 //-Yet another sin/cos--------------------------------------------------------------
205 /** Compute the sine of a function.
206 \param f function
207 \param tol maximum error
208 \param order maximum degree polynomial to use
209 */
210 Piecewise<SBasis> sin( SBasis const &f, double tol, int order){return(cos(-f+M_PI/2,tol,order));}
211 /** Compute the sine of a function.
212 \param f function
213 \param tol maximum error
214 \param order maximum degree polynomial to use
215 */
216 Piecewise<SBasis> sin(Piecewise<SBasis> const &f, double tol, int order){return(cos(-f+M_PI/2,tol,order));}
218 /** Compute the cosine of a function.
219 \param f function
220 \param tol maximum error
221 \param order maximum degree polynomial to use
222 */
223 Piecewise<SBasis> cos(Piecewise<SBasis> const &f, double tol, int order){
224 Piecewise<SBasis> result;
225 for (unsigned i=0; i<f.size(); i++){
226 Piecewise<SBasis> cosfi = cos(f.segs[i],tol,order);
227 cosfi.setDomain(Interval(f.cuts[i],f.cuts[i+1]));
228 result.concat(cosfi);
229 }
230 return result;
231 }
233 /** Compute the cosine of a function.
234 \param f function
235 \param tol maximum error
236 \param order maximum degree polynomial to use
237 */
238 Piecewise<SBasis> cos( SBasis const &f, double tol, int order){
239 double alpha = (f.at0()+f.at1())/2.;
240 SBasis x = f-alpha;
241 double d = x.tailError(0),err=1;
242 //estimate cos(x)-sum_0^order (-1)^k x^2k/2k! by the first neglicted term
243 for (int i=1; i<=2*order; i++) err*=d/i;
245 if (err<tol){
246 SBasis xk=Linear(1), c=Linear(1), s=Linear(0);
247 for (int k=1; k<=2*order; k+=2){
248 xk*=x/k;
249 //take also truncature errors into account...
250 err+=xk.tailError(order);
251 xk.truncate(order);
252 s+=xk;
253 xk*=-x/(k+1);
254 //take also truncature errors into account...
255 err+=xk.tailError(order);
256 xk.truncate(order);
257 c+=xk;
258 }
259 if (err<tol){
260 return Piecewise<SBasis>(std::cos(alpha)*c-std::sin(alpha)*s);
261 }
262 }
263 Piecewise<SBasis> c0,c1;
264 c0 = cos(compose(f,Linear(0.,.5)),tol,order);
265 c1 = cos(compose(f,Linear(.5,1.)),tol,order);
266 c0.setDomain(Interval(0.,.5));
267 c1.setDomain(Interval(.5,1.));
268 c0.concat(c1);
269 return c0;
270 }
272 //--1/x------------------------------------------------------------
273 //TODO: this implementation is just wrong. Remove or redo!
275 void truncateResult(Piecewise<SBasis> &f, int order){
276 if (order>=0){
277 for (unsigned k=0; k<f.segs.size(); k++){
278 f.segs[k].truncate(order);
279 }
280 }
281 }
283 Piecewise<SBasis> reciprocalOnDomain(Interval range, double tol){
284 Piecewise<SBasis> reciprocal_fn;
285 //TODO: deduce R from tol...
286 double R=2.;
287 SBasis reciprocal1_R=reciprocal(Linear(1,R),3);
288 double a=range.min(), b=range.max();
289 if (a*b<0){
290 b=std::max(fabs(a),fabs(b));
291 a=0;
292 }else if (b<0){
293 a=-range.max();
294 b=-range.min();
295 }
297 if (a<=tol){
298 reciprocal_fn.push_cut(0);
299 int i0=(int) floor(std::log(tol)/std::log(R));
300 a=pow(R,i0);
301 reciprocal_fn.push(Linear(1/a),a);
302 }else{
303 int i0=(int) floor(std::log(a)/std::log(R));
304 a=pow(R,i0);
305 reciprocal_fn.cuts.push_back(a);
306 }
308 while (a<b){
309 reciprocal_fn.push(reciprocal1_R/a,R*a);
310 a*=R;
311 }
312 if (range.min()<0 || range.max()<0){
313 Piecewise<SBasis>reciprocal_fn_neg;
314 //TODO: define reverse(pw<sb>);
315 reciprocal_fn_neg.cuts.push_back(-reciprocal_fn.cuts.back());
316 for (unsigned i=0; i<reciprocal_fn.size(); i++){
317 int idx=reciprocal_fn.segs.size()-1-i;
318 reciprocal_fn_neg.push_seg(-reverse(reciprocal_fn.segs.at(idx)));
319 reciprocal_fn_neg.push_cut(-reciprocal_fn.cuts.at(idx));
320 }
321 if (range.max()>0){
322 reciprocal_fn_neg.concat(reciprocal_fn);
323 }
324 reciprocal_fn=reciprocal_fn_neg;
325 }
327 return(reciprocal_fn);
328 }
330 Piecewise<SBasis> reciprocal(SBasis const &f, double tol, int order){
331 Piecewise<SBasis> reciprocal_fn=reciprocalOnDomain(*bounds_fast(f), tol);
332 Piecewise<SBasis> result=compose(reciprocal_fn,f);
333 truncateResult(result,order);
334 return(result);
335 }
336 Piecewise<SBasis> reciprocal(Piecewise<SBasis> const &f, double tol, int order){
337 Piecewise<SBasis> reciprocal_fn=reciprocalOnDomain(*bounds_fast(f), tol);
338 Piecewise<SBasis> result=compose(reciprocal_fn,f);
339 truncateResult(result,order);
340 return(result);
341 }
343 /**
344 * \brief Retruns a Piecewise SBasis with prescribed values at prescribed times.
345 *
346 * \param times: vector of times at which the values are given. Should be sorted in increasing order.
347 * \param values: vector of prescribed values. Should have the same size as times and be sorted accordingly.
348 * \param smoothness: (defaults to 1) regularity class of the result: 0=piecewise linear, 1=continuous derivative, etc...
349 */
350 Piecewise<SBasis> interpolate(std::vector<double> times, std::vector<double> values, unsigned smoothness){
351 assert ( values.size() == times.size() );
352 if ( values.size() == 0 ) return Piecewise<SBasis>();
353 if ( values.size() == 1 ) return Piecewise<SBasis>(values[0]);//what about time??
355 SBasis sk = shift(Linear(1.),smoothness);
356 SBasis bump_in = integral(sk);
357 bump_in -= bump_in.at0();
358 bump_in /= bump_in.at1();
359 SBasis bump_out = reverse( bump_in );
361 Piecewise<SBasis> result;
362 result.cuts.push_back(times[0]);
363 for (unsigned i = 0; i<values.size()-1; i++){
364 result.push(bump_out*values[i]+bump_in*values[i+1],times[i+1]);
365 }
366 return result;
367 }
369 }
371 /*
372 Local Variables:
373 mode:c++
374 c-file-style:"stroustrup"
375 c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
376 indent-tabs-mode:nil
377 fill-column:99
378 End:
379 */
380 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :