1d179a5630ec2dbc7436b536eecf1447b383a47d
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>
38 //#define ZERO 1e-3
41 namespace Geom {
44 #include <stdio.h>
45 #include <math.h>
47 //-|x|-----------------------------------------------------------------------
48 /** Return the absolute value of a function pointwise.
49 \param f function
50 */
51 Piecewise<SBasis> abs(SBasis const &f){
52 return abs(Piecewise<SBasis>(f));
53 }
54 /** Return the absolute value of a function pointwise.
55 \param f function
56 */
57 Piecewise<SBasis> abs(Piecewise<SBasis> const &f){
58 Piecewise<SBasis> absf=partition(f,roots(f));
59 for (unsigned i=0; i<absf.size(); i++){
60 if (absf.segs[i](.5)<0) absf.segs[i]*=-1;
61 }
62 return absf;
63 }
65 //-max(x,y), min(x,y)--------------------------------------------------------
66 /** Return the greater of the two functions pointwise.
67 \param f, g two functions
68 */
69 Piecewise<SBasis> max( SBasis const &f, SBasis const &g){
70 return max(Piecewise<SBasis>(f),Piecewise<SBasis>(g));
71 }
72 /** Return the greater of the two functions pointwise.
73 \param f, g two functions
74 */
75 Piecewise<SBasis> max(Piecewise<SBasis> const &f, SBasis const &g){
76 return max(f,Piecewise<SBasis>(g));
77 }
78 /** Return the greater of the two functions pointwise.
79 \param f, g two functions
80 */
81 Piecewise<SBasis> max( SBasis const &f, Piecewise<SBasis> const &g){
82 return max(Piecewise<SBasis>(f),g);
83 }
84 /** Return the greater of the two functions pointwise.
85 \param f, g two functions
86 */
87 Piecewise<SBasis> max(Piecewise<SBasis> const &f, Piecewise<SBasis> const &g){
88 Piecewise<SBasis> max=partition(f,roots(f-g));
89 Piecewise<SBasis> gg =partition(g,max.cuts);
90 max = partition(max,gg.cuts);
91 for (unsigned i=0; i<max.size(); i++){
92 if (max.segs[i](.5)<gg.segs[i](.5)) max.segs[i]=gg.segs[i];
93 }
94 return max;
95 }
97 /** Return the more negative of the two functions pointwise.
98 \param f, g two functions
99 */
100 Piecewise<SBasis>
101 min( SBasis const &f, SBasis const &g){ return -max(-f,-g); }
102 /** Return the more negative of the two functions pointwise.
103 \param f, g two functions
104 */
105 Piecewise<SBasis>
106 min(Piecewise<SBasis> const &f, SBasis const &g){ return -max(-f,-g); }
107 /** Return the more negative of the two functions pointwise.
108 \param f, g two functions
109 */
110 Piecewise<SBasis>
111 min( SBasis const &f, Piecewise<SBasis> const &g){ return -max(-f,-g); }
112 /** Return the more negative of the two functions pointwise.
113 \param f, g two functions
114 */
115 Piecewise<SBasis>
116 min(Piecewise<SBasis> const &f, Piecewise<SBasis> const &g){ return -max(-f,-g); }
119 //-sign(x)---------------------------------------------------------------
120 /** Return the sign of the two functions pointwise.
121 \param f function
122 */
123 Piecewise<SBasis> signSb(SBasis const &f){
124 return signSb(Piecewise<SBasis>(f));
125 }
126 /** Return the sign of the two functions pointwise.
127 \param f function
128 */
129 Piecewise<SBasis> signSb(Piecewise<SBasis> const &f){
130 Piecewise<SBasis> sign=partition(f,roots(f));
131 for (unsigned i=0; i<sign.size(); i++){
132 sign.segs[i] = (sign.segs[i](.5)<0)? Linear(-1.):Linear(1.);
133 }
134 return sign;
135 }
137 //-Sqrt----------------------------------------------------------
138 static Piecewise<SBasis> sqrt_internal(SBasis const &f,
139 double tol,
140 int order){
141 SBasis sqrtf;
142 if(f.isZero() || order == 0){
143 return Piecewise<SBasis>(sqrtf);
144 }
145 if (f.at0()<-tol*tol && f.at1()<-tol*tol){
146 return sqrt_internal(-f,tol,order);
147 }else if (f.at0()>tol*tol && f.at1()>tol*tol){
148 sqrtf.resize(order+1, Linear(0,0));
149 sqrtf[0] = Linear(std::sqrt(f[0][0]), std::sqrt(f[0][1]));
150 SBasis r = f - multiply(sqrtf, sqrtf); // remainder
151 for(unsigned i = 1; int(i) <= order and i<r.size(); i++) {
152 Linear ci(r[i][0]/(2*sqrtf[0][0]), r[i][1]/(2*sqrtf[0][1]));
153 SBasis cisi = shift(ci, i);
154 r -= multiply(shift((sqrtf*2 + cisi), i), SBasis(ci));
155 r.truncate(order+1);
156 sqrtf[i] = ci;
157 if(r.tailError(i) == 0) // if exact
158 break;
159 }
160 }else{
161 sqrtf = Linear(std::sqrt(fabs(f.at0())), std::sqrt(fabs(f.at1())));
162 }
164 double err = (f - multiply(sqrtf, sqrtf)).tailError(0);
165 if (err<tol){
166 return Piecewise<SBasis>(sqrtf);
167 }
169 Piecewise<SBasis> sqrtf0,sqrtf1;
170 sqrtf0 = sqrt_internal(compose(f,Linear(0.,.5)),tol,order);
171 sqrtf1 = sqrt_internal(compose(f,Linear(.5,1.)),tol,order);
172 sqrtf0.setDomain(Interval(0.,.5));
173 sqrtf1.setDomain(Interval(.5,1.));
174 sqrtf0.concat(sqrtf1);
175 return sqrtf0;
176 }
178 /** Compute the sqrt of a function.
179 \param f function
180 */
181 Piecewise<SBasis> sqrt(SBasis const &f, double tol, int order){
182 return sqrt(max(f,Linear(tol*tol)),tol,order);
183 }
185 /** Compute the sqrt of a function.
186 \param f function
187 */
188 Piecewise<SBasis> sqrt(Piecewise<SBasis> const &f, double tol, int order){
189 Piecewise<SBasis> result;
190 Piecewise<SBasis> zero = Piecewise<SBasis>(Linear(tol*tol));
191 zero.setDomain(f.domain());
192 Piecewise<SBasis> ff=max(f,zero);
194 for (unsigned i=0; i<ff.size(); i++){
195 Piecewise<SBasis> sqrtfi = sqrt_internal(ff.segs[i],tol,order);
196 sqrtfi.setDomain(Interval(ff.cuts[i],ff.cuts[i+1]));
197 result.concat(sqrtfi);
198 }
199 return result;
200 }
202 //-Yet another sin/cos--------------------------------------------------------------
204 /** Compute the sine of a function.
205 \param f function
206 \param tol maximum error
207 \param order maximum degree polynomial to use
208 */
209 Piecewise<SBasis> sin( SBasis const &f, double tol, int order){return(cos(-f+M_PI/2,tol,order));}
210 /** Compute the sine of a function.
211 \param f function
212 \param tol maximum error
213 \param order maximum degree polynomial to use
214 */
215 Piecewise<SBasis> sin(Piecewise<SBasis> const &f, double tol, int order){return(cos(-f+M_PI/2,tol,order));}
217 /** Compute the cosine of a function.
218 \param f function
219 \param tol maximum error
220 \param order maximum degree polynomial to use
221 */
222 Piecewise<SBasis> cos(Piecewise<SBasis> const &f, double tol, int order){
223 Piecewise<SBasis> result;
224 for (unsigned i=0; i<f.size(); i++){
225 Piecewise<SBasis> cosfi = cos(f.segs[i],tol,order);
226 cosfi.setDomain(Interval(f.cuts[i],f.cuts[i+1]));
227 result.concat(cosfi);
228 }
229 return result;
230 }
232 /** Compute the cosine of a function.
233 \param f function
234 \param tol maximum error
235 \param order maximum degree polynomial to use
236 */
237 Piecewise<SBasis> cos( SBasis const &f, double tol, int order){
238 double alpha = (f.at0()+f.at1())/2.;
239 SBasis x = f-alpha;
240 double d = x.tailError(0),err=1;
241 //estimate cos(x)-sum_0^order (-1)^k x^2k/2k! by the first neglicted term
242 for (int i=1; i<=2*order; i++) err*=d/i;
244 if (err<tol){
245 SBasis xk=Linear(1), c=Linear(1), s=Linear(0);
246 for (int k=1; k<=2*order; k+=2){
247 xk*=x/k;
248 //take also truncature errors into account...
249 err+=xk.tailError(order);
250 xk.truncate(order);
251 s+=xk;
252 xk*=-x/(k+1);
253 //take also truncature errors into account...
254 err+=xk.tailError(order);
255 xk.truncate(order);
256 c+=xk;
257 }
258 if (err<tol){
259 return Piecewise<SBasis>(std::cos(alpha)*c-std::sin(alpha)*s);
260 }
261 }
262 Piecewise<SBasis> c0,c1;
263 c0 = cos(compose(f,Linear(0.,.5)),tol,order);
264 c1 = cos(compose(f,Linear(.5,1.)),tol,order);
265 c0.setDomain(Interval(0.,.5));
266 c1.setDomain(Interval(.5,1.));
267 c0.concat(c1);
268 return c0;
269 }
271 //--1/x------------------------------------------------------------
272 //TODO: this implementation is just wrong. Remove or redo!
274 void truncateResult(Piecewise<SBasis> &f, int order){
275 if (order>=0){
276 for (unsigned k=0; k<f.segs.size(); k++){
277 f.segs[k].truncate(order);
278 }
279 }
280 }
282 Piecewise<SBasis> reciprocalOnDomain(Interval range, double tol){
283 Piecewise<SBasis> reciprocal_fn;
284 //TODO: deduce R from tol...
285 double R=2.;
286 SBasis reciprocal1_R=reciprocal(Linear(1,R),3);
287 double a=range.min(), b=range.max();
288 if (a*b<0){
289 b=std::max(fabs(a),fabs(b));
290 a=0;
291 }else if (b<0){
292 a=-range.max();
293 b=-range.min();
294 }
296 if (a<=tol){
297 reciprocal_fn.push_cut(0);
298 int i0=(int) floor(std::log(tol)/std::log(R));
299 a=pow(R,i0);
300 reciprocal_fn.push(Linear(1/a),a);
301 }else{
302 int i0=(int) floor(std::log(a)/std::log(R));
303 a=pow(R,i0);
304 reciprocal_fn.cuts.push_back(a);
305 }
307 while (a<b){
308 reciprocal_fn.push(reciprocal1_R/a,R*a);
309 a*=R;
310 }
311 if (range.min()<0 || range.max()<0){
312 Piecewise<SBasis>reciprocal_fn_neg;
313 //TODO: define reverse(pw<sb>);
314 reciprocal_fn_neg.cuts.push_back(-reciprocal_fn.cuts.back());
315 for (unsigned i=0; i<reciprocal_fn.size(); i++){
316 int idx=reciprocal_fn.segs.size()-1-i;
317 reciprocal_fn_neg.push_seg(-reverse(reciprocal_fn.segs.at(idx)));
318 reciprocal_fn_neg.push_cut(-reciprocal_fn.cuts.at(idx));
319 }
320 if (range.max()>0){
321 reciprocal_fn_neg.concat(reciprocal_fn);
322 }
323 reciprocal_fn=reciprocal_fn_neg;
324 }
326 return(reciprocal_fn);
327 }
329 Piecewise<SBasis> reciprocal(SBasis const &f, double tol, int order){
330 Piecewise<SBasis> reciprocal_fn=reciprocalOnDomain(bounds_fast(f), tol);
331 Piecewise<SBasis> result=compose(reciprocal_fn,f);
332 truncateResult(result,order);
333 return(result);
334 }
335 Piecewise<SBasis> reciprocal(Piecewise<SBasis> const &f, double tol, int order){
336 Piecewise<SBasis> reciprocal_fn=reciprocalOnDomain(bounds_fast(f), tol);
337 Piecewise<SBasis> result=compose(reciprocal_fn,f);
338 truncateResult(result,order);
339 return(result);
340 }
342 }
344 /*
345 Local Variables:
346 mode:c++
347 c-file-style:"stroustrup"
348 c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
349 indent-tabs-mode:nil
350 fill-column:99
351 End:
352 */
353 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :