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
  39 #include <2geom/d2-sbasis.h>
  40 #include <stdio.h>
  41 #include <math.h>
  42 //#define ZERO 1e-3
  43
  44
  45 namespace Geom {
  46
  47
  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 }
  65
  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 }
  97
  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); }
 118
 119
 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 }
 137
 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     }
 164
 165     double err = (f - multiply(sqrtf, sqrtf)).tailError(0);
 166     if (err<tol){
 167         return Piecewise<SBasis>(sqrtf);
 168     }
 169
 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 }
 178
 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 }
 185
 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);
 194
 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 }
 202
 203 //-Yet another sin/cos--------------------------------------------------------------
 204
 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));}
 217
 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 }
 232
 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;
 244
 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 }
 271
 272 //--1/x------------------------------------------------------------
 273 //TODO: this implementation is just wrong. Remove or redo!
 274
 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 }
 282
 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     }
 296
 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     }
 307
 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     }
 326
 327     return(reciprocal_fn);
 328 }
 329
 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 }
 342
 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??
 354
 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 );
 360
 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 }
 368
 369 }
 370
 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:fileencoding=utf-8:textwidth=99 :