FUnifRoots.hpp 4.68 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23
// ===================================================================================
// Copyright ScalFmm 2011 INRIA, Olivier Coulaud, Bérenger Bramas, Matthias Messner
// olivier.coulaud@inria.fr, berenger.bramas@inria.fr
// This software is a computer program whose purpose is to compute the FMM.
//
// This software is governed by the CeCILL-C and LGPL licenses and
// abiding by the rules of distribution of free software.
//
// 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 and CeCILL-C Licenses for more details.
// "http://www.cecill.info".
// "http://www.gnu.org/licenses".
// ===================================================================================
#ifndef FUNIFROOTS_HPP
#define FUNIFROOTS_HPP

#include <cmath>
#include <limits>
#include <cassert>

#include "../../Utils/FNoCopyable.hpp"
24
#include "../../Utils/FMath.hpp"
25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61


/**
 * @author Pierre Blanchard (pierre.blanchard@inria.fr)
 * Please read the license
 */

/**
 * @class FUnifRoots
 *
 * The class @p FUnifRoots provides the equispaced roots of order \f$\ell\f$
 * and the Lagrange polynomials \f$L_n(x)\f$.
 *
 * @tparam ORDER interpolation order \f$\ell\f$
 */
template <int ORDER>
struct FUnifRoots : FNoCopyable
{
  enum {order = ORDER}; //!< interpolation order

  /**
   * Lagrange roots in [-1,1] computed as \f$\bar x_n =
   * -1 + 2\frac{n-1}{\ell}\f$ for \f$n=1,\dots,\ell\f$
   */
  const static double roots[];

  /**
   * Lagrange polynomials \f$ L_n(x) = \Pi_{m=1 \atop m\neq n}^\ell \frac{x-\bar x_m}{\bar x_n-\bar x_m} \f$
   *
   * @param[in] n index
   * @param[in] x coordinate in [-1,1]
   * @return function value
   */
  static FReal L(const unsigned int n, FReal x)
  {
    assert(std::fabs(x)-1.<10.*std::numeric_limits<FReal>::epsilon());
    if (std::fabs(x)>1.) {
62
      //std::cout << "x=" << x << " out of bounds!" << std::endl;
63 64 65 66
      x = (x > FReal( 1.) ? FReal( 1.) : x);
      x = (x < FReal(-1.) ? FReal(-1.) : x);
    }

67 68 69 70 71 72 73 74 75 76
    // Specific precomputation of scale factor 
    // in order to minimize round-off errors
    // NB: scale factor could be hardcoded (just as the roots)
    FReal scale;
    int omn = order-n-1;
    if(omn%2) scale=-1.; // (-1)^(n-1-(k+1)+1)=(-1)^(omn-1)
    else scale=1.;
    scale*=FMath::pow(2.,order-1)*FMath::factorial(n)*FMath::factorial(omn);

    // compute L
77
    FReal L=FReal(1.);
78
    for(unsigned int m=0;m<order;++m){
79 80 81 82 83 84 85 86
      if(m!=n){
        // previous version with risks of round-off error
        //L *= (x-FUnifRoots<order>::roots[m])/(FUnifRoots<order>::roots[n]-FUnifRoots<order>::roots[m]);

        // new version (reducing round-off)
        // regular grid on [-1,1] (h simplifies, only the size of the domain and a remains i.e. 2. and -1.)
        L *= ((order-1)*(x+1.)-2.*m); 
      }
87
    }
88 89
    
    L*=scale;
90

91
    return FReal(L);
92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110
  }


  /**
   * For the derivation of the Lagrange polynomials
   * \f$ \frac{\mathrm{d} L_n(x)}{\mathrm{d}x} = ... \f$
   *
   * @param[in] n index
   * @param[in] x coordinate in [-1,1]
   * @return function value
   */
  static FReal dL(const unsigned int n, FReal x)
  {
    assert(std::fabs(x)-1.<10.*std::numeric_limits<FReal>::epsilon());
    if (std::fabs(x)>1.) {
      x = (x > FReal( 1.) ? FReal( 1.) : x);
      x = (x < FReal(-1.) ? FReal(-1.) : x);
    }

111 112 113 114
    // optimized variant
    FReal NdL=FReal(0.);// init numerator
    FReal DdL=FReal(1.);// init denominator
    FReal tmpNdL;
115 116
    for(unsigned int p=0;p<order;++p){
      if(p!=n){
117 118 119 120 121 122
        tmpNdL=FReal(1.);
        for(unsigned int m=0;m<order;++m)
          if(m!=n && m!=p)
            tmpNdL*=x-FUnifRoots<order>::roots[m];
        NdL+=tmpNdL;
        DdL*=FUnifRoots<order>::roots[n]-FUnifRoots<order>::roots[p];
123 124 125
      }//endif
    }// p

126
    return FReal(NdL/DdL);
127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173

  }
};

// We declare the roots here only once Please look to .cpp for definitions

// order 2
template<> const double FUnifRoots<2>::roots[];

// order 3
template<> const double FUnifRoots<3>::roots[];

// order 4
template<> const double FUnifRoots<4>::roots[];

// order 5
template<> const double FUnifRoots<5>::roots[];

// order 6
template<> const double FUnifRoots<6>::roots[];

// order 7
template<> const double FUnifRoots<7>::roots[];

// order 8
template<> const double FUnifRoots<8>::roots[];

// order 9
template<> const double FUnifRoots<9>::roots[];

// order 10
template<> const double FUnifRoots<10>::roots[];

// order 11
template<> const double FUnifRoots<11>::roots[];

// order 12
template<> const double FUnifRoots<12>::roots[];

// order 13
template<> const double FUnifRoots<13>::roots[];

// order 14
template<> const double FUnifRoots<14>::roots[];


#endif