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