src/2geom/sbasis-math.cpp

   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  */
  32
  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'.
  36
  37 #include <2geom/sbasis-math.h>
  38 //#define ZERO 1e-3
  39
  40 #include <stdio.h>
  41 #include <math.h>
  42
  43 namespace Geom {
  44
  45 #include <2geom/d2-sbasis.h>
  46
  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 }
  64
  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 }
  96
  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); }
 117
 118
 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 }
 136
 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     }
 163
 164     double err = (f - multiply(sqrtf, sqrtf)).tailError(0);
 165     if (err<tol){
 166         return Piecewise<SBasis>(sqrtf);
 167     }
 168
 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 }
 177
 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 }
 184
 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);
 193
 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 }
 201
 202 //-Yet another sin/cos--------------------------------------------------------------
 203
 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));}
 216
 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 }
 231
 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;
 243
 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 }
 270
 271 //--1/x------------------------------------------------------------
 272 //TODO: this implementation is just wrong. Remove or redo!
 273
 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 }
 281
 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     }
 295
 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     }
 306
 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     }
 325
 326     return(reciprocal_fn);
 327 }
 328
 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 }
 341
 342 /**
 343  * \brief Retruns a Piecewise SBasis with prescribed values at prescribed times.
 344  *
 345  * \param times: vector of times at which the values are given. Should be sorted in increasing order.
 346  * \param values: vector of prescribed values. Should have the same size as times and be sorted accordingly.
 347  * \param smoothness: (defaults to 1) regularity class of the result: 0=piecewise linear, 1=continuous derivative, etc...
 348  */
 349 Piecewise<SBasis> interpolate(std::vector<double> times, std::vector<double> values, unsigned smoothness){
 350     assert ( values.size() == times.size() );
 351     if ( values.size() == 0 ) return Piecewise<SBasis>();
 352     if ( values.size() == 1 ) return Piecewise<SBasis>(values[0]);//what about time??
 353
 354     SBasis sk = shift(Linear(1.),smoothness);
 355     SBasis bump_in = integral(sk);
 356     bump_in -= bump_in.at0();
 357     bump_in /= bump_in.at1();
 358     SBasis bump_out = reverse( bump_in );
 359
 360     Piecewise<SBasis> result;
 361     result.cuts.push_back(times[0]);
 362     for (unsigned i = 0; i<values.size()-1; i++){
 363         result.push(bump_out*values[i]+bump_in*values[i+1],times[i+1]);
 364     }
 365     return result;
 366 }
 367
 368 }
 369
 370 /*
 371   Local Variables:
 372   mode:c++
 373   c-file-style:"stroustrup"
 374   c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
 375   indent-tabs-mode:nil
 376   fill-column:99
 377   End:
 378 */
 379 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :