From 55d43e4e27e0ba58a47fad70957dfa989aa173ad Mon Sep 17 00:00:00 2001 From: "Johan B. C. Engelen" Date: Tue, 14 Aug 2007 20:54:48 +0000 Subject: Commit LivePathEffect branch to trunk! (disabled extension/internal/bitmap/*.* in build.xml to fix compilation) (bzr r3472) --- src/2geom/sbasis-math.cpp | 289 ++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 289 insertions(+) create mode 100644 src/2geom/sbasis-math.cpp (limited to 'src/2geom/sbasis-math.cpp') diff --git a/src/2geom/sbasis-math.cpp b/src/2geom/sbasis-math.cpp new file mode 100644 index 000000000..d88c26832 --- /dev/null +++ b/src/2geom/sbasis-math.cpp @@ -0,0 +1,289 @@ +/* + * sbasis-math.cpp - some std functions to work with (pw)s-basis + * + * Authors: + * Jean-Francois Barraud + * + * Copyright (C) 2006-2007 authors + * + * This library is free software; you can redistribute it and/or + * modify it either under the terms of the GNU Lesser General Public + * License version 2.1 as published by the Free Software Foundation + * (the "LGPL") or, at your option, under the terms of the Mozilla + * Public License Version 1.1 (the "MPL"). If you do not alter this + * notice, a recipient may use your version of this file under either + * the MPL or the LGPL. + * + * You should have received a copy of the LGPL along with this library + * in the file COPYING-LGPL-2.1; if not, write to the Free Software + * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + * You should have received a copy of the MPL along with this library + * in the file COPYING-MPL-1.1 + * + * The contents of this file are subject to the Mozilla Public License + * Version 1.1 (the "License"); you may not use this file except in + * compliance with the License. You may obtain a copy of the License at + * http://www.mozilla.org/MPL/ + * + * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY + * OF ANY KIND, either express or implied. See the LGPL or the MPL for + * the specific language governing rights and limitations. + */ + +//this a first try to define sqrt, cos, sin, etc... +//TODO: define a truncated compose(sb,sb, order) and extend it to pw. +//TODO: in all these functions, compute 'order' according to 'tol'. + +#include "sbasis-math.h" +//#define ZERO 1e-3 + + +namespace Geom { + + +#include +#include + +//-|x|----------------------------------------------------------------------- +Piecewise abs(SBasis const &f){ + return abs(Piecewise(f)); +} +Piecewise abs(Piecewise const &f){ + Piecewise absf=partition(f,roots(f)); + for (unsigned i=0; i maxSb( SBasis const &f, SBasis const &g){ + return maxSb(Piecewise(f),Piecewise(g)); +} +Piecewise maxSb(Piecewise const &f, SBasis const &g){ + return maxSb(f,Piecewise(g)); +} +Piecewise maxSb( SBasis const &f, Piecewise const &g){ + return maxSb(Piecewise(f),g); +} +Piecewise maxSb(Piecewise const &f, Piecewise const &g){ + Piecewise maxSb=partition(f,roots(f-g)); + Piecewise gg =partition(g,maxSb.cuts); + maxSb = partition(maxSb,gg.cuts); + for (unsigned i=0; i +minSb( SBasis const &f, SBasis const &g){ return -maxSb(-f,-g); } +Piecewise +minSb(Piecewise const &f, SBasis const &g){ return -maxSb(-f,-g); } +Piecewise +minSb( SBasis const &f, Piecewise const &g){ return -maxSb(-f,-g); } +Piecewise +minSb(Piecewise const &f, Piecewise const &g){ return -maxSb(-f,-g); } + + +//-sign(x)--------------------------------------------------------------- +Piecewise signSb(SBasis const &f){ + return signSb(Piecewise(f)); +} +Piecewise signSb(Piecewise const &f){ + Piecewise sign=partition(f,roots(f)); + for (unsigned i=0; i sqrt_internal(SBasis const &f, + double tol, + int order){ + SBasis sqrtf; + if(f.isZero() || order == 0){ + return Piecewise(sqrtf); + } + if (f.at0()<-tol*tol && f.at1()<-tol*tol){ + return sqrt_internal(-f,tol,order); + }else if (f.at0()>tol*tol && f.at1()>tol*tol){ + sqrtf.resize(order+1, Linear(0,0)); + sqrtf[0] = Linear(std::sqrt(f[0][0]), std::sqrt(f[0][1])); + SBasis r = f - multiply(sqrtf, sqrtf); // remainder + for(unsigned i = 1; int(i) <= order and i(sqrtf); + } + + Piecewise sqrtf0,sqrtf1; + sqrtf0 = sqrt_internal(compose(f,Linear(0.,.5)),tol,order); + sqrtf1 = sqrt_internal(compose(f,Linear(.5,1.)),tol,order); + sqrtf0.setDomain(Interval(0.,.5)); + sqrtf1.setDomain(Interval(.5,1.)); + sqrtf0.concat(sqrtf1); + return sqrtf0; +} + +Piecewise sqrt(SBasis const &f, double tol, int order){ + return sqrt(maxSb(f,Linear(tol*tol)),tol,order); +} + +Piecewise sqrt(Piecewise const &f, double tol, int order){ + Piecewise result; + Piecewise ff=maxSb(f,Linear(tol*tol)); + + for (unsigned i=0; i sqrtfi = sqrt_internal(ff.segs[i],tol,order); + sqrtfi.setDomain(Interval(ff.cuts[i],ff.cuts[i+1])); + result.concat(sqrtfi); + } + return result; +} + +//-Yet another sin/cos-------------------------------------------------------------- + +Piecewise sin( SBasis const &f, double tol, int order){return(cos(-f+M_PI/2,tol,order));} +Piecewise sin(Piecewise const &f, double tol, int order){return(cos(-f+M_PI/2,tol,order));} + +Piecewise cos(Piecewise const &f, double tol, int order){ + Piecewise result; + for (unsigned i=0; i cosfi = cos(f.segs[i],tol,order); + cosfi.setDomain(Interval(f.cuts[i],f.cuts[i+1])); + result.concat(cosfi); + } + return result; +} + +Piecewise cos( SBasis const &f, double tol, int order){ + double alpha = (f.at0()+f.at1())/2.; + SBasis x = f-alpha; + double d = x.tailError(0),err=1; + //estimate cos(x)-sum_0^order (-1)^k x^2k/2k! by the first neglicted term + for (int i=1; i<=2*order; i++) err*=d/i; + + if (err(std::cos(alpha)*c-std::sin(alpha)*s); + } + } + Piecewise c0,c1; + c0 = cos(compose(f,Linear(0.,.5)),tol,order); + c1 = cos(compose(f,Linear(.5,1.)),tol,order); + c0.setDomain(Interval(0.,.5)); + c1.setDomain(Interval(.5,1.)); + c0.concat(c1); + return c0; +} + +//--1/x------------------------------------------------------------ +//TODO: this implementation is just wrong. Remove or redo! + +void truncateResult(Piecewise &f, int order){ + if (order>=0){ + for (unsigned k=0; k reciprocalOnDomain(Interval range, double tol){ + Piecewise reciprocal_fn; + //TODO: deduce R from tol... + double R=2.; + SBasis reciprocal1_R=reciprocal(Linear(1,R),3); + double a=range.min(), b=range.max(); + if (a*b<0){ + b=std::max(fabs(a),fabs(b)); + a=0; + }else if (b<0){ + a=-range.max(); + b=-range.min(); + } + + if (a<=tol){ + reciprocal_fn.push_cut(0); + int i0=(int) floor(log(tol)/log(R)); + a=pow(R,i0); + reciprocal_fn.push(Linear(1/a),a); + }else{ + int i0=(int) floor(log(a)/log(R)); + a=pow(R,i0); + reciprocal_fn.cuts.push_back(a); + } + + while (areciprocal_fn_neg; + //TODO: define reverse(pw); + reciprocal_fn_neg.cuts.push_back(-reciprocal_fn.cuts.back()); + for (unsigned i=0; i0){ + reciprocal_fn_neg.concat(reciprocal_fn); + } + reciprocal_fn=reciprocal_fn_neg; + } + + return(reciprocal_fn); +} + +Piecewise reciprocal(SBasis const &f, double tol, int order){ + Piecewise reciprocal_fn=reciprocalOnDomain(bounds_fast(f), tol); + Piecewise result=compose(reciprocal_fn,f); + truncateResult(result,order); + return(result); +} +Piecewise reciprocal(Piecewise const &f, double tol, int order){ + Piecewise reciprocal_fn=reciprocalOnDomain(bounds_fast(f), tol); + Piecewise result=compose(reciprocal_fn,f); + truncateResult(result,order); + return(result); +} + +} + +/* + Local Variables: + mode:c++ + c-file-style:"stroustrup" + c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +)) + indent-tabs-mode:nil + fill-column:99 + End: +*/ +// vim: filetype = cpp:expandtab:shiftwidth = 4:tabstop = 8:softtabstop = 4:encoding = utf-8:textwidth = 99 : -- cgit v1.2.3