sbasis.h revision e6bdf746e2d9e775704a475a29cc1bb167ec271c
a63dbee6633e3ed991bb2b34cbed7c9f02c1a839Jon A. Cruz * sbasis.h - S-power basis function class
a63dbee6633e3ed991bb2b34cbed7c9f02c1a839Jon A. Cruz * Nathan Hurst <njh@mail.csse.monash.edu.au>
a63dbee6633e3ed991bb2b34cbed7c9f02c1a839Jon A. Cruz * Michael Sloan <mgsloan@gmail.com>
a63dbee6633e3ed991bb2b34cbed7c9f02c1a839Jon A. Cruz * Copyright (C) 2006-2007 authors
a63dbee6633e3ed991bb2b34cbed7c9f02c1a839Jon A. Cruz * This library is free software; you can redistribute it and/or
00339f176618a0eebcbb2261306a9cd841ad43ecJon A. Cruz * modify it either under the terms of the GNU Lesser General Public
a63dbee6633e3ed991bb2b34cbed7c9f02c1a839Jon A. Cruz * License version 2.1 as published by the Free Software Foundation
a63dbee6633e3ed991bb2b34cbed7c9f02c1a839Jon A. Cruz * (the "LGPL") or, at your option, under the terms of the Mozilla
a63dbee6633e3ed991bb2b34cbed7c9f02c1a839Jon A. Cruz * Public License Version 1.1 (the "MPL"). If you do not alter this
a63dbee6633e3ed991bb2b34cbed7c9f02c1a839Jon A. Cruz * notice, a recipient may use your version of this file under either
a63dbee6633e3ed991bb2b34cbed7c9f02c1a839Jon A. Cruz * the MPL or the LGPL.
a63dbee6633e3ed991bb2b34cbed7c9f02c1a839Jon A. Cruz * You should have received a copy of the LGPL along with this library
a63dbee6633e3ed991bb2b34cbed7c9f02c1a839Jon A. Cruz * in the file COPYING-LGPL-2.1; if not, write to the Free Software
a63dbee6633e3ed991bb2b34cbed7c9f02c1a839Jon A. Cruz * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
a63dbee6633e3ed991bb2b34cbed7c9f02c1a839Jon A. Cruz * You should have received a copy of the MPL along with this library
a63dbee6633e3ed991bb2b34cbed7c9f02c1a839Jon A. Cruz * in the file COPYING-MPL-1.1
a63dbee6633e3ed991bb2b34cbed7c9f02c1a839Jon A. Cruz * The contents of this file are subject to the Mozilla Public License
47badd0035ae8c9135c51444f3770b17ae504ddaAlex Valavanis * Version 1.1 (the "License"); you may not use this file except in
47badd0035ae8c9135c51444f3770b17ae504ddaAlex Valavanis * compliance with the License. You may obtain a copy of the License at
a63dbee6633e3ed991bb2b34cbed7c9f02c1a839Jon A. Cruz * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
a63dbee6633e3ed991bb2b34cbed7c9f02c1a839Jon A. Cruz * OF ANY KIND, either express or implied. See the LGPL or the MPL for
17238dd97c572210f1036dd8a7a8657e67bfa9c3Campbell Barton * the specific language governing rights and limitations.
2c312e5f8adf7306ce4b41fb6639475df9c3ea87Jon A. Cruz/*** An empty SBasis is identically 0. */
c864652e7fca4ed7605835d5cb6661dc422cabaaJon A. Cruz //IMPL: FragmentConcept
c864652e7fca4ed7605835d5cb6661dc422cabaaJon A. Cruz typedef double output_type;
c864652e7fca4ed7605835d5cb6661dc422cabaaJon A. Cruz inline bool isZero() const {
c864652e7fca4ed7605835d5cb6661dc422cabaaJon A. Cruz if(empty()) return true;
c864652e7fca4ed7605835d5cb6661dc422cabaaJon A. Cruz for(unsigned i = 0; i < size(); i++) {
c864652e7fca4ed7605835d5cb6661dc422cabaaJon A. Cruz return true;
c864652e7fca4ed7605835d5cb6661dc422cabaaJon A. Cruz bool isFinite() const;
c864652e7fca4ed7605835d5cb6661dc422cabaaJon A. Cruz inline double at0() const {
c864652e7fca4ed7605835d5cb6661dc422cabaaJon A. Cruz if(empty()) return 0; else return (*this)[0][0];
c864652e7fca4ed7605835d5cb6661dc422cabaaJon A. Cruz inline double at1() const{
c864652e7fca4ed7605835d5cb6661dc422cabaaJon A. Cruz if(empty()) return 0; else return (*this)[0][1];
c864652e7fca4ed7605835d5cb6661dc422cabaaJon A. Cruz double valueAt(double t) const {
c864652e7fca4ed7605835d5cb6661dc422cabaaJon A. Cruz double s = t*(1-t);
c864652e7fca4ed7605835d5cb6661dc422cabaaJon A. Cruz//TODO: rewrite as horner
c864652e7fca4ed7605835d5cb6661dc422cabaaJon A. Cruz for(unsigned k = 0; k < size(); k++) {
c864652e7fca4ed7605835d5cb6661dc422cabaaJon A. Cruz double valueAndDerivative(double t, double &der) const {
c864652e7fca4ed7605835d5cb6661dc422cabaaJon A. Cruz double s = t*(1-t);
c864652e7fca4ed7605835d5cb6661dc422cabaaJon A. Cruz//TODO: rewrite as horner
c864652e7fca4ed7605835d5cb6661dc422cabaaJon A. Cruz for(unsigned k = 0; k < size(); k++) {
c864652e7fca4ed7605835d5cb6661dc422cabaaJon A. Cruz // p0 and p1 at this point form a linear approximation at t
c864652e7fca4ed7605835d5cb6661dc422cabaaJon A. Cruz double operator()(double t) const {
a63dbee6633e3ed991bb2b34cbed7c9f02c1a839Jon A. Cruz std::vector<double> valueAndDerivatives(double /*t*/, unsigned /*n*/) const {
a63dbee6633e3ed991bb2b34cbed7c9f02c1a839Jon A. Cruz SBasis toSBasis() const { return SBasis(*this); }
a63dbee6633e3ed991bb2b34cbed7c9f02c1a839Jon A. Cruz// compute f(g)
a63dbee6633e3ed991bb2b34cbed7c9f02c1a839Jon A. Cruz//MUTATOR PRISON
a63dbee6633e3ed991bb2b34cbed7c9f02c1a839Jon A. Cruz Linear& operator[](unsigned i) { return this->at(i); }
a63dbee6633e3ed991bb2b34cbed7c9f02c1a839Jon A. Cruz //remove extra zeros
a63dbee6633e3ed991bb2b34cbed7c9f02c1a839Jon A. Cruz while(!empty() && 0 == back()[0] && 0 == back()[1])
a63dbee6633e3ed991bb2b34cbed7c9f02c1a839Jon A. Cruz void truncate(unsigned k) { if(k < size()) resize(k); }
a63dbee6633e3ed991bb2b34cbed7c9f02c1a839Jon A. Cruz//TODO: figure out how to stick this in linear, while not adding an sbasis dep
a63dbee6633e3ed991bb2b34cbed7c9f02c1a839Jon A. Cruzinline SBasis Linear::toSBasis() const { return SBasis(*this); }
a63dbee6633e3ed991bb2b34cbed7c9f02c1a839Jon A. Cruz//implemented in sbasis-roots.cpp
a63dbee6633e3ed991bb2b34cbed7c9f02c1a839Jon A. CruzInterval bounds_fast(SBasis const &a, int order = 0);
a63dbee6633e3ed991bb2b34cbed7c9f02c1a839Jon A. CruzInterval bounds_local(SBasis const &a, const Interval &t, int order = 0);
a63dbee6633e3ed991bb2b34cbed7c9f02c1a839Jon A. Cruz for(unsigned k = 0; k < a.size(); k++)
a63dbee6633e3ed991bb2b34cbed7c9f02c1a839Jon A. Cruz//IMPL: ScalableConcept
a63dbee6633e3ed991bb2b34cbed7c9f02c1a839Jon A. Cruz for(unsigned i = 0; i < p.size(); i++) {
a63dbee6633e3ed991bb2b34cbed7c9f02c1a839Jon A. Cruzinline SBasis operator*(double k, SBasis const &a) { return a*k; }
a63dbee6633e3ed991bb2b34cbed7c9f02c1a839Jon A. Cruzinline SBasis operator/(SBasis const &a, double k) { return a*(1./k); }
a63dbee6633e3ed991bb2b34cbed7c9f02c1a839Jon A. Cruzinline SBasis& operator/=(SBasis& a, double b) { return (a*=(1./b)); }
a63dbee6633e3ed991bb2b34cbed7c9f02c1a839Jon A. Cruz//IMPL: AddableConcept
a63dbee6633e3ed991bb2b34cbed7c9f02c1a839Jon A. CruzSBasis operator+(const SBasis& a, const SBasis& b);
a63dbee6633e3ed991bb2b34cbed7c9f02c1a839Jon A. CruzSBasis operator-(const SBasis& a, const SBasis& b);
a63dbee6633e3ed991bb2b34cbed7c9f02c1a839Jon A. Cruz//TODO: remove?
a63dbee6633e3ed991bb2b34cbed7c9f02c1a839Jon A. Cruzinline SBasis operator+(const SBasis & a, Linear const & b) {
a63dbee6633e3ed991bb2b34cbed7c9f02c1a839Jon A. Cruz if(b.isZero()) return a;
a63dbee6633e3ed991bb2b34cbed7c9f02c1a839Jon A. Cruz if(a.isZero()) return b;
a63dbee6633e3ed991bb2b34cbed7c9f02c1a839Jon A. Cruzinline SBasis operator-(const SBasis & a, Linear const & b) {
a63dbee6633e3ed991bb2b34cbed7c9f02c1a839Jon A. Cruz if(b.isZero()) return a;
a63dbee6633e3ed991bb2b34cbed7c9f02c1a839Jon A. Cruzinline SBasis& operator+=(SBasis& a, const Linear& b) {
a63dbee6633e3ed991bb2b34cbed7c9f02c1a839Jon A. Cruzinline SBasis& operator-=(SBasis& a, const Linear& b) {
a63dbee6633e3ed991bb2b34cbed7c9f02c1a839Jon A. Cruz//IMPL: OffsetableConcept
a63dbee6633e3ed991bb2b34cbed7c9f02c1a839Jon A. Cruzinline SBasis operator+(const SBasis & a, double b) {
a63dbee6633e3ed991bb2b34cbed7c9f02c1a839Jon A. Cruzinline SBasis operator-(const SBasis & a, double b) {
a63dbee6633e3ed991bb2b34cbed7c9f02c1a839Jon A. Cruzinline SBasis truncate(SBasis const &a, unsigned terms) {
a63dbee6633e3ed991bb2b34cbed7c9f02c1a839Jon A. Cruz c.insert(c.begin(), a.begin(), a.begin() + std::min(terms, (unsigned)a.size()));
a63dbee6633e3ed991bb2b34cbed7c9f02c1a839Jon A. CruzSBasis multiply(SBasis const &a, SBasis const &b);
a63dbee6633e3ed991bb2b34cbed7c9f02c1a839Jon A. Cruz// return a kth order approx to 1/a)
a63dbee6633e3ed991bb2b34cbed7c9f02c1a839Jon A. CruzSBasis divide(SBasis const &a, SBasis const &b, int k);
2c312e5f8adf7306ce4b41fb6639475df9c3ea87Jon A. Cruzinline SBasis operator*(SBasis const & a, SBasis const & b) {
00339f176618a0eebcbb2261306a9cd841ad43ecJon A. Cruz return multiply(a, b);
00339f176618a0eebcbb2261306a9cd841ad43ecJon A. Cruzinline SBasis& operator*=(SBasis& a, SBasis const & b) {
00339f176618a0eebcbb2261306a9cd841ad43ecJon A. Cruz//valuation: degree of the first non zero coefficient.
a63dbee6633e3ed991bb2b34cbed7c9f02c1a839Jon A. Cruzinline unsigned
00339f176618a0eebcbb2261306a9cd841ad43ecJon A. CruzSBasis compose(SBasis const &a, SBasis const &b);
00339f176618a0eebcbb2261306a9cd841ad43ecJon A. CruzSBasis compose(SBasis const &a, SBasis const &b, unsigned k);
a63dbee6633e3ed991bb2b34cbed7c9f02c1a839Jon A. Cruz//compose_inverse(f,g)=compose(f,inverse(g)), but is numerically more stable in some good cases...
a63dbee6633e3ed991bb2b34cbed7c9f02c1a839Jon A. Cruz//TODO: requires g(0)=0 & g(1)=1 atm. generalization should be obvious.
a63dbee6633e3ed991bb2b34cbed7c9f02c1a839Jon A. CruzSBasis compose_inverse(SBasis const &f, SBasis const &g, unsigned order=2, double tol=1e-3);
a63dbee6633e3ed991bb2b34cbed7c9f02c1a839Jon A. Cruzinline SBasis portion(const SBasis &t, double from, double to) { return compose(t, Linear(from, to)); }
a63dbee6633e3ed991bb2b34cbed7c9f02c1a839Jon A. Cruz// compute f(g)
a63dbee6633e3ed991bb2b34cbed7c9f02c1a839Jon A. Cruzinline std::ostream &operator<< (std::ostream &out_file, const Linear &bo) {
a63dbee6633e3ed991bb2b34cbed7c9f02c1a839Jon A. Cruz out_file << "{" << bo[0] << ", " << bo[1] << "}";
a63dbee6633e3ed991bb2b34cbed7c9f02c1a839Jon A. Cruzinline std::ostream &operator<< (std::ostream &out_file, const SBasis & p) {
a63dbee6633e3ed991bb2b34cbed7c9f02c1a839Jon A. Cruz for(unsigned i = 0; i < p.size(); i++) {
a63dbee6633e3ed991bb2b34cbed7c9f02c1a839Jon A. Cruzstd::vector<std::vector<double> > multi_roots(SBasis const &f,
a63dbee6633e3ed991bb2b34cbed7c9f02c1a839Jon A. Cruz double b=1);
a63dbee6633e3ed991bb2b34cbed7c9f02c1a839Jon A. Cruz Local Variables:
a63dbee6633e3ed991bb2b34cbed7c9f02c1a839Jon A. Cruz c-file-style:"stroustrup"
a63dbee6633e3ed991bb2b34cbed7c9f02c1a839Jon A. Cruz c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
a63dbee6633e3ed991bb2b34cbed7c9f02c1a839Jon A. Cruz indent-tabs-mode:nil
a63dbee6633e3ed991bb2b34cbed7c9f02c1a839Jon A. Cruz fill-column:99
a63dbee6633e3ed991bb2b34cbed7c9f02c1a839Jon A. Cruz// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :