diff options
Diffstat (limited to 'include/2geom/orphan-code/sbasisN.h')
| -rw-r--r-- | include/2geom/orphan-code/sbasisN.h | 1123 |
1 files changed, 1123 insertions, 0 deletions
diff --git a/include/2geom/orphan-code/sbasisN.h b/include/2geom/orphan-code/sbasisN.h new file mode 100644 index 0000000..0b5a48f --- /dev/null +++ b/include/2geom/orphan-code/sbasisN.h @@ -0,0 +1,1123 @@ +/**
+ * \file
+ * \brief Multi-dimensional symmetric power basis function.
+ * A SBasisN<n> is a polynomial f of n variables (t0,...,tn-1),
+ * written in a particular form. Let si = ti(1-t_i). f is written as
+ *
+ * f = sum_p s^p a_{p}(t0,...,t_{n-1})
+ *
+ * where p=(p0,...,p_{n-1}) is a multi index (called MultiDegree<n> in the code)
+ * s^p = prod_i si^pi, and a_p is a LinearN<n>.
+ * Recall a LinearN<n> is sum over all choices xi = ti or (1-ti) of terms of form
+ * a * x0*...*x_{n-1}
+ *
+ * Caution: degrees are expressed as degrees of s=t*(1-t). The real degree
+ * (with respect to t) of the polynomial is twice that + 0 or 1 depending
+ * whether the relevant LinearN<n> coeff is constant or not.
+ *//*
+ *
+ * Authors:
+ * JF Barraud <jf.barraud@gmail.com>
+ * Nathan Hurst <njh@mail.csse.monash.edu.au>
+ * Michael Sloan <mgsloan@gmail.com>
+ *
+ * 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_SBASISN_H
+#define SEEN_SBASISN_H
+#include <vector>
+#include <cassert>
+#include <iostream>
+
+#include <2geom/orphan-code/linearN.h>
+#include <2geom/linear.h>//for conversion purpose
+#include <2geom/sbasis.h>//for conversion purpose
+#include <2geom/interval.h>
+#include <2geom/utils.h>
+#include <2geom/exception.h>
+
+
+namespace Geom{
+
+/** MultiDegree:
+ * \brief multi-degree (p0,...,p^{n-1}) of a s0^p0...s_{n-1}p^{n-1} monomial.
+ *
+ * a "Multi_deg" is a sequence p={p0,...,p_n-1} representing the monomial
+ * s^p = s_0^{p_0}*...*s_{n-1}^{p_{n-1}}.
+ * Caution: the degrees are expressed with respect to si! ( in SBasis code
+ * below, si = ti*(1-ti) ).
+ */
+
+template<unsigned n>
+class MultiDegree{
+public:
+ unsigned p[n];
+ MultiDegree(){
+ for (unsigned i = 0; i <n; i++) {
+ p[i]=0;
+ }
+ }
+ MultiDegree( unsigned const other_p[] ){
+ p = other_p;
+ }
+ MultiDegree(unsigned const idx, unsigned const sizes[]){
+ unsigned q = idx;
+ for (unsigned i = n-1; i >0; i--) {
+ div_t d = std::div(int(q), int(sizes[i]));
+ p[i] = d.rem;
+ q = d.quot;
+ }
+ p[0] = q;
+ }
+ unsigned operator[](const unsigned i) const {
+ assert(i < n); return p[i];
+ }
+ unsigned& operator[](const unsigned i) {
+ assert(i < n); return p[i];
+ }
+
+ unsigned asIdx(unsigned const sizes[]) const{
+ unsigned ret = p[0];
+ bool in_range = (p[0]<sizes[0]);
+ for (unsigned i = 1; i < n; i++) {
+ in_range = in_range && (p[i]<sizes[i]);
+ ret = ret*sizes[i] + p[i];
+ }
+ if (in_range) return ret;
+ //TODO: find a better warning than returning out of range idx!
+ ret =1;
+ for (unsigned i = 0; i < n; i++) {
+ ret *= sizes[i];
+ }
+ return ret;
+ }
+ bool stepUp(unsigned const sizes[], unsigned frozen_mask = 0){
+ unsigned i = 0;
+ while ( i < n && ( (1<<i) & frozen_mask ) ) i++;
+ while ( i <n && p[i] == sizes[i]-1 ) {
+ i++;
+ while (i<n && ( (1<<i) & frozen_mask ) ) i++;
+ }
+ if (i<n){
+ p[i]+=1;
+ for (unsigned j = 0; j < i; j++) {
+ if ( !( (1<<j) & frozen_mask ) ) p[j] = 0;
+ }
+ return true;
+ }else{
+ return false;
+ }
+ }
+ bool stepDown(unsigned const sizes[], unsigned frozen_mask = 0){
+ int i = n-1;
+ while (i>=0 && ( (1<<i) & frozen_mask ) ) i--;
+ while ( i >= 0 && p[i] == 0 ) {
+ i--;
+ while (i>=0 && ( (1<<i) & frozen_mask ) ) i--;
+ }
+ if ( i >= 0 ){
+ p[i]-=1;
+ for (unsigned j = i+1; j < n; j++) {
+ if ( !( (1<<j) & frozen_mask ) ) p[j] = sizes[j]-1;
+ }
+ return true;
+ }else{
+ return false;
+ }
+ }
+};
+
+/**
+ * Returns the maximal degree appearing in the two arguments for each variables.
+ */
+template <unsigned n>
+MultiDegree<n> max(MultiDegree<n> const &p, MultiDegree<n> const &q){
+ MultiDegree<n> ret;
+ for (unsigned i = 0; i <n; i++) {
+ ret.p[i] = (p[i]>q[i] ? p[i] : q[i]);
+ }
+ return ret;
+}
+
+template <unsigned n>
+MultiDegree<n> operator + (MultiDegree<n> const &p, MultiDegree<n> const &q){
+ MultiDegree<n> ret;
+ for (unsigned i = 0; i <n; i++) {
+ ret.p[i] = p[i] + q[i];
+ }
+ return ret;
+}
+template <unsigned n>
+MultiDegree<n> operator += (MultiDegree<n> const &p, MultiDegree<n> const &q){
+ for (unsigned i = 0; i <n; i++) {
+ p[i] += q[i];
+ }
+ return p;
+}
+
+/**
+ * \brief MultiDegree comparison.
+ * A MultiDegree \param p is smaller than another \param q
+ * if all it's smaller for all variables.
+ *
+ * In particular, p<=q and q<=p can both be false!
+ */
+template<unsigned n>
+bool operator<=(MultiDegree<n> const &p, MultiDegree<n> const &q){
+ for (unsigned i = 0; i <n; i++) {
+ if (p[i]>q[i]) return false;
+ }
+ return true;
+}
+
+
+/**
+ * \brief Polynomials of n variables, written in SBasis form.
+ * An SBasisN<n> f is internaly represented as a vector of LinearN<n>.
+ * It should be thought of as an n-dimensional vector: the coef of s0^p0...s_{n-1}p^{n-1}
+ * is soterd in f[p0,...,p_{n-1}]. The sizes of each dimension is stored in "sizes".
+ * Note: an empty SBasis is identically 0.
+ */
+template<unsigned n>
+class SBasisN : public std::vector<LinearN<n> >{
+public:
+ unsigned sizes[n];
+ SBasisN() {
+ for (unsigned i = 0; i < n; i++) {
+ sizes[i] = 0;
+ }
+ }
+ explicit SBasisN(double a) {
+ for (unsigned i = 0; i < n; i++) {
+ sizes[i] = 1;
+ }
+ this->push_back(LinearN<n>(a));
+ }
+ SBasisN(SBasisN<n> const & a) : std::vector<LinearN<n> >(a){
+ //TODO: efficient array copy??
+ for (unsigned i = 0; i < n; i++) {
+ sizes[i] = a.sizes[i];
+ }
+ }
+ SBasisN(LinearN<n> const & bo) {
+ for (unsigned i = 0; i < n; i++) {
+ sizes[i] = 1;
+ }
+ this->push_back(bo);
+ }
+ SBasisN(LinearN<n>* bo) {
+ for (unsigned i = 0; i < n; i++) {
+ sizes[i] = 1;
+ }
+ this->push_back(*bo);
+ }
+
+//----------------------------------------------
+//-- Degree/Sizing facilities ------------------
+//----------------------------------------------
+/**
+ * Internal recursive function used to compute partial degrees.
+ */
+ bool find_non_empty_level(unsigned var, MultiDegree<n> &fixed_degrees)const{
+ if (this->size()==0){
+ for (unsigned i = 0; i < n; i++) {
+ fixed_degrees[i] = 0;//FIXME this should be -infty!!
+ }
+ return false;
+ }
+ if ( !((*this)[fixed_degrees.asIdx(sizes)].isZero()) ) return true;
+
+ unsigned frozen = (1<<var);
+ if ( fixed_degrees.stepDown(sizes, frozen) ){
+ if ( find_non_empty_level(var, fixed_degrees) ) return true;
+ }
+ if ( fixed_degrees[var] > 0 ){
+ fixed_degrees[var] -= 1;
+ for (unsigned i = 0; i < n; i++) {
+ if (i!=var) fixed_degrees[i] = sizes[i]-1;
+ }
+ if (find_non_empty_level(var, fixed_degrees)) return true;
+ }
+ return false;//FIXME: this should return -infty in all variables!
+ }
+
+/**
+ * Returns the degree of an SBasisN<n> with respect to a given variable form its sizes.
+ * All terms are taken into account, even eventual trailing zeros.
+ * Note: degree is expressed with respect to s = t*(1-t), not t itself.
+ */
+ unsigned quick_degree(unsigned var) const{
+ return ( sizes[var] > 0 ? sizes[var]-1 : 0 );//this should be -infty.
+ }
+/**
+ * Computes the multi degree of the SBasis from it's sizes.
+ * All terms are taken into account, even eventual trailing zeros.
+ * Note: degrees are expressed with respect to s = t*(1-t), not t itself.
+ */
+ MultiDegree<n> quick_multi_degree() const{
+ MultiDegree<n> ret;
+ if (this->size()==0) return ret;//should be -infty for all vars.
+ for (unsigned i = 0; i < n; i++) {
+ assert( sizes[i]>0 );
+ ret.p[i] = sizes[i]-1;
+ }
+ return ret;
+ }
+/**
+ * Returns the degree of an SBasisN<n> with respect to a given variable.
+ * Trailing zeros are not taken into account.
+ * Note: degree is expressed with respect to s = t*(1-t), not t itself.
+ */
+ unsigned degree(unsigned var)const{
+ MultiDegree<n> degrees;
+ for(unsigned i = 0; i < n; i++) {
+ degrees[i] = sizes[i]-1;
+ }
+ if ( find_non_empty_level(var, degrees) ) return degrees[var];
+ else return 0;//should be -infty.
+ }
+/**
+ * Returns the *real* degree of an SBasisN<n> with respect to a given variable.
+ * Trailing zeros are not taken into account.
+ * Note: degree is expressed with respect to t itself, not s = t*(1-t).
+ * In particular: real_t_degree() = 2*degree() + 0 or 1.
+ */
+ unsigned real_t_degree(unsigned var)const{
+ unsigned deg = 0;
+ bool even = true;
+ bool notfound = true;
+ unsigned frozen = (1<<var);
+ MultiDegree<n> degrees;
+ for(unsigned i = 0; i < n; i++) {
+ degrees[i] = sizes[i]-1;
+ }
+ while( notfound ){
+ if ( find_non_empty_level(var, degrees) && degrees[var]>= deg ){
+ deg = degrees[var];
+ even = (*this)[degrees.asIdx(sizes)].isConstant(var);
+ }
+ notfound = even && degrees.stepDown(sizes, frozen);
+ }
+ return 2*deg + ( even ? 0 : 1 );
+ }
+/**
+ * Returns the *real* degrees of an SBasisN<n>.
+ * Trailing zeros are not taken into account.
+ * Note: degree is expressed with respect to t itself, not s = t*(1-t).
+ * In particular: real_t_degree() = 2*degree() + 0 or 1.
+ */
+ MultiDegree<n> real_t_degrees()const{
+ MultiDegree<n>res;
+ for(unsigned i = 0; i < n; i++) {
+ res[i] = real_t_degree(i);
+ }
+ return res;
+ }
+/**
+ * Computes the multi degree of the SBasis.
+ * Trailing zeros are not taken into account.
+ * Note: degree is expressed with respect to s = t*(1-t), not t itself.
+ */
+ MultiDegree<n> multi_degree() const{
+ MultiDegree<n> ret;
+ if (this->size()==0) return ret;//should be -infty for all vars.
+ for (unsigned i = 0; i < n; i++) {
+ ret[i] = this->degree(i);
+ }
+ return ret;
+ }
+/**
+ * Returns the highest degree over all variables.
+ * Note: degree is expressed with respect to s = t*(1-t), not t itself.
+ */
+ unsigned max_degree() const {
+ if (this->size()==0) return 0;//should be -infty!
+ unsigned d=0;
+ for (unsigned i = 0; i < n; i++) {
+ assert( sizes[i]>0 );
+ if (d < sizes[i]-1) d = sizes[i]-1;
+ }
+ return d;
+ }
+
+/**
+ * Resize an SBasisN<n> to match new sizes.
+ *
+ * Caution: if a new size is smaller, the corresponding coefficients are discarded.
+ */
+ void multi_resize(unsigned new_sizes[], LinearN<n> def_value = LinearN<n>(0.)){
+ SBasisN<n> result;
+ bool nothing_todo = true;
+ unsigned tot_size = 1;
+ for(unsigned i = 0; i < n; i++) {
+ nothing_todo = nothing_todo && (sizes[i] == new_sizes[i]);
+ result.sizes[i] = new_sizes[i];
+ tot_size *= new_sizes[i];
+ }
+ if (nothing_todo) return;
+ result.resize(tot_size, def_value);
+ for(unsigned i = 0; i < tot_size; i++) {
+ MultiDegree<n> d( i, result.sizes );
+ unsigned j = d.asIdx(sizes);
+ if ( j < this->size() ){
+ result[i] = (*this)[j];
+ }
+ }
+ *this = result;
+ }
+
+ //remove extra zeros
+ void normalize() {
+ MultiDegree<n> max_p = multi_degree();
+ unsigned new_sizes[n];
+ for (unsigned i=0; i<n; i++){
+ new_sizes[i] = max_p[i]+1;
+ }
+ multi_resize(new_sizes);
+ }
+
+//-----------------------------
+//-- Misc. --------------------
+//-----------------------------
+
+/**
+ * Returns the number of variables this function takes as input: n.
+ */
+ unsigned input_dim(){return n;};
+
+ //IMPL: FragmentConcept
+ typedef double output_type;
+
+ inline bool isZero() const {
+ if(this->size()==0) return true;
+ for(unsigned i = 0; i < this->size(); i++) {
+ if(!(*this)[i].isZero()) return false;
+ }
+ return true;
+ }
+ inline bool isConstant() const {
+ if (this->size()==0) return true;
+ if(!(*this)[0].isConstant()) return false;
+ for (unsigned i = 1; i < this->size(); i++) {
+ if(!(*this)[i].isZero()) return false;
+ }
+ return true;
+ }
+
+ bool isFinite() const{
+ for (unsigned i = 0; i < this->size(); i++) {
+ if(!(*this)[i].isFinite()) return false;
+ }
+ return true;
+ }
+
+
+//------------------------------------------
+//-- Evaluation methods --------------------
+//------------------------------------------
+/**
+ * Returns the value of the SBasis at a given corner of [0,1]^n.
+ * \param k describes the corner: if i-th bit is 0, ti=0, otherwise ti=1.
+ */
+ inline double atCorner(unsigned k) const {
+ if(this->size()==0) return 0.;
+ return (*this)[0].atCorner(k);
+ }
+/**
+ * Returns the value of the SBasis at a given corner of [0,1]^n.
+ * \param t[n] describes the corner: the values should be 0's and 1's.
+ */
+ inline double atCorner(double t[]) const {
+ if(this->size()==0) return 0.;
+ return (*this)[0].atCorner(t);
+ }
+/**
+ * Returns a "slice" of the array.
+ * Returns an SBasis containing all the coeff of (s-)degree \param deg in variable \param var
+ */
+ //TODO: move by bigger blocks (?) but they are broken into pieces in general...
+ SBasisN<n> slice(unsigned const var, unsigned const deg) const{
+ if (deg >= sizes[var] ) return SBasisN<n>();
+ SBasisN<n> res;
+ unsigned tot_size = 1;
+ for (unsigned i = 0; i < n; i++) {
+ res.sizes[i] = (i==var ? 1 : sizes[i]);
+ tot_size *= res.sizes[i];
+ }
+ res.resize( tot_size, LinearN<n>(0.));
+ for (unsigned i = 0; i < tot_size; i++) {
+ MultiDegree<n> d(i,res.sizes);
+ d.p[var] = deg;
+ res[i] = (*this)[d.asIdx(sizes)];
+ }
+ return res;
+ }
+/**
+ * Returns a the SBasisN<n-1> obtained by setting variable \param var to 0.
+ */
+ inline SBasisN<n-1> at0(unsigned var=0, unsigned deg=0) const {
+ SBasisN<n> sl = slice(var,deg);
+ SBasisN<n-1> res;
+ res.reserve(sl.size());
+ for (unsigned i = 0; i < n-1; i++) {
+ res.sizes[i] = sizes[ ( i<var ? i : i+1 ) ];
+ }
+ for (unsigned i = 0; i < sl.size(); i++) {
+ res.push_back( sl[i].at0(var) );
+ }
+ return res;
+ }
+/**
+ * Returns a the SBasisN<n-1> obtained by setting variable \param var to 1.
+ */
+ inline SBasisN<n-1> at1(unsigned var=0, unsigned deg=0) const {
+ SBasisN<n> sl = slice(var,deg);
+ SBasisN<n-1> res;
+ res.reserve(sl.size());
+ for (unsigned i = 0; i < n-1; i++) {
+ res.sizes[i] = sizes[ ( i<var ? i : i+1 ) ];
+ }
+ for (unsigned i = 0; i < sl.size(); i++) {
+ res.push_back( sl[i].at1(var) );
+ }
+ return res;
+ }
+/**
+ * Returns a the SBasisN<n-1> obtained by setting variable \param var to \param t.
+ */
+ inline SBasisN<n-1> partialEval(double t, unsigned var=0 ) const {
+ SBasisN<n> sl;
+ double s = t*(1-t);
+ double si = 1;
+ for (unsigned i = 0; i <sizes[var]; i++) {
+ sl = sl + slice(var, i)*si;
+ si *= s;
+ }
+ SBasisN<n-1> res;
+ res.resize(sl.size(), LinearN<n-1>(0.));
+ for (unsigned i = 0; i < n-1; i++) {
+ res.sizes[i] = sizes[ ( i<var ? i : i+1 ) ];
+ }
+ for (unsigned i = 0; i < sl.size(); i++) {
+ res[i] = sl[i].partialEval(t,var);
+ }
+ return res;
+ }
+
+/**
+ * \brief Internal recursive function.
+ * Replace each variable by it's value in the 's=t*(1-t)' factor
+ * but not in the LinearN<n> coeffs. Then sum up all coefficients.
+ * \param t[n]: values of the variables.
+ */
+ LinearN<n> sumCoefs( double t[], unsigned const k, unsigned const idx) const{
+ LinearN<n> a;
+ if (k == n){
+ a = (*this)[idx];
+ return (*this)[idx];
+ }
+ double s = t[k]*(1-t[k]);
+ double si=1;
+ for (unsigned i=0; i<sizes[k]; i++){
+ a += sumCoefs(t,k+1,idx*sizes[k]+i)*si;;
+ si *= s;
+ }
+ return a;
+ }
+/**
+ * Evaluate at given n-dimensional point.
+ * \param t[n]: values of the variables.
+ */
+ double valueAt(double t[]) const {
+ LinearN<n> a = sumCoefs(t,0,0);
+ return a.valueAt(t);
+ }
+
+ double operator()(double t[]) const {
+ return valueAt(t);
+ }
+
+ //double valueAndDerivative(double t, double &der) const;
+ //std::vector<double> valueAndDerivatives(double t, unsigned n) const;
+ //SBasisN toSBasisN() const { return SBasisN(*this); }
+ //double tailError(unsigned tail) const;
+
+
+//--------------------------------------------------
+//-- Coeff. manipulation ---------------------------
+//--------------------------------------------------
+
+/**
+ * Accessing the SBasisN<n> coefficients.
+ */
+ LinearN<n> operator[](unsigned i) const {
+ assert(i < this->size());
+ return std::vector<LinearN<n> >::operator[](i);
+ }
+ LinearN<n> operator[](MultiDegree<n> const &p) const {
+ unsigned i = p.asIdx(sizes);
+ assert(i < this->size());
+ return std::vector<LinearN<n> >::operator[](i);
+ }
+
+//MUTATOR PRISON
+ LinearN<n>& operator[](unsigned i) { return this->at(i); }
+// LinearN<n>& operator[](MultiDegree const &p) {
+// unsigned i = p.asIdx(sizes);
+// return this->at(i);
+// }
+
+ void appendCoef(const SBasisN<n-1> &a, const SBasisN<n-1> &b, unsigned var=0){
+ unsigned new_sizes[n];
+ MultiDegree<n-1> deg_a = a.multi_degree(), deg_b = b.multi_degree();
+ MultiDegree<n-1> dcoef = max( deg_a, deg_b );
+ for (unsigned i=0; i<n; i++){
+ if ( i == var ){
+ new_sizes[var] = sizes[var] + 1;
+ }else{
+ unsigned coef_size = dcoef[(i<var?i:i-1)] + 1;
+ new_sizes[i] = ( sizes[i]>coef_size ? sizes[i] : coef_size );
+ }
+ }
+ multi_resize(new_sizes);
+
+ MultiDegree<n> d;
+ d[var] = sizes[var]-1;
+ unsigned frozen_mask = (1<<var);
+ do{
+ for (unsigned i=0; i<n-1; i++){
+ dcoef.p[i] = d.p[ ( i<var ? i : i+1) ];
+ }
+ LinearN<n-1> a_d,b_d;
+ unsigned ia = dcoef.asIdx(a.sizes);
+ if ( ia < a.size() ) a_d = a[ia];
+ unsigned ib = dcoef.asIdx(b.sizes);
+ if ( ib < b.size() ) b_d = b[ib];
+ (*this)[d.asIdx(sizes)] = LinearN<n>(a_d,b_d);
+ }while (d.stepUp(sizes,frozen_mask));
+ }
+
+//private:
+ //void derive(); // in place version
+};
+
+//SBasisN<0> is a double. Specialize it out.
+template<>
+class SBasisN<0>{
+public:
+ double d;
+ SBasisN () {}
+ SBasisN(double d) :d(d) {}
+ operator double() const { return d; }
+};
+
+
+//SBasisN<1> are usual SBasis. Allow conversion.
+SBasis toSBasis(SBasisN<1> f){
+ SBasis res(f.size(), Linear());
+ for (unsigned i = 0; i < f.size(); i++) {
+ res[i] = toLinear(f[i]);
+ }
+ return res;
+}
+
+//TODO: figure out how to stick this in linear, while not adding an sbasis dep
+template<unsigned n>
+inline SBasisN<n> LinearN<n>::toSBasisN() const { return SBasisN<n>(*this); }
+
+
+
+
+//implemented in sbasis-roots.cpp
+//OptInterval bounds_exact(SBasisN const &a);
+//OptInterval bounds_fast(SBasisN const &a, int order = 0);
+//OptInterval bounds_local(SBasisN const &a, const OptInterval &t, int order = 0);
+
+/** Returns a function which reverses the domain of a.
+ \param a sbasis function
+
+useful for reversing a parameteric curve.
+*/
+//template<unsigned n>
+//inline SBasisN<n> reverse(SBasisN<n> const &a);
+
+//IMPL: ScalableConcept
+template<unsigned n>
+inline SBasisN<n> operator-(const SBasisN<n>& p) {
+ if(p.isZero()) return SBasisN<n>();
+ SBasisN<n> result;
+ for(unsigned i = 0; i < n; i++) {
+ result.sizes[i] = p.sizes[i];
+ }
+ result.reserve(p.size());
+ for(unsigned i = 0; i < p.size(); i++) {
+ result.push_back(-p[i]);
+ }
+ return result;
+}
+template<unsigned n>
+SBasisN<n> operator*(SBasisN<n> const &a, double c){
+ if(a.isZero()) return SBasisN<n>();
+ SBasisN<n> result;
+ for(unsigned i = 0; i < n; i++) {
+ result.sizes[i] = a.sizes[i];
+ }
+ result.reserve(a.size());
+ for(unsigned i = 0; i < a.size(); i++) {
+ result.push_back(a[i] * c);
+ }
+ return result;
+}
+template<unsigned n>
+inline SBasisN<n> operator*(double k, SBasisN<n> const &a) { return a*k; }
+template<unsigned n>
+inline SBasisN<n> operator/(SBasisN<n> const &a, double k) { return a*(1./k); }
+template<unsigned n>
+SBasisN<n>& operator*=(SBasisN<n>& a, double c){
+ for(unsigned i = 0; i < a.size(); i++) a[i] *= c;
+ return a;
+}
+template<unsigned n>
+inline SBasisN<n>& operator/=(SBasisN<n>& a, double b) { return (a*=(1./b)); }
+
+//IMPL: AddableConcept
+template<unsigned n>
+SBasisN<n> operator + (const SBasisN<n>& a, const SBasisN<n>& b){
+ if( a.isZero() ) return b;
+ if( b.isZero() ) return a;
+ SBasisN<n> result;
+ MultiDegree<n> deg = max(a.quick_multi_degree(),b.quick_multi_degree());
+ unsigned max_size = 1;
+ for(unsigned i = 0; i < n; i++) {
+ result.sizes[i] = deg[i]+1;
+ max_size *= result.sizes[i];
+ }
+ result.resize( max_size, LinearN<n>(0.) );
+ for(unsigned i = 0; i < result.size(); i++) {
+ MultiDegree<n> p(i,result.sizes);
+ unsigned ia = p.asIdx(a.sizes);
+ unsigned ib = p.asIdx(b.sizes);
+ if (ia<a.size()) {
+ result[i] += a[ia];
+ }
+ if (ib<b.size()) {
+ result[i] += b[ib];
+ }
+ }
+ return result;
+}
+template<unsigned n>
+SBasisN<n> operator-(const SBasisN<n>& a, const SBasisN<n>& b){return a+(-b);}
+template<unsigned n>
+SBasisN<n>& operator+=(SBasisN<n>& a, const SBasisN<n>& b){
+ if(b.isZero()) return a;
+ a = a + b;
+ return a;
+}
+template<unsigned n>
+SBasisN<n>& operator-=(SBasisN<n>& a, const SBasisN<n>& b){
+ a += -b;
+ return a;
+}
+
+//TODO: remove?
+template<unsigned n>
+inline SBasisN<n> operator+(const SBasisN<n> & a, LinearN<n> const & b) {
+ if(b.isZero()) return a;
+ if(a.isZero()) return b;
+ SBasisN<n> result(a);
+ result[0] += b;
+ return result;
+}
+template<unsigned n>
+
+inline SBasisN<n> operator-(const SBasisN<n> & a, LinearN<n> const & b) {
+ if(b.isZero()) return a;
+ if(a.isZero()) return -b;
+ SBasisN<n> result(a);
+ result[0] -= b;
+ return result;
+}
+template<unsigned n>
+inline SBasisN<n>& operator+=(SBasisN<n>& a, const LinearN<n>& b) {
+ if(a.size()==0)
+ a.push_back(b);
+ else
+ a[0] += b;
+ return a;
+}
+template<unsigned n>
+inline SBasisN<n>& operator-=(SBasisN<n>& a, const LinearN<n>& b) {
+ if(a.size()==0)
+ a.push_back(-b);
+ else
+ a[0] -= b;
+ return a;
+}
+
+//IMPL: OffsetableConcept
+template<unsigned n>
+inline SBasisN<n> operator+(const SBasisN<n> & a, double b) {
+ if(a.isZero()) return LinearN<n>(b);
+ SBasisN<n> result(a);
+ result[0] += b;
+ return result;
+}
+template<unsigned n>
+inline SBasisN<n> operator-(const SBasisN<n> & a, double b) {
+ if(a.isZero()) return LinearN<n>(-b);
+ SBasisN<n> result(a);
+ result[0] -= b;
+ return result;
+}
+template<unsigned n>
+inline SBasisN<n>& operator+=(SBasisN<n>& a, double b) {
+ if(a.size()==0)
+ a.push_back(LinearN<n>(b));
+ else
+ a[0] += b;
+ return a;
+}
+template<unsigned n>
+inline SBasisN<n>& operator-=(SBasisN<n>& a, double b) {
+ if(a.size()==0)
+ a.push_back(LinearN<n>(-b));
+ else
+ a[0] -= b;
+ return a;
+}
+
+template<unsigned n>
+SBasisN<n> shift(SBasisN<n> const &a, MultiDegree<n> sh){
+ SBasisN<n> result;
+ MultiDegree<n> deg = a.quick_multi_degree() + sh;
+ for(unsigned i = 0; i < n; i++) {
+ result.sizes[i] = deg[i]+1;
+ }
+ unsigned max_size = deg.asIdx(result.sizes);
+ result.resize( max_size, LinearN<n>(0.) );
+ for(unsigned i = 0; i < a.size(); i++) {
+ MultiDegree<n> p(i,a.sizes);
+ p+=sh;
+ result[p.asIdx(result.sizes)]=a[i];
+ }
+ return result;
+}
+template<unsigned n>
+SBasisN<n> shift(LinearN<n> const &a, MultiDegree<n> sh){
+ SBasisN<n> result;
+ for(unsigned i = 0; i < n; i++) {
+ result.sizes[i] = sh[i]+1;
+ }
+ unsigned max_size = sh.asIdx(result.sizes);
+ result.resize( max_size, LinearN<n>(0.) );
+ result[max_size-1]=a;
+ return result;
+}
+//shift only in one variable
+template<unsigned n>
+SBasisN<n> shift(LinearN<n> const &a, unsigned sh, unsigned var){
+ assert( var < n );
+ SBasisN<n> result;
+ for(unsigned i = 0; i < n; i++) {
+ result.sizes[i] = 1;
+ }
+ result.sizes[var] = sh+1;
+ result.resize( sh+1, LinearN<n>(0.) );
+ result[sh]=a;
+ return result;
+}
+
+//truncate only in first variable
+template<unsigned n>
+inline SBasisN<n> truncate(SBasisN<n> const &a, unsigned first_size) {
+ if ( first_size <= a.sizes[0] ) return a;
+ SBasisN<n> c;
+ for (unsigned i = 0; i < n; i++) {
+ c.sizes[i] = a.sizes[i];
+ }
+ c.sizes[0] = first_size;
+ unsigned tot_size = 1;
+ for(unsigned i = 0; i < n; i++) {
+ tot_size*=c.sizes[i];
+ }
+ c.insert(c.begin(), a.begin(), a.begin() + tot_size);
+ return c;
+}
+
+template<unsigned n>
+SBasisN<n> multiply(SBasisN<n> const &a, SBasisN<n> const &b){
+ SBasisN<n> c;
+ MultiDegree<n> d;
+ MultiDegree<n> t_deg = a.real_t_degrees() + b.real_t_degrees();
+ for(unsigned i = 0; i < n; i++) {
+ d[i] = ( t_deg[i]%2 == 0 ? t_deg[i]/2 : (t_deg[i]-1)/2 ) ;
+ }
+ unsigned new_sizes[n], tot_size = 1;
+ for(unsigned i = 0; i < n; i++) {
+ //c.sizes[i] = d[i] + 1+1;//product of linears might give 1 more s in each dir!!
+ new_sizes[i] = d[i] + 1;
+ tot_size*=new_sizes[i];
+ }
+ c.resize( tot_size, LinearN<n>(0.) );
+ for(unsigned i = 0; i < n; i++) {
+ c.sizes[i] = new_sizes[i];
+ }
+
+ for(unsigned ai = 0; ai < a.size(); ai++) {
+ for(unsigned bj = 0; bj < b.size(); bj++) {
+ MultiDegree<n> di( ai, a.sizes );
+ MultiDegree<n> dj( bj, b.sizes );
+ //compute a[ai]*b[bj]:
+ for(unsigned p = 0; p < (1<<n); p++) {
+ for(unsigned q = 0; q < (1<<n); q++) {
+
+ //compute a[ai][p]*b[bj][q]:
+ unsigned m = p^q;//m has ones for factors s, 0 for (t-s) or ((1-t)-s).
+ for(unsigned r = 0; r < (1<<n); r++) {
+ if (!(r&m)) {// a 1 in r means take t (or (1-t)), otherwise take -s.
+ int sign = 1;
+ MultiDegree<n> dr;
+ unsigned t0 = 0, t1 = 0;
+ for (unsigned var = 0; var < n; var++) {
+ //if var is in mask m, no choice, take s
+ if ( m & (1<<var) ){
+ dr.p[var] = 1;
+ }//if var is in mask r, take t or (1-t)
+ else if ( r & (1<<var) ){
+ dr.p[var] = 0;
+ if ( p&(1<<var) ) {
+ t0 = t0 | (1<<var);
+ }else{
+ t1 = t1 | (1<<var);
+ }
+ }//ohterwise take -s
+ else{
+ dr.p[var] = 1;
+ sign *= -1;
+ }
+ }
+ unsigned idx = (di+dj+dr).asIdx(c.sizes);
+ if (idx < c.size()){
+ for(unsigned s = 0; s < (1<<n); s++) {
+ if ( (t0 & ~s) || (t1 & s) ){
+ c[idx][s] += 0;
+ }else{
+ c[idx][s] += sign * a[ai][p] * b[bj][q];
+ }
+ }
+ }
+ }
+ }//r loop: all choices have been expanded in the product a[ai][p]*b[bj][q]
+ }//q loop
+ }//p loop: all products a[ai][p]*b[bj][q] have been computed.
+ }//bj loop
+ }//ai loop: all a[ai]b[bj] have been computed.
+
+ //TODO: normalize c, or even better, compute with the right size from scratch
+ return c;
+}
+
+
+template<unsigned n>
+inline SBasisN<n> operator*(SBasisN<n> const & a, SBasisN<n> const & b) {
+ return multiply(a, b);
+}
+
+template<unsigned n>
+inline SBasisN<n>& operator*=(SBasisN<n>& a, SBasisN<n> const & b) {
+ a = multiply(a, b);
+ return a;
+}
+
+template<unsigned m,unsigned n>
+SBasisN<m> compose(LinearN<n> const &f, std::vector<SBasisN<m> > const &t, unsigned fixed=0, unsigned flags=0 ){
+ assert (t.size() == n );
+ if (fixed == n) {
+ return SBasisN<m>(1.) * f[flags];
+ }else{
+ SBasisN<m> a0 = compose(f, t, fixed+1, flags);
+ SBasisN<m> a1 = compose(f, t, fixed+1, flags|(1<<fixed));
+ return (-t[fixed]+1) * a0 + t[fixed] * a1;
+ }
+}
+
+template<unsigned m,unsigned n>
+SBasisN<m> compose(SBasisN<n> const &f, std::vector<SBasisN<m> > const &t, unsigned const k=0, unsigned const idx = 0){
+ assert (t.size() == n );
+ if (k == n){
+ return compose( f[idx], t);
+ }
+ SBasisN<m> a;
+ SBasisN<m> s = multiply( t[k], (-t[k]+1.) );
+ SBasisN<m> si= SBasisN<m>(1.);
+ for (unsigned i=0; i<f.sizes[k]; i++){
+ a += compose(f, t,k+1,idx*f.sizes[k]+i)*si;;
+ si *= s;
+ }
+ return a;
+}
+
+template <unsigned n>
+inline std::ostream &operator<< (std::ostream &out_file, const MultiDegree<n> & d) {
+ out_file << "s^{";
+ for(unsigned i = 0; i < n; i++) {
+ out_file << d[i] << (i == n-1 ? "}" : ",");
+ }
+ return out_file;
+}
+template <unsigned n>
+inline std::ostream &operator<< (std::ostream &out_file, const SBasisN<n> & p) {
+ for(unsigned i = 0; i < p.size(); i++) {
+ MultiDegree<n> d(i, p.sizes);
+ out_file << d << " " << p[i] << " + ";
+ }
+ return out_file;
+}
+
+
+//--------------------------------------------------
+//--------------------------------------------------
+//--------------------------------------------------
+//--------------------------------------------------
+//--------------------------------------------------
+
+#if 0
+
+
+// This performs a multiply and accumulate operation in about the same time as multiply. return a*b + c
+template<unsigned n>
+SBasisN<n> multiply_add(SBasisN<n> const &a, SBasisN<n> const &b, SBasisN<n> c);
+
+template<unsigned n>
+SBasisN<n> integral(SBasisN<n> const &c);
+template<unsigned n>
+SBasisN<n> derivative(SBasisN<n> const &a);
+
+template<unsigned n>
+SBasisN<n> sqrt(SBasisN<n> const &a, int k);
+
+// return a kth order approx to 1/a)
+template<unsigned n>
+SBasisN<n> reciprocal(LinearN<n> const &a, int k);
+template<unsigned n>
+SBasisN<n> divide(SBasisN<n> const &a, SBasisN<n> const &b, int k);
+
+
+/** Returns the degree of the first non zero coefficient.
+ \param a sbasis function
+ \param tol largest abs val considered 0
+ \returns first non zero coefficient
+*/
+template<unsigned n>
+inline unsigned
+valuation(SBasisN<n> const &a, double tol=0){
+ unsigned val=0;
+ while( val<a.size() &&
+ fabs(a[val][0])<tol &&
+ fabs(a[val][1])<tol )
+ val++;
+ return val;
+}
+
+// a(b(t))
+template<unsigned n>
+SBasisN<n> compose(SBasisN<n> const &a, SBasisN<n> const &b);
+template<unsigned n>
+SBasisN<n> compose(SBasisN<n> const &a, SBasisN<n> const &b, unsigned k);
+template<unsigned n>
+SBasisN<n> inverse(SBasisN<n> a, int k);
+//compose_inverse(f,g)=compose(f,inverse(g)), but is numerically more stable in some good cases...
+//TODO: requires g(0)=0 & g(1)=1 atm. generalization should be obvious.
+template<unsigned n>
+SBasisN<n> compose_inverse(SBasisN<n> const &f, SBasisN<n> const &g, unsigned order=2, double tol=1e-3);
+
+/** Returns the sbasis on domain [0,1] that was t on [from, to]
+ \param a sbasis function
+ \param from,to interval
+ \returns sbasis
+
+*/
+template<unsigned n>
+inline SBasisN<n> portion(const SBasisN<n> &t, double from, double to) { return compose(t, LinearN<n>(from, to)); }
+
+// compute f(g)
+template<unsigned n>
+inline SBasisN<n>
+SBasisN<n>::operator()(SBasisN<n> const & g) const {
+ return compose(*this, g);
+}
+
+template<unsigned n>
+inline std::ostream &operator<< (std::ostream &out_file, const LinearN<n> &bo) {
+ out_file << "{" << bo[0] << ", " << bo[1] << "}";
+ return out_file;
+}
+
+template<unsigned n>
+inline std::ostream &operator<< (std::ostream &out_file, const SBasisN<n> & p) {
+ for(unsigned i = 0; i < p.size(); i++) {
+ out_file << p[i] << "s^" << i << " + ";
+ }
+ return out_file;
+}
+
+// These are deprecated, use sbasis-math.h versions if possible
+template<unsigned n>
+SBasisN<n> sin(LinearN<n> bo, int k);
+template<unsigned n>
+SBasisN<n> cos(LinearN<n> bo, int k);
+
+template<unsigned n>
+std::vector<double> roots(SBasisN<n> const & s);
+template<unsigned n>
+std::vector<std::vector<double> > multi_roots(SBasisN<n> const &f,
+ std::vector<double> const &levels,
+ double htol=1e-7,
+ double vtol=1e-7,
+ double a=0,
+ double b=1);
+
+#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 :
+#endif
|