summaryrefslogtreecommitdiffstats
path: root/scribus/plugins/tools/2geomtools/lib2geom/sbasis-math.cpp
diff options
context:
space:
mode:
authorcraig <craig@11d20701-8431-0410-a711-e3c959e3b870>2012-01-01 11:40:09 +0000
committercraig <craig@11d20701-8431-0410-a711-e3c959e3b870>2012-01-01 11:40:09 +0000
commit7ed83b6c6666eb8b6b104c211ae7e52907350372 (patch)
tree4430b556abac0ad660a0aacf1887d77f85d8be02 /scribus/plugins/tools/2geomtools/lib2geom/sbasis-math.cpp
downloadscribus-7ed83b6c6666eb8b6b104c211ae7e52907350372.tar.gz
scribus-7ed83b6c6666eb8b6b104c211ae7e52907350372.tar.xz
scribus-7ed83b6c6666eb8b6b104c211ae7e52907350372.zip
Branch 1.3.5 tree to 1.4.x tree, goodbye 1.3.x
git-svn-id: svn://scribus.net/branches/Version14x/Scribus@17163 11d20701-8431-0410-a711-e3c959e3b870
Diffstat (limited to 'scribus/plugins/tools/2geomtools/lib2geom/sbasis-math.cpp')
-rw-r--r--scribus/plugins/tools/2geomtools/lib2geom/sbasis-math.cpp291
1 files changed, 291 insertions, 0 deletions
diff --git a/scribus/plugins/tools/2geomtools/lib2geom/sbasis-math.cpp b/scribus/plugins/tools/2geomtools/lib2geom/sbasis-math.cpp
new file mode 100644
index 0000000..1166a29
--- /dev/null
+++ b/scribus/plugins/tools/2geomtools/lib2geom/sbasis-math.cpp
@@ -0,0 +1,291 @@
+/*
+ * 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<sb>.
+//TODO: in all these functions, compute 'order' according to 'tol'.
+
+#include "sbasis-math.h"
+//#define ZERO 1e-3
+
+#include <cstdio>
+#include <cmath>
+#include "angle.h"
+
+namespace Geom {
+
+//-|x|-----------------------------------------------------------------------
+Piecewise<SBasis> abs(SBasis const &f){
+ return abs(Piecewise<SBasis>(f));
+}
+Piecewise<SBasis> abs(Piecewise<SBasis> const &f){
+ Piecewise<SBasis> absf=partition(f,roots(f));
+ for (unsigned i=0; i<absf.size(); i++){
+ if (absf.segs[i](.5)<0) absf.segs[i]*=-1;
+ }
+ return absf;
+}
+
+//-max(x,y), min(x,y)--------------------------------------------------------
+Piecewise<SBasis> max( SBasis const &f, SBasis const &g){
+ return max(Piecewise<SBasis>(f),Piecewise<SBasis>(g));
+}
+Piecewise<SBasis> max(Piecewise<SBasis> const &f, SBasis const &g){
+ return max(f,Piecewise<SBasis>(g));
+}
+Piecewise<SBasis> max( SBasis const &f, Piecewise<SBasis> const &g){
+ return max(Piecewise<SBasis>(f),g);
+}
+Piecewise<SBasis> max(Piecewise<SBasis> const &f, Piecewise<SBasis> const &g){
+ Piecewise<SBasis> max=partition(f,roots(f-g));
+ Piecewise<SBasis> gg =partition(g,max.cuts);
+ max = partition(max,gg.cuts);
+ for (unsigned i=0; i<max.size(); i++){
+ if (max.segs[i](.5)<gg.segs[i](.5)) max.segs[i]=gg.segs[i];
+ }
+ return max;
+}
+
+Piecewise<SBasis>
+min( SBasis const &f, SBasis const &g){ return -max(-f,-g); }
+Piecewise<SBasis>
+min(Piecewise<SBasis> const &f, SBasis const &g){ return -max(-f,-g); }
+Piecewise<SBasis>
+min( SBasis const &f, Piecewise<SBasis> const &g){ return -max(-f,-g); }
+Piecewise<SBasis>
+min(Piecewise<SBasis> const &f, Piecewise<SBasis> const &g){ return -max(-f,-g); }
+
+
+//-sign(x)---------------------------------------------------------------
+Piecewise<SBasis> signSb(SBasis const &f){
+ return signSb(Piecewise<SBasis>(f));
+}
+Piecewise<SBasis> signSb(Piecewise<SBasis> const &f){
+ Piecewise<SBasis> sign=partition(f,roots(f));
+ for (unsigned i=0; i<sign.size(); i++){
+ sign.segs[i] = (sign.segs[i](.5)<0)? Linear(-1.):Linear(1.);
+ }
+ return sign;
+}
+
+//-Sqrt----------------------------------------------------------
+static Piecewise<SBasis> sqrt_internal(SBasis const &f,
+ double tol,
+ int order){
+ SBasis sqrtf;
+ if(f.isZero() || order == 0){
+ return Piecewise<SBasis>(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 && i<r.size(); ++i) {
+ Linear ci(r[i][0]/(2*sqrtf[0][0]), r[i][1]/(2*sqrtf[0][1]));
+ SBasis cisi = shift(ci, i);
+ r -= multiply(shift((sqrtf*2 + cisi), i), SBasis(ci));
+ r.truncate(order+1);
+ sqrtf[i] = ci;
+ if(r.tailError(i) == 0) // if exact
+ break;
+ }
+ }else{
+ sqrtf = Linear(std::sqrt(fabs(f.at0())), std::sqrt(fabs(f.at1())));
+ }
+
+ double err = (f - multiply(sqrtf, sqrtf)).tailError(0);
+ if (err<tol){
+ return Piecewise<SBasis>(sqrtf);
+ }
+
+ Piecewise<SBasis> 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<SBasis> sqrt(SBasis const &f, double tol, int order){
+ return sqrt(max(f,Linear(tol*tol)),tol,order);
+}
+
+Piecewise<SBasis> sqrt(Piecewise<SBasis> const &f, double tol, int order){
+ Piecewise<SBasis> result;
+ Piecewise<SBasis> zero = Piecewise<SBasis>(Linear(tol*tol));
+ zero.setDomain(f.domain());
+ Piecewise<SBasis> ff=max(f,zero);
+
+ for (unsigned i=0; i<ff.size(); i++){
+ Piecewise<SBasis> 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<SBasis> sin( SBasis const &f, double tol, int order){return cos(-f+M_PI_2,tol,order);}
+Piecewise<SBasis> sin(Piecewise<SBasis> const &f, double tol, int order){return cos(-f+M_PI_2,tol,order);}
+
+Piecewise<SBasis> cos(Piecewise<SBasis> const &f, double tol, int order){
+ Piecewise<SBasis> result;
+ for (unsigned i=0; i<f.size(); i++){
+ Piecewise<SBasis> cosfi = cos(f.segs[i],tol,order);
+ cosfi.setDomain(Interval(f.cuts[i],f.cuts[i+1]));
+ result.concat(cosfi);
+ }
+ return result;
+}
+
+Piecewise<SBasis> 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<tol){
+ SBasis xk=Linear(1), c=Linear(1), s=Linear(0);
+ for (int k=1; k<=2*order; k+=2){
+ xk*=x/k;
+ //take also truncature errors into account...
+ err+=xk.tailError(order);
+ xk.truncate(order);
+ s+=xk;
+ xk*=-x/(k+1);
+ //take also truncature errors into account...
+ err+=xk.tailError(order);
+ xk.truncate(order);
+ c+=xk;
+ }
+ if (err<tol){
+ return Piecewise<SBasis>(std::cos(alpha)*c-std::sin(alpha)*s);
+ }
+ }
+ Piecewise<SBasis> 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<SBasis> &f, int order){
+ if (order>=0){
+ for (unsigned k=0; k<f.segs.size(); k++){
+ f.segs[k].truncate(order);
+ }
+ }
+}
+
+Piecewise<SBasis> reciprocalOnDomain(Interval range, double tol){
+ Piecewise<SBasis> 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(std::log(tol)/std::log(R));
+ a=pow(R,i0);
+ reciprocal_fn.push(Linear(1/a),a);
+ }else{
+ int i0=(int) floor(std::log(a)/std::log(R));
+ a=pow(R,i0);
+ reciprocal_fn.cuts.push_back(a);
+ }
+
+ while (a<b){
+ reciprocal_fn.push(reciprocal1_R/a,R*a);
+ a*=R;
+ }
+ if (range.min()<0 || range.max()<0){
+ Piecewise<SBasis>reciprocal_fn_neg;
+ //TODO: define reverse(pw<sb>);
+ reciprocal_fn_neg.cuts.push_back(-reciprocal_fn.cuts.back());
+ for (unsigned i=0; i<reciprocal_fn.size(); i++){
+ int idx=reciprocal_fn.segs.size()-1-i;
+ reciprocal_fn_neg.push_seg(-reverse(reciprocal_fn.segs.at(idx)));
+ reciprocal_fn_neg.push_cut(-reciprocal_fn.cuts.at(idx));
+ }
+ if (range.max()>0){
+ reciprocal_fn_neg.concat(reciprocal_fn);
+ }
+ reciprocal_fn=reciprocal_fn_neg;
+ }
+
+ return(reciprocal_fn);
+}
+
+Piecewise<SBasis> reciprocal(SBasis const &f, double tol, int order){
+ Piecewise<SBasis> reciprocal_fn=reciprocalOnDomain(bounds_fast(f), tol);
+ Piecewise<SBasis> result=compose(reciprocal_fn,f);
+ truncateResult(result,order);
+ return(result);
+}
+Piecewise<SBasis> reciprocal(Piecewise<SBasis> const &f, double tol, int order){
+ Piecewise<SBasis> reciprocal_fn=reciprocalOnDomain(bounds_fast(f), tol);
+ Piecewise<SBasis> 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 :