sbasis-math.cpp
Go to the documentation of this file.
00001 /* 00002 * sbasis-math.cpp - some std functions to work with (pw)s-basis 00003 * 00004 * Authors: 00005 * Jean-Francois Barraud 00006 * 00007 * Copyright (C) 2006-2007 authors 00008 * 00009 * This library is free software; you can redistribute it and/or 00010 * modify it either under the terms of the GNU Lesser General Public 00011 * License version 2.1 as published by the Free Software Foundation 00012 * (the "LGPL") or, at your option, under the terms of the Mozilla 00013 * Public License Version 1.1 (the "MPL"). If you do not alter this 00014 * notice, a recipient may use your version of this file under either 00015 * the MPL or the LGPL. 00016 * 00017 * You should have received a copy of the LGPL along with this library 00018 * in the file COPYING-LGPL-2.1; if not, write to the Free Software 00019 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 00020 * You should have received a copy of the MPL along with this library 00021 * in the file COPYING-MPL-1.1 00022 * 00023 * The contents of this file are subject to the Mozilla Public License 00024 * Version 1.1 (the "License"); you may not use this file except in 00025 * compliance with the License. You may obtain a copy of the License at 00026 * http://www.mozilla.org/MPL/ 00027 * 00028 * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY 00029 * OF ANY KIND, either express or implied. See the LGPL or the MPL for 00030 * the specific language governing rights and limitations. 00031 */ 00032 00033 //this a first try to define sqrt, cos, sin, etc... 00034 //TODO: define a truncated compose(sb,sb, order) and extend it to pw<sb>. 00035 //TODO: in all these functions, compute 'order' according to 'tol'. 00036 00037 #include <2geom/sbasis-math.h> 00038 00039 #include <2geom/d2-sbasis.h> 00040 #include <stdio.h> 00041 #include <math.h> 00042 //#define ZERO 1e-3 00043 00044 00045 namespace Geom { 00046 00047 00048 //-|x|----------------------------------------------------------------------- 00052 Piecewise<SBasis> abs(SBasis const &f){ 00053 return abs(Piecewise<SBasis>(f)); 00054 } 00058 Piecewise<SBasis> abs(Piecewise<SBasis> const &f){ 00059 Piecewise<SBasis> absf=partition(f,roots(f)); 00060 for (unsigned i=0; i<absf.size(); i++){ 00061 if (absf.segs[i](.5)<0) absf.segs[i]*=-1; 00062 } 00063 return absf; 00064 } 00065 00066 //-max(x,y), min(x,y)-------------------------------------------------------- 00070 Piecewise<SBasis> max( SBasis const &f, SBasis const &g){ 00071 return max(Piecewise<SBasis>(f),Piecewise<SBasis>(g)); 00072 } 00076 Piecewise<SBasis> max(Piecewise<SBasis> const &f, SBasis const &g){ 00077 return max(f,Piecewise<SBasis>(g)); 00078 } 00082 Piecewise<SBasis> max( SBasis const &f, Piecewise<SBasis> const &g){ 00083 return max(Piecewise<SBasis>(f),g); 00084 } 00088 Piecewise<SBasis> max(Piecewise<SBasis> const &f, Piecewise<SBasis> const &g){ 00089 Piecewise<SBasis> max=partition(f,roots(f-g)); 00090 Piecewise<SBasis> gg =partition(g,max.cuts); 00091 max = partition(max,gg.cuts); 00092 for (unsigned i=0; i<max.size(); i++){ 00093 if (max.segs[i](.5)<gg.segs[i](.5)) max.segs[i]=gg.segs[i]; 00094 } 00095 return max; 00096 } 00097 00101 Piecewise<SBasis> 00102 min( SBasis const &f, SBasis const &g){ return -max(-f,-g); } 00106 Piecewise<SBasis> 00107 min(Piecewise<SBasis> const &f, SBasis const &g){ return -max(-f,-g); } 00111 Piecewise<SBasis> 00112 min( SBasis const &f, Piecewise<SBasis> const &g){ return -max(-f,-g); } 00116 Piecewise<SBasis> 00117 min(Piecewise<SBasis> const &f, Piecewise<SBasis> const &g){ return -max(-f,-g); } 00118 00119 00120 //-sign(x)--------------------------------------------------------------- 00124 Piecewise<SBasis> signSb(SBasis const &f){ 00125 return signSb(Piecewise<SBasis>(f)); 00126 } 00130 Piecewise<SBasis> signSb(Piecewise<SBasis> const &f){ 00131 Piecewise<SBasis> sign=partition(f,roots(f)); 00132 for (unsigned i=0; i<sign.size(); i++){ 00133 sign.segs[i] = (sign.segs[i](.5)<0)? Linear(-1.):Linear(1.); 00134 } 00135 return sign; 00136 } 00137 00138 //-Sqrt---------------------------------------------------------- 00139 static Piecewise<SBasis> sqrt_internal(SBasis const &f, 00140 double tol, 00141 int order){ 00142 SBasis sqrtf; 00143 if(f.isZero() || order == 0){ 00144 return Piecewise<SBasis>(sqrtf); 00145 } 00146 if (f.at0()<-tol*tol && f.at1()<-tol*tol){ 00147 return sqrt_internal(-f,tol,order); 00148 }else if (f.at0()>tol*tol && f.at1()>tol*tol){ 00149 sqrtf.resize(order+1, Linear(0,0)); 00150 sqrtf[0] = Linear(std::sqrt(f[0][0]), std::sqrt(f[0][1])); 00151 SBasis r = f - multiply(sqrtf, sqrtf); // remainder 00152 for(unsigned i = 1; int(i) <= order && i<r.size(); i++) { 00153 Linear ci(r[i][0]/(2*sqrtf[0][0]), r[i][1]/(2*sqrtf[0][1])); 00154 SBasis cisi = shift(ci, i); 00155 r -= multiply(shift((sqrtf*2 + cisi), i), SBasis(ci)); 00156 r.truncate(order+1); 00157 sqrtf[i] = ci; 00158 if(r.tailError(i) == 0) // if exact 00159 break; 00160 } 00161 }else{ 00162 sqrtf = Linear(std::sqrt(fabs(f.at0())), std::sqrt(fabs(f.at1()))); 00163 } 00164 00165 double err = (f - multiply(sqrtf, sqrtf)).tailError(0); 00166 if (err<tol){ 00167 return Piecewise<SBasis>(sqrtf); 00168 } 00169 00170 Piecewise<SBasis> sqrtf0,sqrtf1; 00171 sqrtf0 = sqrt_internal(compose(f,Linear(0.,.5)),tol,order); 00172 sqrtf1 = sqrt_internal(compose(f,Linear(.5,1.)),tol,order); 00173 sqrtf0.setDomain(Interval(0.,.5)); 00174 sqrtf1.setDomain(Interval(.5,1.)); 00175 sqrtf0.concat(sqrtf1); 00176 return sqrtf0; 00177 } 00178 00182 Piecewise<SBasis> sqrt(SBasis const &f, double tol, int order){ 00183 return sqrt(max(f,Linear(tol*tol)),tol,order); 00184 } 00185 00189 Piecewise<SBasis> sqrt(Piecewise<SBasis> const &f, double tol, int order){ 00190 Piecewise<SBasis> result; 00191 Piecewise<SBasis> zero = Piecewise<SBasis>(Linear(tol*tol)); 00192 zero.setDomain(f.domain()); 00193 Piecewise<SBasis> ff=max(f,zero); 00194 00195 for (unsigned i=0; i<ff.size(); i++){ 00196 Piecewise<SBasis> sqrtfi = sqrt_internal(ff.segs[i],tol,order); 00197 sqrtfi.setDomain(Interval(ff.cuts[i],ff.cuts[i+1])); 00198 result.concat(sqrtfi); 00199 } 00200 return result; 00201 } 00202 00203 //-Yet another sin/cos-------------------------------------------------------------- 00204 00210 Piecewise<SBasis> sin( SBasis const &f, double tol, int order){return(cos(-f+M_PI/2,tol,order));} 00216 Piecewise<SBasis> sin(Piecewise<SBasis> const &f, double tol, int order){return(cos(-f+M_PI/2,tol,order));} 00217 00223 Piecewise<SBasis> cos(Piecewise<SBasis> const &f, double tol, int order){ 00224 Piecewise<SBasis> result; 00225 for (unsigned i=0; i<f.size(); i++){ 00226 Piecewise<SBasis> cosfi = cos(f.segs[i],tol,order); 00227 cosfi.setDomain(Interval(f.cuts[i],f.cuts[i+1])); 00228 result.concat(cosfi); 00229 } 00230 return result; 00231 } 00232 00238 Piecewise<SBasis> cos( SBasis const &f, double tol, int order){ 00239 double alpha = (f.at0()+f.at1())/2.; 00240 SBasis x = f-alpha; 00241 double d = x.tailError(0),err=1; 00242 //estimate cos(x)-sum_0^order (-1)^k x^2k/2k! by the first neglicted term 00243 for (int i=1; i<=2*order; i++) err*=d/i; 00244 00245 if (err<tol){ 00246 SBasis xk=Linear(1), c=Linear(1), s=Linear(0); 00247 for (int k=1; k<=2*order; k+=2){ 00248 xk*=x/k; 00249 //take also truncature errors into account... 00250 err+=xk.tailError(order); 00251 xk.truncate(order); 00252 s+=xk; 00253 xk*=-x/(k+1); 00254 //take also truncature errors into account... 00255 err+=xk.tailError(order); 00256 xk.truncate(order); 00257 c+=xk; 00258 } 00259 if (err<tol){ 00260 return Piecewise<SBasis>(std::cos(alpha)*c-std::sin(alpha)*s); 00261 } 00262 } 00263 Piecewise<SBasis> c0,c1; 00264 c0 = cos(compose(f,Linear(0.,.5)),tol,order); 00265 c1 = cos(compose(f,Linear(.5,1.)),tol,order); 00266 c0.setDomain(Interval(0.,.5)); 00267 c1.setDomain(Interval(.5,1.)); 00268 c0.concat(c1); 00269 return c0; 00270 } 00271 00272 //--1/x------------------------------------------------------------ 00273 //TODO: this implementation is just wrong. Remove or redo! 00274 00275 void truncateResult(Piecewise<SBasis> &f, int order){ 00276 if (order>=0){ 00277 for (unsigned k=0; k<f.segs.size(); k++){ 00278 f.segs[k].truncate(order); 00279 } 00280 } 00281 } 00282 00283 Piecewise<SBasis> reciprocalOnDomain(Interval range, double tol){ 00284 Piecewise<SBasis> reciprocal_fn; 00285 //TODO: deduce R from tol... 00286 double R=2.; 00287 SBasis reciprocal1_R=reciprocal(Linear(1,R),3); 00288 double a=range.min(), b=range.max(); 00289 if (a*b<0){ 00290 b=std::max(fabs(a),fabs(b)); 00291 a=0; 00292 }else if (b<0){ 00293 a=-range.max(); 00294 b=-range.min(); 00295 } 00296 00297 if (a<=tol){ 00298 reciprocal_fn.push_cut(0); 00299 int i0=(int) floor(std::log(tol)/std::log(R)); 00300 a=pow(R,i0); 00301 reciprocal_fn.push(Linear(1/a),a); 00302 }else{ 00303 int i0=(int) floor(std::log(a)/std::log(R)); 00304 a=pow(R,i0); 00305 reciprocal_fn.cuts.push_back(a); 00306 } 00307 00308 while (a<b){ 00309 reciprocal_fn.push(reciprocal1_R/a,R*a); 00310 a*=R; 00311 } 00312 if (range.min()<0 || range.max()<0){ 00313 Piecewise<SBasis>reciprocal_fn_neg; 00314 //TODO: define reverse(pw<sb>); 00315 reciprocal_fn_neg.cuts.push_back(-reciprocal_fn.cuts.back()); 00316 for (unsigned i=0; i<reciprocal_fn.size(); i++){ 00317 int idx=reciprocal_fn.segs.size()-1-i; 00318 reciprocal_fn_neg.push_seg(-reverse(reciprocal_fn.segs.at(idx))); 00319 reciprocal_fn_neg.push_cut(-reciprocal_fn.cuts.at(idx)); 00320 } 00321 if (range.max()>0){ 00322 reciprocal_fn_neg.concat(reciprocal_fn); 00323 } 00324 reciprocal_fn=reciprocal_fn_neg; 00325 } 00326 00327 return(reciprocal_fn); 00328 } 00329 00330 Piecewise<SBasis> reciprocal(SBasis const &f, double tol, int order){ 00331 Piecewise<SBasis> reciprocal_fn=reciprocalOnDomain(*bounds_fast(f), tol); 00332 Piecewise<SBasis> result=compose(reciprocal_fn,f); 00333 truncateResult(result,order); 00334 return(result); 00335 } 00336 Piecewise<SBasis> reciprocal(Piecewise<SBasis> const &f, double tol, int order){ 00337 Piecewise<SBasis> reciprocal_fn=reciprocalOnDomain(*bounds_fast(f), tol); 00338 Piecewise<SBasis> result=compose(reciprocal_fn,f); 00339 truncateResult(result,order); 00340 return(result); 00341 } 00342 00350 Piecewise<SBasis> interpolate(std::vector<double> times, std::vector<double> values, unsigned smoothness){ 00351 assert ( values.size() == times.size() ); 00352 if ( values.size() == 0 ) return Piecewise<SBasis>(); 00353 if ( values.size() == 1 ) return Piecewise<SBasis>(values[0]);//what about time?? 00354 00355 SBasis sk = shift(Linear(1.),smoothness); 00356 SBasis bump_in = integral(sk); 00357 bump_in -= bump_in.at0(); 00358 bump_in /= bump_in.at1(); 00359 SBasis bump_out = reverse( bump_in ); 00360 00361 Piecewise<SBasis> result; 00362 result.cuts.push_back(times[0]); 00363 for (unsigned i = 0; i<values.size()-1; i++){ 00364 result.push(bump_out*values[i]+bump_in*values[i+1],times[i+1]); 00365 } 00366 return result; 00367 } 00368 00369 } 00370 00371 /* 00372 Local Variables: 00373 mode:c++ 00374 c-file-style:"stroustrup" 00375 c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +)) 00376 indent-tabs-mode:nil 00377 fill-column:99 00378 End: 00379 */ 00380 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :