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 :
Inkscape

sbasis-math.cpp
