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> ff=max(f,Linear(tol*tol));
 149
 150     for (unsigned i=0; i<ff.size(); i++){
 151         Piecewise<SBasis> sqrtfi = sqrt_internal(ff.segs[i],tol,order);
 152         sqrtfi.setDomain(Interval(ff.cuts[i],ff.cuts[i+1]));
 153         result.concat(sqrtfi);
 154     }
 155     return result;
 156 }
 157
 158 //-Yet another sin/cos--------------------------------------------------------------
 159
 160 Piecewise<SBasis> sin(          SBasis  const &f, double tol, int order){return(cos(-f+M_PI/2,tol,order));}
 161 Piecewise<SBasis> sin(Piecewise<SBasis> const &f, double tol, int order){return(cos(-f+M_PI/2,tol,order));}
 162
 163 Piecewise<SBasis> cos(Piecewise<SBasis> const &f, double tol, int order){
 164     Piecewise<SBasis> result;
 165     for (unsigned i=0; i<f.size(); i++){
 166         Piecewise<SBasis> cosfi = cos(f.segs[i],tol,order);
 167         cosfi.setDomain(Interval(f.cuts[i],f.cuts[i+1]));
 168         result.concat(cosfi);
 169     }
 170     return result;
 171 }
 172
 173 Piecewise<SBasis> cos(          SBasis  const &f, double tol, int order){
 174     double alpha = (f.at0()+f.at1())/2.;
 175     SBasis x = f-alpha;
 176     double d = x.tailError(0),err=1;
 177     //estimate cos(x)-sum_0^order (-1)^k x^2k/2k! by the first neglicted term
 178     for (int i=1; i<=2*order; i++) err*=d/i;
 179
 180     if (err<tol){
 181         SBasis xk=Linear(1), c=Linear(1), s=Linear(0);
 182         for (int k=1; k<=2*order; k+=2){
 183             xk*=x/k;
 184             //take also truncature errors into account...
 185             err+=xk.tailError(order);
 186             xk.truncate(order);
 187             s+=xk;
 188             xk*=-x/(k+1);
 189             //take also truncature errors into account...
 190             err+=xk.tailError(order);
 191             xk.truncate(order);
 192             c+=xk;
 193         }
 194         if (err<tol){
 195             return Piecewise<SBasis>(std::cos(alpha)*c-std::sin(alpha)*s);
 196         }
 197     }
 198     Piecewise<SBasis> c0,c1;
 199     c0 = cos(compose(f,Linear(0.,.5)),tol,order);
 200     c1 = cos(compose(f,Linear(.5,1.)),tol,order);
 201     c0.setDomain(Interval(0.,.5));
 202     c1.setDomain(Interval(.5,1.));
 203     c0.concat(c1);
 204     return c0;
 205 }
 206
 207 //--1/x------------------------------------------------------------
 208 //TODO: this implementation is just wrong. Remove or redo!
 209
 210 void truncateResult(Piecewise<SBasis> &f, int order){
 211     if (order>=0){
 212         for (unsigned k=0; k<f.segs.size(); k++){
 213             f.segs[k].truncate(order);
 214         }
 215     }
 216 }
 217
 218 Piecewise<SBasis> reciprocalOnDomain(Interval range, double tol){
 219     Piecewise<SBasis> reciprocal_fn;
 220     //TODO: deduce R from tol...
 221     double R=2.;
 222     SBasis reciprocal1_R=reciprocal(Linear(1,R),3);
 223     double a=range.min(), b=range.max();
 224     if (a*b<0){
 225         b=std::max(fabs(a),fabs(b));
 226         a=0;
 227     }else if (b<0){
 228         a=-range.max();
 229         b=-range.min();
 230     }
 231
 232     if (a<=tol){
 233         reciprocal_fn.push_cut(0);
 234         int i0=(int) floor(std::log(tol)/std::log(R));
 235         a=pow(R,i0);
 236         reciprocal_fn.push(Linear(1/a),a);
 237     }else{
 238         int i0=(int) floor(std::log(a)/std::log(R));
 239         a=pow(R,i0);
 240         reciprocal_fn.cuts.push_back(a);
 241     }
 242
 243     while (a<b){
 244         reciprocal_fn.push(reciprocal1_R/a,R*a);
 245         a*=R;
 246     }
 247     if (range.min()<0 || range.max()<0){
 248         Piecewise<SBasis>reciprocal_fn_neg;
 249         //TODO: define reverse(pw<sb>);
 250         reciprocal_fn_neg.cuts.push_back(-reciprocal_fn.cuts.back());
 251         for (unsigned i=0; i<reciprocal_fn.size(); i++){
 252             int idx=reciprocal_fn.segs.size()-1-i;
 253             reciprocal_fn_neg.push_seg(-reverse(reciprocal_fn.segs.at(idx)));
 254             reciprocal_fn_neg.push_cut(-reciprocal_fn.cuts.at(idx));
 255         }
 256         if (range.max()>0){
 257             reciprocal_fn_neg.concat(reciprocal_fn);
 258         }
 259         reciprocal_fn=reciprocal_fn_neg;
 260     }
 261
 262     return(reciprocal_fn);
 263 }
 264
 265 Piecewise<SBasis> reciprocal(SBasis const &f, double tol, int order){
 266     Piecewise<SBasis> reciprocal_fn=reciprocalOnDomain(bounds_fast(f), tol);
 267     Piecewise<SBasis> result=compose(reciprocal_fn,f);
 268     truncateResult(result,order);
 269     return(result);
 270 }
 271 Piecewise<SBasis> reciprocal(Piecewise<SBasis> const &f, double tol, int order){
 272     Piecewise<SBasis> reciprocal_fn=reciprocalOnDomain(bounds_fast(f), tol);
 273     Piecewise<SBasis> result=compose(reciprocal_fn,f);
 274     truncateResult(result,order);
 275     return(result);
 276 }
 277
 278 }
 279
 280 /*
 281   Local Variables:
 282   mode:c++
 283   c-file-style:"stroustrup"
 284   c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
 285   indent-tabs-mode:nil
 286   fill-column:99
 287   End:
 288 */
 289 // vim: filetype = cpp:expandtab:shiftwidth = 4:tabstop = 8:softtabstop = 4:encoding = utf-8:textwidth = 99 :