numerics/recursion.h

// -*- c++ -*-

// +--------------------------------------------------------------------+
// | This file is part of MathTL - the Mathematical Template Library    |
// |                                                                    |
// | Copyright (c) 2002-2009                                            |
// | Thorsten Raasch, Manuel Werner                                     |
// +--------------------------------------------------------------------+

#ifndef _MATHTL_RECURSION_H
#define _MATHTL_RECURSION_H

#include <cmath>
#include <algebra/polynomial.h>
#include <algebra/vector.h>

namespace MathTL
{
  // some tools for three-term recursions
  // reference: Deuflhard/Bornemann, Numerische Mathematik 1

  class Recursion
  {
  public:
    virtual ~Recursion() {}

    virtual double a(const unsigned int k) const = 0;
    virtual double b(const unsigned int k) const = 0;
    virtual double c(const unsigned int k) const = 0;

    double operator () (const unsigned int n) const;

    double forwardSummation(const Vector<double>& coeffs) const;


    double adjointSummation(const Vector<double>& coeffs) const;
    
  protected:
    double p0_, p1_;
  };

  class HomogeneousRecursion : public Recursion
  {
  public:
    double c(const unsigned int k) const { return 0.0; }
  };

  class MonomialRecursion: public HomogeneousRecursion
  {
  public:
    MonomialRecursion(const double x) : x_(x) { p0_ = 1.0; p1_ = x_; }

    double a(const unsigned int k) const { return x_; }
    double b(const unsigned int k) const { return 0.0; }

  protected:
    double x_;
  };

  class ChebyshevRecursion : public HomogeneousRecursion
  {
  public:
    ChebyshevRecursion(const double x) : x_(x) { p0_ = 1.0; p1_ = x_; }

    double a(const unsigned int k) const { return 2*x_; }
    double b(const unsigned int k) const { return -1.0; }

  protected:
    double x_;
  };

  class CosineRecursion : public HomogeneousRecursion
  {
  public:
    CosineRecursion(const double x) : x_(x) { p0_ = 1.0; p1_ = cos(x_); }

    double a(const unsigned int k) const { return 2*cos(x_); }
    double b(const unsigned int k) const { return -1.0; }

  protected:
    double x_;
  };

  class SineRecursion : public HomogeneousRecursion
  {
  public:
    SineRecursion(const double x) : x_(x) { p0_ = 0.0; p1_ = sin(x_); }

    double a(const unsigned int k) const { return 2*cos(x_); }
    double b(const unsigned int k) const { return -1.0; }

  protected:
    double x_;
  };
  
  class JacobiRecursion : public HomogeneousRecursion
  {
  public:
    JacobiRecursion(const double a, const double b, const double x);
  
    double a(const unsigned int k) const;
    double b(const unsigned int k) const;
    
  protected:
    double a_, b_, x_;
  };

  class LegendreRecursion : public JacobiRecursion
  {
  public:
    LegendreRecursion(const double x) : JacobiRecursion(0.0,0.0,x) {}
  };
}

// include implementation of inline functions
#include <numerics/recursion.cpp>

#endif