src/2geom/sbasis.h

	sbasis.h revision e6bdf746e2d9e775704a475a29cc1bb167ec271c
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b/*
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b *  sbasis.h - S-power basis function class
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b *
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b *  Authors:
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b *   Nathan Hurst <njh@mail.csse.monash.edu.au>
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b *   Michael Sloan <mgsloan@gmail.com>
2eefc362ae6a0a94b84ee5bc9e7844ef45c3642cAlex Valavanis *
2eefc362ae6a0a94b84ee5bc9e7844ef45c3642cAlex Valavanis * Copyright (C) 2006-2007 authors
2eefc362ae6a0a94b84ee5bc9e7844ef45c3642cAlex Valavanis *
2eefc362ae6a0a94b84ee5bc9e7844ef45c3642cAlex Valavanis * This library is free software; you can redistribute it and/or
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b * modify it either under the terms of the GNU Lesser General Public
2eefc362ae6a0a94b84ee5bc9e7844ef45c3642cAlex Valavanis * License version 2.1 as published by the Free Software Foundation
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b * (the "LGPL") or, at your option, under the terms of the Mozilla
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b * Public License Version 1.1 (the "MPL"). If you do not alter this
2eefc362ae6a0a94b84ee5bc9e7844ef45c3642cAlex Valavanis * notice, a recipient may use your version of this file under either
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b * the MPL or the LGPL.
2eefc362ae6a0a94b84ee5bc9e7844ef45c3642cAlex Valavanis *
2eefc362ae6a0a94b84ee5bc9e7844ef45c3642cAlex Valavanis * You should have received a copy of the LGPL along with this library
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b * in the file COPYING-LGPL-2.1; if not, write to the Free Software
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b * You should have received a copy of the MPL along with this library
2eefc362ae6a0a94b84ee5bc9e7844ef45c3642cAlex Valavanis * in the file COPYING-MPL-1.1
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b *
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b * The contents of this file are subject to the Mozilla Public License
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b * Version 1.1 (the "License"); you may not use this file except in
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b * compliance with the License. You may obtain a copy of the License at
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b * http://www.mozilla.org/MPL/
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b *
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b * OF ANY KIND, either express or implied. See the LGPL or the MPL for
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b * the specific language governing rights and limitations.
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b */
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b#ifndef SEEN_SBASIS_H
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b#define SEEN_SBASIS_H
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b#include <vector>
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b#include <cassert>
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b#include <iostream>
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b#include "linear.h"
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b#include "interval.h"
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b#include "utils.h"
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_bnamespace Geom {
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b/*** An empty SBasis is identically 0. */
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_bclass SBasis : public std::vector<Linear>{
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_bpublic:
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b    SBasis() {}
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b    explicit SBasis(double a) {
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b        push_back(Linear(a,a));
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b    }
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b    SBasis(SBasis const & a) :
f6d8c168c2d1b1e168545b090dfbb22a3f4923a3Alex Valavanis        std::vector<Linear>(a)
f6d8c168c2d1b1e168545b090dfbb22a3f4923a3Alex Valavanis    {}
7a10406c7610c65f3600986ded11320fe870ccb9Alex Valavanis    SBasis(Linear const & bo) {
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b        push_back(bo);
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b    }
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b    //IMPL: FragmentConcept
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b    typedef double output_type;
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b    inline bool isZero() const {
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b        if(empty()) return true;
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b        for(unsigned i = 0; i < size(); i++) {
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b            if(!(*this)[i].isZero()) return false;
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b        }
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b        return true;
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b    }
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b    bool isFinite() const;
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b    inline double at0() const {
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b        if(empty()) return 0; else return (*this)[0][0];
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b    }
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b    inline double at1() const{
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b        if(empty()) return 0; else return (*this)[0][1];
f6d8c168c2d1b1e168545b090dfbb22a3f4923a3Alex Valavanis    }
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b    double valueAt(double t) const {
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b        double s = t*(1-t);
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b        double p0 = 0, p1 = 0;
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b        double sk = 1;
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b//TODO: rewrite as horner
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b        for(unsigned k = 0; k < size(); k++) {
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b            p0 += sk*(*this)[k][0];
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b            p1 += sk*(*this)[k][1];
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b            sk *= s;
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b        }
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b        return (1-t)*p0 + t*p1;
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b    }
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b    double valueAndDerivative(double t, double &der) const {
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b        double s = t*(1-t);
d19de0ff8b0793e74e98429ded50868ce7e9f75cAlex Valavanis        double p0 = 0, p1 = 0;
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b        double sk = 1;
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b//TODO: rewrite as horner
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b        for(unsigned k = 0; k < size(); k++) {
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b            p0 += sk*(*this)[k][0];
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b            p1 += sk*(*this)[k][1];
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b            sk *= s;
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b        }
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b        // p0 and p1 at this point form a linear approximation at t
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b        der = p1 - p0;
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b        return (1-t)*p0 + t*p1;
f6d8c168c2d1b1e168545b090dfbb22a3f4923a3Alex Valavanis    }
f6d8c168c2d1b1e168545b090dfbb22a3f4923a3Alex Valavanis    double operator()(double t) const {
f6d8c168c2d1b1e168545b090dfbb22a3f4923a3Alex Valavanis        return valueAt(t);
f6d8c168c2d1b1e168545b090dfbb22a3f4923a3Alex Valavanis    }
f6d8c168c2d1b1e168545b090dfbb22a3f4923a3Alex Valavanis
f6d8c168c2d1b1e168545b090dfbb22a3f4923a3Alex Valavanis    std::vector<double> valueAndDerivatives(double /*t*/, unsigned /*n*/) const {
f6d8c168c2d1b1e168545b090dfbb22a3f4923a3Alex Valavanis        //TODO
f6d8c168c2d1b1e168545b090dfbb22a3f4923a3Alex Valavanis        throw NotImplemented();
f6d8c168c2d1b1e168545b090dfbb22a3f4923a3Alex Valavanis    }
f6d8c168c2d1b1e168545b090dfbb22a3f4923a3Alex Valavanis
f6d8c168c2d1b1e168545b090dfbb22a3f4923a3Alex Valavanis    SBasis toSBasis() const { return SBasis(*this); }
f6d8c168c2d1b1e168545b090dfbb22a3f4923a3Alex Valavanis
f6d8c168c2d1b1e168545b090dfbb22a3f4923a3Alex Valavanis    double tailError(unsigned tail) const;
f6d8c168c2d1b1e168545b090dfbb22a3f4923a3Alex Valavanis
f6d8c168c2d1b1e168545b090dfbb22a3f4923a3Alex Valavanis// compute f(g)
f6d8c168c2d1b1e168545b090dfbb22a3f4923a3Alex Valavanis    SBasis operator()(SBasis const & g) const;
f6d8c168c2d1b1e168545b090dfbb22a3f4923a3Alex Valavanis
f6d8c168c2d1b1e168545b090dfbb22a3f4923a3Alex Valavanis    Linear operator[](unsigned i) const {
f6d8c168c2d1b1e168545b090dfbb22a3f4923a3Alex Valavanis        assert(i < size());
f6d8c168c2d1b1e168545b090dfbb22a3f4923a3Alex Valavanis        return std::vector<Linear>::operator[](i);
f6d8c168c2d1b1e168545b090dfbb22a3f4923a3Alex Valavanis    }
f6d8c168c2d1b1e168545b090dfbb22a3f4923a3Alex Valavanis
f6d8c168c2d1b1e168545b090dfbb22a3f4923a3Alex Valavanis//MUTATOR PRISON
f6d8c168c2d1b1e168545b090dfbb22a3f4923a3Alex Valavanis    Linear& operator[](unsigned i) { return this->at(i); }
f6d8c168c2d1b1e168545b090dfbb22a3f4923a3Alex Valavanis
f6d8c168c2d1b1e168545b090dfbb22a3f4923a3Alex Valavanis    //remove extra zeros
f6d8c168c2d1b1e168545b090dfbb22a3f4923a3Alex Valavanis    void normalize() {
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b        while(!empty() && 0 == back()[0] && 0 == back()[1])
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b            pop_back();
f6d8c168c2d1b1e168545b090dfbb22a3f4923a3Alex Valavanis    }
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b    void truncate(unsigned k) { if(k < size()) resize(k); }
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b};
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b//TODO: figure out how to stick this in linear, while not adding an sbasis dep
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_binline SBasis Linear::toSBasis() const { return SBasis(*this); }
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b//implemented in sbasis-roots.cpp
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_bInterval bounds_exact(SBasis const &a);
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_bInterval bounds_fast(SBasis const &a, int order = 0);
f6d8c168c2d1b1e168545b090dfbb22a3f4923a3Alex ValavanisInterval bounds_local(SBasis const &a, const Interval &t, int order = 0);
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b
f6d8c168c2d1b1e168545b090dfbb22a3f4923a3Alex Valavanisinline SBasis reverse(SBasis const &a) {
f6d8c168c2d1b1e168545b090dfbb22a3f4923a3Alex Valavanis    SBasis result;
f6d8c168c2d1b1e168545b090dfbb22a3f4923a3Alex Valavanis    result.reserve(a.size());
f6d8c168c2d1b1e168545b090dfbb22a3f4923a3Alex Valavanis    for(unsigned k = 0; k < a.size(); k++)
f6d8c168c2d1b1e168545b090dfbb22a3f4923a3Alex Valavanis       result.push_back(reverse(a[k]));
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b    return result;
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b}
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b//IMPL: ScalableConcept
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_binline SBasis operator-(const SBasis& p) {
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b    if(p.isZero()) return SBasis();
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b    SBasis result;
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b    result.reserve(p.size());
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b    for(unsigned i = 0; i < p.size(); i++) {
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b        result.push_back(-p[i]);
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b    }
f6d8c168c2d1b1e168545b090dfbb22a3f4923a3Alex Valavanis    return result;
f6d8c168c2d1b1e168545b090dfbb22a3f4923a3Alex Valavanis}
f6d8c168c2d1b1e168545b090dfbb22a3f4923a3Alex ValavanisSBasis operator*(SBasis const &a, double k);
f6d8c168c2d1b1e168545b090dfbb22a3f4923a3Alex Valavanisinline SBasis operator*(double k, SBasis const &a) { return a*k; }
f6d8c168c2d1b1e168545b090dfbb22a3f4923a3Alex Valavanisinline SBasis operator/(SBasis const &a, double k) { return a*(1./k); }
f6d8c168c2d1b1e168545b090dfbb22a3f4923a3Alex ValavanisSBasis& operator*=(SBasis& a, double b);
f6d8c168c2d1b1e168545b090dfbb22a3f4923a3Alex Valavanisinline SBasis& operator/=(SBasis& a, double b) { return (a*=(1./b)); }
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b//IMPL: AddableConcept
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_bSBasis operator+(const SBasis& a, const SBasis& b);
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_bSBasis operator-(const SBasis& a, const SBasis& b);
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_bSBasis& operator+=(SBasis& a, const SBasis& b);
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_bSBasis& operator-=(SBasis& a, const SBasis& b);
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b//TODO: remove?
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_binline SBasis operator+(const SBasis & a, Linear const & b) {
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b    if(b.isZero()) return a;
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b    if(a.isZero()) return b;
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b    SBasis result(a);
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b    result[0] += b;
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b    return result;
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b}
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_binline SBasis operator-(const SBasis & a, Linear const & b) {
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b    if(b.isZero()) return a;
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b    SBasis result(a);
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b    result[0] -= b;
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b    return result;
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b}
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_binline SBasis& operator+=(SBasis& a, const Linear& b) {
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b    if(a.isZero())
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b        a.push_back(b);
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b    else
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b        a[0] += b;
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b    return a;
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b}
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_binline SBasis& operator-=(SBasis& a, const Linear& b) {
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b    if(a.isZero())
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b        a.push_back(-b);
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b    else
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b        a[0] -= b;
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b    return a;
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b}
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b//IMPL: OffsetableConcept
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_binline SBasis operator+(const SBasis & a, double b) {
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b    if(a.isZero()) return Linear(b, b);
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b    SBasis result(a);
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b    result[0] += b;
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b    return result;
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b}
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_binline SBasis operator-(const SBasis & a, double b) {
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b    if(a.isZero()) return Linear(-b, -b);
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b    SBasis result(a);
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b    result[0] -= b;
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b    return result;
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b}
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_binline SBasis& operator+=(SBasis& a, double b) {
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b    if(a.isZero())
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b        a.push_back(Linear(b,b));
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b    else
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b        a[0] += b;
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b    return a;
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b}
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_binline SBasis& operator-=(SBasis& a, double b) {
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b    if(a.isZero())
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b        a.push_back(Linear(-b,-b));
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b    else
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b        a[0] -= b;
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b    return a;
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b}
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_bSBasis shift(SBasis const &a, int sh);
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_bSBasis shift(Linear const &a, int sh);
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_binline SBasis truncate(SBasis const &a, unsigned terms) {
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b    SBasis c;
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b    c.insert(c.begin(), a.begin(), a.begin() + std::min(terms, (unsigned)a.size()));
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b    return c;
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b}
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_bSBasis multiply(SBasis const &a, SBasis const &b);
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_bSBasis integral(SBasis const &c);
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_bSBasis derivative(SBasis const &a);
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_bSBasis sqrt(SBasis const &a, int k);
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b// return a kth order approx to 1/a)
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_bSBasis reciprocal(Linear const &a, int k);
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_bSBasis divide(SBasis const &a, SBasis const &b, int k);
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_binline SBasis operator*(SBasis const & a, SBasis const & b) {
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b    return multiply(a, b);
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b}
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_binline SBasis& operator*=(SBasis& a, SBasis const & b) {
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b    a = multiply(a, b);
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b    return a;
b465906c7b234e8f2413c12a1a92d1dd44077405Alex Valavanis}
d19de0ff8b0793e74e98429ded50868ce7e9f75cAlex Valavanis
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b//valuation: degree of the first non zero coefficient.
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_binline unsigned
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_bvaluation(SBasis const &a, double tol=0){
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b    unsigned val=0;
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b    while( val<a.size() &&
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b           fabs(a[val][0])<tol &&
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b           fabs(a[val][1])<tol )
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b        val++;
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b    return val;
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b}
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b// a(b(t))
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_bSBasis compose(SBasis const &a, SBasis const &b);
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_bSBasis compose(SBasis const &a, SBasis const &b, unsigned k);
d19de0ff8b0793e74e98429ded50868ce7e9f75cAlex ValavanisSBasis inverse(SBasis a, int k);
d19de0ff8b0793e74e98429ded50868ce7e9f75cAlex Valavanis//compose_inverse(f,g)=compose(f,inverse(g)), but is numerically more stable in some good cases...
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b//TODO: requires g(0)=0 & g(1)=1 atm. generalization should be obvious.
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_bSBasis compose_inverse(SBasis const &f, SBasis const &g, unsigned order=2, double tol=1e-3);
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b
d19de0ff8b0793e74e98429ded50868ce7e9f75cAlex Valavanisinline SBasis portion(const SBasis &t, double from, double to) { return compose(t, Linear(from, to)); }
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b
d19de0ff8b0793e74e98429ded50868ce7e9f75cAlex Valavanis// compute f(g)
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_binline SBasis
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_bSBasis::operator()(SBasis const & g) const {
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b    return compose(*this, g);
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b}
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_binline std::ostream &operator<< (std::ostream &out_file, const Linear &bo) {
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b    out_file << "{" << bo[0] << ", " << bo[1] << "}";
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b    return out_file;
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b}
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_binline std::ostream &operator<< (std::ostream &out_file, const SBasis & p) {
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b    for(unsigned i = 0; i < p.size(); i++) {
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b        out_file << p[i] << "s^" << i << " + ";
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b    }
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b    return out_file;
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b}
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_bSBasis sin(Linear bo, int k);
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_bSBasis cos(Linear bo, int k);
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_bstd::vector<double> roots(SBasis const & s);
d19de0ff8b0793e74e98429ded50868ce7e9f75cAlex Valavanisstd::vector<std::vector<double> > multi_roots(SBasis const &f,
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b                                 std::vector<double> const &levels,
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b                                 double htol=1e-7,
d19de0ff8b0793e74e98429ded50868ce7e9f75cAlex Valavanis                                 double vtol=1e-7,
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b                                 double a=0,
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b                                 double b=1);
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b}
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b/*
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b  Local Variables:
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b  mode:c++
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b  c-file-style:"stroustrup"
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b  c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b  indent-tabs-mode:nil
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b  fill-column:99
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b  End:
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b*/
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b#endif
a0df1b8dd5b14367c583ce2f72a2ca6bf1cde799gustav_b