poly.hpp

Go to the documentation of this file.

/*
 * $Id: poly.hpp 309 2013-11-17 16:26:16Z jdl3 $
 * Copyright (C) 2012 John D Lamb
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 * 
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */
 
#ifndef CCGSL_POLY_HPP
#define CCGSL_POLY_HPP
 
#include<stdexcept>
#include<gsl/gsl_poly.h>
#include<gsl/gsl_errno.h>
#include"complex.hpp"
 
namespace gsl {
  namespace poly {
    class complex_workspace {
    public:
      complex_workspace(){
    ccgsl_pointer = 0;
    count = 0; // initially nullptr will do
      }
      // Refines random access container
      // Refines assignable
      explicit complex_workspace( size_t const n ){
    ccgsl_pointer = gsl_poly_complex_workspace_alloc( n );
    // just plausibly we could allocate complex_workspace but not count
    try { count = new size_t; } catch( std::bad_alloc& e ){
      // try to tidy up before rethrowing
      gsl_poly_complex_workspace_free( ccgsl_pointer );
      throw e;
    }
    *count = 1; // initially there is just one reference to ccgsl_pointer
      }
      explicit complex_workspace( gsl_poly_complex_workspace* v ){
    ccgsl_pointer = v;
    // just plausibly we could fail to allocate count: no further action needed.
    count = new size_t;
    *count = 1; // initially there is just one reference to ccgsl_pointer
      }
      // copy constructor
      complex_workspace( complex_workspace const& v ){ ccgsl_pointer = v.ccgsl_pointer;
    count = v.count; if( count != 0 ) ++*count; }
      // assignment operator
      complex_workspace& operator=( complex_workspace const& v ){
    // first, possibly delete anything pointed to by this
    if( count == 0 or --*count == 0 ){
      if( ccgsl_pointer != 0 ) gsl_poly_complex_workspace_free( ccgsl_pointer );
      delete count;
    } // Then copy
    ccgsl_pointer = v.ccgsl_pointer; count = v.count; if( count != 0 ) ++*count; return *this;
      }
      // destructor
      ~complex_workspace(){
    if( count == 0 or --*count == 0 ){
      // could have allocated null pointer
      if( ccgsl_pointer != 0 ) gsl_poly_complex_workspace_free( ccgsl_pointer );
      delete count;
    }
      }
#ifdef __GXX_EXPERIMENTAL_CXX0X__
    complex_workspace( complex_workspace&& v ) : ccgsl_pointer( v.ccgsl_pointer ), count( nullptr ){
      std::swap( count, v.count );
      v.ccgsl_pointer = nullptr;
    }
    complex_workspace& operator=( complex_workspace&& v ){
      complex_workspace( std::move( v ) ).swap( *this );
      return *this;
    }
#endif
      // Refines equality comparable
      // == operator
      bool operator==( complex_workspace const& v ) const { return ccgsl_pointer == v.ccgsl_pointer; }
      // != operator
      bool operator!=( complex_workspace const& v ) const { return not operator==( v ); }
      // Refines forward container
      // Refines less than comparable
      // operator<
      bool operator<( complex_workspace const& v ) const { return ccgsl_pointer < v.ccgsl_pointer; }
      // operator>
      bool operator>( complex_workspace const& v ) const { return ccgsl_pointer > v.ccgsl_pointer; }
      // operator<=
      bool operator<=( complex_workspace const& v ) const { return ccgsl_pointer <= v.ccgsl_pointer; }
      // operator>=
      bool operator>=( complex_workspace const& v ) const { return ccgsl_pointer >= v.ccgsl_pointer; }
      bool empty() const { return ccgsl_pointer == 0; }
      // swap() --- should work even if sizes don't match
      void swap( complex_workspace& v ){
    std::swap( ccgsl_pointer, v.ccgsl_pointer );
    std::swap( count, v.count );
      }
    private:
      gsl_poly_complex_workspace* ccgsl_pointer;
      size_t* count;
    public:
      // shared reference functions
      gsl_poly_complex_workspace* get() const { return ccgsl_pointer; }
      bool unique() const { return count != 0 and *count == 1; }
      size_t use_count() const { return count == 0 ? 0 : *count; }
#ifdef __GXX_EXPERIMENTAL_CXX0X__
      explicit
#endif
      operator bool() const { return ccgsl_pointer != 0; }
    };
 
#ifndef DOXYGEN_SKIP
    inline double eval( double const c[], int const len, double const x ){
      return gsl_poly_eval( c, len, x ); }
#endif // DOXYGEN_SKIP
 
    template<typename T>
    inline double eval( T const& c, double const x ){
      return gsl_poly_eval( c.data(), c.size(), x ); }
 
#ifndef DOXYGEN_SKIP
    inline complex complex_eval( double const c[], int const len, gsl::complex const z ){
      return gsl_poly_complex_eval( c, len, z.get() ); }
#endif // DOXYGEN_SKIP
    
    template<typename T>
    inline complex complex_eval( T const& c, gsl::complex const z ){
      return gsl_poly_complex_eval( c.data(), c.size() , z.get() ); }
 
#ifndef DOXYGEN_SKIP
    inline int eval_derivs( double const c[], size_t const lenc, double const x,
              double res[], size_t const lenres ){
      return gsl_poly_eval_derivs( c, lenc, x, res, lenres ); }
#endif // DOXYGEN_SKIP
    
    template<typename T, typename U>
    inline int eval_derivs( T const& c, double const x, U& res ){
      return gsl_poly_eval_derivs( c.data(), c.size(), x, res.data(), res.size() ); }
 
#ifndef DOXYGEN_SKIP
    inline int dd_init( double dd[], double const x[], double const y[], size_t size ){
      return gsl_poly_dd_init( dd, x, y, size ); }
#endif // DOXYGEN_SKIP
    
    template<typename T, typename U, typename V>
    inline int dd_init( T& dd, U const& x, V const& y ){
      size_t const size = dd.size() < x.size() ?
    ( dd.size() < y.size() ? dd.size() : y.size() ) : x.size();
      return gsl_poly_dd_init( dd.data(), x.data(), y.data(), size ); }
    
#ifndef DOXYGEN_SKIP
    inline int
    dd_hermite_init( double dd[], double z[], double const xa[], double const ya[],
             double const dya[], size_t const size ){
      return gsl_poly_dd_hermite_init( dd, z, xa, ya, dya, size ); }
#endif // DOXYGEN_SKIP
 
    template<typename DD, typename Z, typename XA, typename YA, typename DYA>
    inline int
    dd_hermite_init( DD& dd, Z& z, XA const& xa, YA const& ya,
             DYA const& dya){
      size_t const size {xa.size()};
      if(z.size() != 2*size)
    throw std::length_error("gsl::poly::dd_hermite_init(): "
                "z must have twice the length of xa.");
      if(dd.size() != size)
    throw std::length_error("gsl::poly::dd_hermite_init(): xa and dd must have same length.");
      if(ya.size() != size)
    throw std::length_error("gsl::poly::dd_hermite_init(): xa and ya must have same length.");
      if(dya.size() != size)
    throw std::length_error("gsl::poly::dd_hermite_init(): "
                "xa and dya must have same length.");
      return gsl_poly_dd_hermite_init( dd.data(), z.data(), xa.data(), ya.data(), dya.data(),
                       size );
    }
 
 
#ifndef DOXYGEN_SKIP
    inline double dd_eval( double const dd[], double const xa[], size_t const size,
                double const x ){
      return gsl_poly_dd_eval( dd, xa, size, x ); }
#endif // DOXYGEN_SKIP
 
    template<typename T, typename U>
    inline double dd_eval( T const& dd, U const& xa, double const x ){
      size_t const size = dd.size() < xa.size() ? dd.size() : xa.size();
      return gsl_poly_dd_eval( dd.data(), xa.data(), size, x ); }
 
#ifndef DOXYGEN_SKIP
    inline int dd_taylor( double c[], double xp, double const dd[], double const x[],
            size_t size, double w[] ){
      return gsl_poly_dd_taylor( c, xp, dd, x, size, w ); }
#endif // DOXYGEN_SKIP
 
    template<typename C, typename DD, typename X, typename W>
    inline int dd_taylor( C& c, double xp, DD const& dd, X const& x, W& w ){
      size_t const size = c.size() < dd.size() ?
    ( c.size() < x.size() ? c.size() : x.size() ) : dd.size();
      if( w.size() < size )
    gsl_error( "workspace too small", __FILE__, __LINE__, GSL_EBADLEN );
      return gsl_poly_dd_taylor( c.data(), xp, dd.data(), x.data(), size, w.data() ); }
 
#ifndef DOXYGEN_SKIP
    inline int solve_quadratic( double a, double b, double c, double* x0, double* x1 ){
      return gsl_poly_solve_quadratic( a, b, c, x0, x1 ); }
#endif // DOXYGEN_SKIP
    inline int solve_quadratic( double a, double b, double c, double& x0, double& x1 ){
      return gsl_poly_solve_quadratic( a, b, c, &x0, &x1 ); }
 
    inline int complex_solve_quadratic( double a, double b, double c,
                    complex& z0, complex& z1 ){
      return gsl_poly_complex_solve_quadratic( a, b, c, &z0.get(), &z1.get() ); }
 
 
#ifndef DOXYGEN_SKIP
    inline int solve_cubic( double a, double b, double c, double* x0,
                 double* x1, double* x2 ){
      return gsl_poly_solve_cubic( a, b, c, x0, x1, x2 ); }
#endif // DOXYGEN_SKIP
    inline int solve_cubic( double a, double b, double c, double& x0,
                 double& x1, double& x2 ){
      return gsl_poly_solve_cubic( a, b, c, &x0, &x1, &x2 ); }
 
    inline int complex_solve_cubic( double a, double b, double c, complex& z0,
                    complex& z1, complex& z2 ){
      return gsl_poly_complex_solve_cubic( a, b, c, &z0.get(), &z1.get(), &z2.get() ); }
 
#ifndef DOXYGEN_SKIP
    inline int complex_solve( double const* a, size_t n, complex_workspace& w,
                  complex_packed_ptr z ){
      return gsl_poly_complex_solve( a, n, w.get(), z ); }
#endif // DOXYGEN_SKIP
 
    template<typename A, typename Z>
    inline int complex_solve( A const& a, complex_workspace& w, Z& z ){
      size_t const n = a.size();
      size_t const n_2 = z.size();
      if( 2 * n != n_2 + 2 )
    gsl_error( "mismatch in sizes of coefficients and roots",
           __FILE__, __LINE__, GSL_EBADLEN );
      return gsl_poly_complex_solve( a.data(), n, w.get(), z.data() ); }
  }
#ifdef DOXYGEN_SKIP
  namespace poly::complex_poly {
#else
  namespace complex_poly {
#endif
    inline complex complex_eval( gsl_complex c[], size_t const len, gsl::complex const z ){
      return gsl_complex_poly_complex_eval( c, len, z.get() ); }
  }
}
#endif