From 61f3ab8f23f4c924d455757bf3e65f8487521b5a Mon Sep 17 00:00:00 2001 From: Daniel Baumann Date: Sat, 13 Apr 2024 13:57:42 +0200 Subject: Adding upstream version 1.3. Signed-off-by: Daniel Baumann --- include/2geom/orphan-code/sbasis-of.h | 638 ++++++++++++++++++++++++++++++++++ 1 file changed, 638 insertions(+) create mode 100644 include/2geom/orphan-code/sbasis-of.h (limited to 'include/2geom/orphan-code/sbasis-of.h') diff --git a/include/2geom/orphan-code/sbasis-of.h b/include/2geom/orphan-code/sbasis-of.h new file mode 100644 index 0000000..e5b76d6 --- /dev/null +++ b/include/2geom/orphan-code/sbasis-of.h @@ -0,0 +1,638 @@ +/** + * \file + * \brief Defines S-power basis function class + * with coefficient in arbitrary ring + * + * Authors: + * Nathan Hurst + * Michael Sloan + * + * 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. + */ + +#ifndef SEEN_SBASIS_OF_H +#define SEEN_SBASIS_OF_H +#include +#include +#include + +#include <2geom/interval.h> +#include <2geom/utils.h> +#include <2geom/exception.h> + +#include <2geom/orphan-code/linear-of.h> + +namespace Geom { + +template +class SBasisOf; + +#ifdef USE_SBASIS_OF +typedef Geom::SBasisOf SBasis; +#endif + +/*** An empty SBasisOf is identically 0. */ +template +class SBasisOf : public std::vector >{ +public: + SBasisOf() {} + explicit SBasisOf(T a) { + this->push_back(LinearOf(a,a)); + } + SBasisOf(SBasisOf const & a) : + std::vector >(a) + {} + SBasisOf(LinearOf const & bo) { + this->push_back(bo); + } + SBasisOf(LinearOf* bo) { + this->push_back(*bo); + } + //static unsigned input_dim(){return T::input_dim()+1;} + + //IMPL: FragmentConcept + typedef T output_type; + inline bool isZero() const { + if(this->empty()) return true; + for(unsigned i = 0; i < this->size(); i++) { + if(!(*this)[i].isZero()) return false; + } + return true; + } + inline bool isConstant() const { + if (this->empty()) return true; + for (unsigned i = 0; i < this->size(); i++) { + if(!(*this)[i].isConstant()) return false; + } + return true; + } + + //TODO: code this... + bool isFinite() const; + + inline T at0() const { + if(this->empty()) return T(0); else return (*this)[0][0]; + } + inline T at1() const{ + if(this->empty()) return T(0); else return (*this)[0][1]; + } + + T valueAt(double t) const { + double s = t*(1-t); + T p0 = T(0.), p1 = T(0.); + for(unsigned k = this->size(); k > 0; k--) { + const LinearOf &lin = (*this)[k-1]; + p0 = p0*s + lin[0]; + p1 = p1*s + lin[1]; + } + return p0*(1-t) + p1*t; + } + + T operator()(double t) const { + return valueAt(t); + } + + /** + * The size of the returned vector equals n+1. + */ + std::vector valueAndDerivatives(double t, unsigned n) const{ + std::vector ret(n+1); + ret[0] = valueAt(t); + SBasisOf tmp = *this; + for(unsigned i = 0; i < n; i++) { + tmp.derive(); + ret[i+1] = tmp.valueAt(t); + } + return ret; + } + + //The following lines only makes sens if T=double! + SBasisOf toSBasis() const { return SBasisOf(*this); } + double tailError(unsigned tail) const{ + Interval bs = *bounds_fast(*this, tail); + return std::max(fabs(bs.min()),fabs(bs.max())); + } + +// compute f(g) + SBasisOf operator()(SBasisOf const & g) const; + + LinearOf operator[](unsigned i) const { + assert(i < this->size()); + return std::vector >::operator[](i); + } + +//MUTATOR PRISON + LinearOf& operator[](unsigned i) { return this->at(i); } + + //remove extra zeros + void normalize() { + while(!this->empty() && this->back().isZero()) + this->pop_back(); + } + + void truncate(unsigned k) { if(k < this->size()) this->resize(k); } +private: + void derive(); // in place version + unsigned dim; +}; + +//template<> +//inline unsigned SBasisOf::input_dim() { return 1; } + +//-------------------------------------------------------------------------- +#ifdef USE_SBASIS_OF + +//implemented in sbasis-roots.cpp +OptInterval bounds_exact(SBasis const &a); +OptInterval bounds_fast(SBasis const &a, int order = 0); +OptInterval bounds_local(SBasis const &a, const OptInterval &t, int order = 0); + +std::vector roots(SBasis const & s); +std::vector > multi_roots(SBasis const &f, + std::vector const &levels, + double htol=1e-7, + double vtol=1e-7, + double a=0, + double b=1); +#endif +//-------------------------------------------------------------------------- + + +//TODO: figure out how to stick this in linear, while not adding an sbasis dep +template +inline SBasisOf LinearOf::toSBasis() const { return SBasisOf(*this); } + +template +inline SBasisOf reverse(SBasisOf const &a) { + SBasisOf result; + result.reserve(a.size()); + for(unsigned k = 0; k < a.size(); k++) + result.push_back(reverse(a[k])); + return result; +} + +//IMPL: ScalableConcept +template +inline SBasisOf operator-(const SBasisOf& p) { + if(p.isZero()) return SBasisOf(); + SBasisOf result; + result.reserve(p.size()); + + for(unsigned i = 0; i < p.size(); i++) { + result.push_back(-p[i]); + } + return result; +} + +template +SBasisOf operator*(SBasisOf const &a, double k){ + SBasisOf c; + //TODO: what does this mean for vectors of vectors?? + //c.reserve(a.size()); + for(unsigned i = 0; i < a.size(); i++) + c.push_back(a[i] * k); + return c; +} + +template +inline SBasisOf operator*(double k, SBasisOf const &a) { return a*k; } +template +inline SBasisOf operator/(SBasisOf const &a, double k) { return a*(1./k); } +template +SBasisOf& operator*=(SBasisOf& a, double b){ + if (a.isZero()) return a; + if (b == 0) + a.clear(); + else + for(unsigned i = 0; i < a.size(); i++) + a[i] *= b; + return a; +} + +template +inline SBasisOf& operator/=(SBasisOf& a, double b) { return (a*=(1./b)); } + +/* +//We can also multiply by element of ring coeff T: +template +SBasisOf operator*(SBasisOf const &a, T k){ + SBasisOf c; + //TODO: what does this mean for vectors of vectors?? + //c.reserve(a.size()); + for(unsigned i = 0; i < a.size(); i++) + c.push_back(a[i] * k); + return c; +} + +template +inline SBasisOf operator*(T k, SBasisOf const &a) { return a*k; } +template +SBasisOf& operator*=(SBasisOf& a, T b){ + if (a.isZero()) return a; + if (b == 0) + a.clear(); + else + for(unsigned i = 0; i < a.size(); i++) + a[i] *= b; + return a; +} +*/ + +//IMPL: AddableConcept +template +inline SBasisOf operator+(const SBasisOf& a, const SBasisOf& b){ + SBasisOf result; + const unsigned out_size = std::max(a.size(), b.size()); + const unsigned min_size = std::min(a.size(), b.size()); + //TODO: what does this mean for vector; + //result.reserve(out_size); + + for(unsigned i = 0; i < min_size; i++) { + result.push_back(a[i] + b[i]); + } + for(unsigned i = min_size; i < a.size(); i++) + result.push_back(a[i]); + for(unsigned i = min_size; i < b.size(); i++) + result.push_back(b[i]); + + assert(result.size() == out_size); + return result; +} + +template +SBasisOf operator-(const SBasisOf& a, const SBasisOf& b){ + SBasisOf result; + const unsigned out_size = std::max(a.size(), b.size()); + const unsigned min_size = std::min(a.size(), b.size()); + //TODO: what does this mean for vector; + //result.reserve(out_size); + + for(unsigned i = 0; i < min_size; i++) { + result.push_back(a[i] - b[i]); + } + for(unsigned i = min_size; i < a.size(); i++) + result.push_back(a[i]); + for(unsigned i = min_size; i < b.size(); i++) + result.push_back(-b[i]); + + assert(result.size() == out_size); + return result; +} + +template +SBasisOf& operator+=(SBasisOf& a, const SBasisOf& b){ + const unsigned out_size = std::max(a.size(), b.size()); + const unsigned min_size = std::min(a.size(), b.size()); + //TODO: what does this mean for vectors of vectors + //a.reserve(out_size); + for(unsigned i = 0; i < min_size; i++) + a[i] += b[i]; + for(unsigned i = min_size; i < b.size(); i++) + a.push_back(b[i]); + + assert(a.size() == out_size); + return a; +} + +template +SBasisOf& operator-=(SBasisOf& a, const SBasisOf& b){ + const unsigned out_size = std::max(a.size(), b.size()); + const unsigned min_size = std::min(a.size(), b.size()); + //TODO: what does this mean for vectors of vectors + //a.reserve(out_size); + for(unsigned i = 0; i < min_size; i++) + a[i] -= b[i]; + for(unsigned i = min_size; i < b.size(); i++) + a.push_back(-b[i]); + + assert(a.size() == out_size); + return a; +} + +//TODO: remove? +template +inline SBasisOf operator+(const SBasisOf & a, LinearOf const & b) { + if(b.isZero()) return a; + if(a.isZero()) return b; + SBasisOf result(a); + result[0] += b; + return result; +} +template +inline SBasisOf operator-(const SBasisOf & a, LinearOf const & b) { + if(b.isZero()) return a; + SBasisOf result(a); + result[0] -= b; + return result; +} +template +inline SBasisOf& operator+=(SBasisOf& a, const LinearOf& b) { + if(a.isZero()) + a.push_back(b); + else + a[0] += b; + return a; +} +template +inline SBasisOf& operator-=(SBasisOf& a, const LinearOf& b) { + if(a.isZero()) + a.push_back(-b); + else + a[0] -= b; + return a; +} + +//IMPL: OffsetableConcept +/* +template +inline SBasisOf operator+(const SBasisOf & a, double b) { + if(a.isZero()) return LinearOf(b, b); + SBasisOf result(a); + result[0] += b; + return result; +} +template +inline SBasisOf operator-(const SBasisOf & a, double b) { + if(a.isZero()) return LinearOf(-b, -b); + SBasisOf result(a); + result[0] -= b; + return result; +} +template +inline SBasisOf& operator+=(SBasisOf& a, double b) { + if(a.isZero()) + a.push_back(LinearOf(b,b)); + else + a[0] += b; + return a; +} +template +inline SBasisOf& operator-=(SBasisOf& a, double b) { + if(a.isZero()) + a.push_back(LinearOf(-b,-b)); + else + a[0] -= b; + return a; +} +*/ +//We can also offset by elements of coeff ring T +template +inline SBasisOf operator+(const SBasisOf & a, T b) { + if(a.isZero()) return LinearOf(b, b); + SBasisOf result(a); + result[0] += b; + return result; +} +template +inline SBasisOf operator-(const SBasisOf & a, T b) { + if(a.isZero()) return LinearOf(-b, -b); + SBasisOf result(a); + result[0] -= b; + return result; +} +template +inline SBasisOf& operator+=(SBasisOf& a, T b) { + if(a.isZero()) + a.push_back(LinearOf(b,b)); + else + a[0] += b; + return a; +} +template +inline SBasisOf& operator-=(SBasisOf& a, T b) { + if(a.isZero()) + a.push_back(LinearOf(-b,-b)); + else + a[0] -= b; + return a; +} + + +template +SBasisOf shift(SBasisOf const &a, int sh){ + SBasisOf c = a; + if(sh > 0) { + c.insert(c.begin(), sh, LinearOf(0,0)); + } else { + //TODO: truncate + } + return c; +} + +template +SBasisOf shift(LinearOf const &a, int sh) { + SBasisOf c; + if(sh > 0) { + c.insert(c.begin(), sh, LinearOf(0,0)); + c.push_back(a); + } + return c; +} + +template +inline SBasisOf truncate(SBasisOf const &a, unsigned terms) { + SBasisOf c; + c.insert(c.begin(), a.begin(), a.begin() + std::min(terms, (unsigned)a.size())); + return c; +} + +template +SBasisOf multiply_add(SBasisOf const &a, SBasisOf const &b, SBasisOf c) { + if(a.isZero() || b.isZero()) + return c; + c.resize(a.size() + b.size(), LinearOf(T(0.),T(0.))); + for(unsigned j = 0; j < b.size(); j++) { + for(unsigned i = j; i < a.size()+j; i++) { + T tri = (b[j][1]-b[j][0])*(a[i-j][1]-a[i-j][0]); + c[i+1/*shift*/] += LinearOf(-tri); + } + } + for(unsigned j = 0; j < b.size(); j++) { + for(unsigned i = j; i < a.size()+j; i++) { + for(unsigned dim = 0; dim < 2; dim++) + c[i][dim] += b[j][dim]*a[i-j][dim]; + } + } + c.normalize(); + //assert(!(0 == c.back()[0] && 0 == c.back()[1])); + return c; +} + +template +SBasisOf multiply(SBasisOf const &a, SBasisOf const &b) { + SBasisOf c; + if(a.isZero() || b.isZero()) + return c; + return multiply_add(a, b, c); +} + +template +SBasisOf integral(SBasisOf const &c){ + SBasisOf a; + T aTri = T(0.); + for(int k = c.size()-1; k >= 0; k--) { + aTri = (HatOf(c[k]).d + (k+1)*aTri/2)/(2*k+1); + a[k][0] -= aTri/2; + a[k][1] += aTri/2; + } + a.normalize(); + return a; +} + +template +SBasisOf derivative(SBasisOf const &a){ + SBasisOf c; + c.resize(a.size(), LinearOf()); + if(a.isZero()) + return c; + + for(unsigned k = 0; k < a.size()-1; k++) { + T d = (2*k+1)*(a[k][1] - a[k][0]); + + c[k][0] = d + (k+1)*a[k+1][0]; + c[k][1] = d - (k+1)*a[k+1][1]; + } + int k = a.size()-1; + T d = (2*k+1)*(a[k][1] - a[k][0]); + //TODO: do a real test to know if d==0! + if(d == T(0.0)) + c.pop_back(); + else { + c[k][0] = d; + c[k][1] = d; + } + + return c; +} + +template +void SBasisOf::derive() { // in place version + if(isZero()) return; + for(unsigned k = 0; k < this->size()-1; k++) { + T d = (2*k+1)*((*this)[k][1] - (*this)[k][0]); + + (*this)[k][0] = d + (k+1)*(*this)[k+1][0]; + (*this)[k][1] = d - (k+1)*(*this)[k+1][1]; + } + int k = this->size()-1; + T d = (2*k+1)*((*this)[k][1] - (*this)[k][0]); + if(d == 0)//TODO: give this a meaning for general coeff ring. + this->pop_back(); + else { + (*this)[k][0] = d; + (*this)[k][1] = d; + } +} + + +template +inline SBasisOf operator*(SBasisOf const & a, SBasisOf const & b) { + return multiply(a, b); +} + +template +inline SBasisOf& operator*=(SBasisOf& a, SBasisOf const & b) { + a = multiply(a, b); + return a; +} + +// a(b(t)) +//TODO: compose all compatibles types! +template +SBasisOf compose(SBasisOf const &a, SBasisOf const &b){ + SBasisOf s = multiply((SBasisOf(LinearOf(1,1))-b), b); + SBasisOf r; + + for(int i = a.size()-1; i >= 0; i--) { + r = multiply_add(r, s, SBasisOf(LinearOf(HatOf(a[i][0]))) - b*a[i][0] + b*a[i][1]); + } + return r; +} + +template +SBasisOf compose(SBasisOf const &a, SBasisOf const &b, unsigned k){ + SBasisOf s = multiply((SBasisOf(LinearOf(1,1))-b), b); + SBasisOf r; + + for(int i = a.size()-1; i >= 0; i--) { + r = multiply_add(r, s, SBasisOf(LinearOf(HatOf(a[i][0]))) - b*a[i][0] + b*a[i][1]); + } + r.truncate(k); + return r; +} +template +SBasisOf compose(LinearOf const &a, SBasisOf const &b){ + return compose(SBasisOf(a),b); +} +template +SBasisOf compose(SBasisOf const &a, LinearOf const &b){ + return compose(a,SBasisOf(b)); +} +template//TODO: don't be so lazy!! +SBasisOf compose(LinearOf const &a, LinearOf const &b){ + return compose(SBasisOf(a),SBasisOf(b)); +} + + + +template +inline SBasisOf portion(const SBasisOf &t, double from, double to) { return compose(t, LinearOf(from, to)); } + +// compute f(g) +template +inline SBasisOf +SBasisOf::operator()(SBasisOf const & g) const { + return compose(*this, g); +} + +template +inline std::ostream &operator<< (std::ostream &out_file, const LinearOf &bo) { + out_file << "{" << bo[0] << ", " << bo[1] << "}"; + return out_file; +} + +template +inline std::ostream &operator<< (std::ostream &out_file, const SBasisOf & p) { + for(unsigned i = 0; i < p.size(); i++) { + out_file << p[i] << "s^" << i << " + "; + } + return out_file; +} + +}; +#endif + + +/* + 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:fileencoding=utf-8:textwidth=99 : -- cgit v1.2.3