#ifndef FILE_RECURSIVE_POL
#define FILE_RECURSIVE_POL

/*********************************************************************/
/* File:   recursive_pol.hpp                                         */
/* Author: Start                                                     */
/* Date:   6. Feb. 2003                                              */
/*********************************************************************/

/*
  Recursive Polynomials
*/


inline double Mult (double x, double y) { return x*y; }
inline double Minus (double x, double y) { return x-y; }



/**
   Legendre polynomials P_i on (-1,1), 
   up to order n --> n+1 values
   
   l P_l = (2 l - 1) x P_{l-1} - (l-1) P_{l-2}

   P_0 = 1
   P_1 = x
   P_2 = -1/2 + 3/2 x^2
 */
template <class S, class T>
inline void LegendrePolynomial (int n, S x, T & values)
{
  S p1 = 1.0, p2 = 0.0, p3;

  if (n >= 0)
    values[0] = 1.0;

  for (int j=1; j<=n; j++)
    {
      p3 = p2; p2 = p1;
      p1 = ((2.0*j-1.0)*x*p2 - (j-1.0)*p3) / j;
      //      p1 = Minus ( ((2.0*j-1.0)/j) * Mult (x,p2) ,
      //		   ((j-1.0)/j) * p3);
      values[j] = p1;
    }
}


template <class S, class T>
inline void GegenbauerPolynomial (int n, S x, double lam, T & values)
{
  S p1 = 1.0, p2 = 0.0, p3;

  if (n >= 0)
    values[0] = 1.0;

  for (int j=1; j<=n; j++)
    {
      p3=p2; p2=p1;
      p1=( 2.0*(j+lam-1.0) * x * p2 - (j+2*lam-2.0) * p3) / j;
      values[j] = p1;
    }  
}


/**
   Integrated Legendre polynomials on (-1,1)

   value[0] = -1
   value[1] = x
   value[i] (x) = \int_{-1}^x P_{i-1} (s) ds  for i >= 2

   WARNING: is not \int P_i
 */
template <class S, class T>
inline void IntegratedLegendrePolynomial (int n, S x, T & values)
{
  S p1 = -1.0, p2 = 0.0, p3;

  if (n >= 0)
    values[0] = -1.0;

  for (int j=1; j<=n; j++)
    {
      p3=p2; p2=p1;
      p1=( (2*j-3) * x * p2 - (j-3) * p3) / j;
      values[j] = p1;
    }
}


/**
   Derived Legendre polynomials on (-1,1)
   value[i] (x) = d/dx P_{i+1}
 */
template <class S, class T>
inline void DerivedLegendrePolynomial (int n, S x, T & values)
{
  GegenbauerPolynomial<S,T> (n, x, 1.5, values);
}




/**
   Compute triangle edge-shape functions

   functions vanish on upper two edges

   x,y: coordinates in triangle (-1, 0), (1, 0), (0, 1)

   f_i (x, 0) = IntegratedLegendrePol_i (x)

   f_i ... pol of order i


   Monomial extension:
*/
template <class Sx, class Sy, class T>
inline void TriangleExtensionMonomial (int n, Sx x, Sy y, T & values)
{
  Sx p1 = -1.0, p2 = 0.0, p3;
  values[0] = -1.0;
  Sy fy = (1-y)*(1-y);
  for (int j=1; j<=n; j++)
    {
      p3=p2; p2=p1;
      p1=( (2*j-3) * x * p2 - (j-3) * fy * p3) / j;
      values[j] = p1;
    }    
}

template <class Sx, class Sy, class T>
inline void DiffTriangleExtensionMonomial (int n, Sx x, Sy y, T & values)
{
  ARRAY<AutoDiff<2> > ad_values(n+1);
  AutoDiff<2> ad_x(x, 0);
  AutoDiff<2> ad_y(y, 1);

  TriangleExtensionMonomial (n, ad_x, ad_y, ad_values);

  for (int i = 0; i <= n; i++)
    for (int j = 0; j < 2; j++)
      values(i,j) = ad_values[i].DValue(j);
}



/**
   Extension is the optimal averaging extension:
*/
template <class Sx, class Sy, class T>
inline void TriangleExtensionJacobi (int n, Sx x, Sy y, T & values)
{
  if ( (1-y) != 0.0)
    {
      int j;

      JacobiPolynomial (n-2, x / (1-y), 2, 2, values);
      Sy fac = (1.-x-y) * (1.+x-y);
      for (j = 0; j <= n-2; j++)
	{
	  values[j] *= fac;
	  fac *= 1-y;
	}
      for (j = n; j >= 2; j--)
	values[j] = values[j-2];
      if (n >= 0) values[0] = 0;
      if (n >= 1) values[1] = 0;
    }
  else
    {
      for (int j = 0; j <= n; j++)
	values[j] = 0;
    }
}

template <class Sx, class Sy, class T>
inline void DiffTriangleExtensionJacobi (int n, Sx x, Sy y, T & values)
{
  ARRAY<AutoDiff<2> > ad_values(n+1);
  AutoDiff<2> ad_x(x, 0);
  AutoDiff<2> ad_y(y, 1);

  TriangleExtensionJacobi (n, ad_x, ad_y, ad_values);
  for (int i = 0; i <= n; i++)
    for (int j = 0; j < 2; j++)
      values(i,j) = ad_values[i].DValue(j);
}







/**
   Extension is the optimal averaging extension:
*/
template <class Sx, class Sy, class T>
inline void TriangleExtensionOpt (int n, Sx x, Sy y, T & values)
{
  if (y < 1e-10)
    {
      IntegratedLegendrePolynomial (n, x, values);
    }
  else
    {
      ARRAY<Sx> ge1(n+2);
      ARRAY<Sx> ge2(n+2);
      ARRAY<Sx> ge3(n+2);
      ARRAY<Sx> ge4(n+2);

      GegenbauerPolynomial (n+1, Sx(-1.0), -1.5, ge1);
      GegenbauerPolynomial (n+1, x-y, -1.5, ge2);
      GegenbauerPolynomial (n+1, x+y, -1.5, ge3);
      GegenbauerPolynomial (n+1, Sx(1.0), -1.5, ge4);
 
      for (int i = 0; i <= n; i++)
	values[i] = 1.0/3.0 *
	  (  (2*y/(1+x+y)/(1+x+y) - y/2) * ge1[i+1]  +
	     (-1/(2*y) + 2*y/(1-x+y)/(1-x+y)) * ge2[i+1] +
	     (1/(2*y) - 2*y/(1+x+y)/(1+x+y) ) * ge3[i+1] +
	     (-2*y/(1-x+y)/(1-x+y) + y/2 ) * ge4[i+1] );
    }
}

template <class Sx, class Sy, class T>
inline void DiffTriangleExtensionOpt (int n, Sx x, Sy y, T & values)
{
  ARRAY<AutoDiff<2> > ad_values(n+1);
  AutoDiff<2> ad_x(x, 0);
  AutoDiff<2> ad_y(y, 1);

   if (y < 1e-10)
   {
       values = 0.;
   }
   else
   {
       TriangleExtensionOpt (n, ad_x, ad_y, ad_values);

       for (int i = 0; i <= n; i++)
         for (int j = 0; j < 2; j++)
           values(i,j) = ad_values[i].DValue(j);
  }
}

/* 
   E_i(x,y) = P_i(x/t) * t^i 
*/ 
template <class Sx, class St, class T>
inline void ScaledLegendrePolynomial (int n, Sx x, St t, T & values)
{
  Sx p1 = 1.0, p2 = 0.0, p3;
  St tt = t*t;
  
  if (n>=0) values[0] = 1.0;

  for (int j=1; j<=n; j++)
    {
      p3=p2; p2=p1;
      p1=((2.0*j-1.0) * x*p2 - tt*(j-1.0)*p3)/j;
      // p1= Minus (((2.0*j-1.0)/j) * Mult (x,p2) , ((j-1.0)/j) * Mult (tt,p3));
      values[j] = p1;
    }
}

template <class T> 
inline void DiffScaledLegendrePolynomial (int n, double x, double t, T & values)
{
  ArrayMem<AutoDiff<2>,10> ad_values(n+1);
  AutoDiff<2> ad_x(x, 0);
  AutoDiff<2> ad_t(t, 1);

  ScaledLegendrePolynomial(n, ad_x, ad_t, ad_values);

  for (int i = 0; i <= n; i++)
    for (int j = 0; j < 2; j++)
      values(i,j) = ad_values[i].DValue(j);
}


template <class Sx, class St, class T>
inline void ScaledIntegratedLegendrePolynomial (int n, Sx x, St t, T & values)
{
  Sx p1 = -1.0, p2 = 0.0, p3;
  St tt = t*t;
  if (n >= 0)
    values[0] = -1.0;

  for (int j=1; j<=n; j++)
    {
      p3=p2; p2=p1;
      p1=((2.0*j-3.0) * x*p2 - t*t*(j-3.0)*p3)/j;
      values[j] = p1;
    }
}



template <class S, class T>
inline void JacobiPolynomial (int n, S x, double alpha, double beta, T & values)
{
  S p1 = 1.0, p2 = 0.0, p3;

  if (n >= 0) 
    p2 = values[0] = 1.0;
  if (n >= 1) 
    p1 = values[1] = 0.5 * (2*(alpha+1)+(alpha+beta+2)*(x-1));

  for (int i  = 1; i < n; i++)
    {
      p3 = p2; p2=p1;
      p1 =
	1.0 / ( 2 * (i+1) * (i+alpha+beta+1) * (2*i+alpha+beta) ) *
	( 
	 ( (2*i+alpha+beta+1)*(alpha*alpha-beta*beta) + 
	   (2*i+alpha+beta)*(2*i+alpha+beta+1)*(2*i+alpha+beta+2) * x) 
	 * p2
	 - 2*(i+alpha)*(i+beta) * (2*i+alpha+beta+2) * p3
	 );
      values[i+1] = p1;
    }
}



template <class S, class St, class T>
inline void ScaledJacobiPolynomial (int n, S x, St t, double alpha, double beta, T & values)
{
  S p1 = 1.0, p2 = 0.0, p3;

  if (n >= 0) values[0] = 1.0;
  for (int i=0; i < n; i++)
    {
      p3 = p2; p2=p1;
      p1 =
	1.0 / ( 2 * (i+1) * (i+alpha+beta+1) * (2*i+alpha+beta) ) *
	( 
	 ( (2*i+alpha+beta+1)*(alpha*alpha-beta*beta) * t + 
	   (2*i+alpha+beta)*(2*i+alpha+beta+1)*(2*i+alpha+beta+2) * x) 
	 * p2
	 - 2*(i+alpha)*(i+beta) * (2*i+alpha+beta+2) * t * t * p3
	 );
      values[i+1] = p1;
    }
}

template <class T> 
inline void ScaledLegendrePolynomialandDiff(int n, double x, double t, T & P, T & Px, T & Pt)  
{
  /*
  ArrayMem<AutoDiff<2>,10> ad_values(n+1);
  AutoDiff<2> ad_x(x, 0);
  AutoDiff<2> ad_t(t, 1);

  ScaledLegendrePolynomial(n, ad_x, ad_t, ad_values);

  for (int i = 0; i <= n; i++)
    {
      P[i] = ad_values[i].Value();
      Px[i] = ad_values[i].DValue(0);
      Pt[i] = ad_values[i].DValue(1);
    }
  */
  if(n>=0) 
    {
      P[0] = 1.; 
      Px[0] = 0.; 
      Pt[0] = 0.; 
      if(n>=1) 
	{
	  P[1] = x;
	  Px[1] = 1.; 
	  Pt[1] = 0.; 
	} 

      double px0 = 0., px1 = 0., px2 =1.;
      double sqt = t*t; 
      for(int l = 2; l<=n; l++) 
	{ 
	  px0=px1; 
	  px1=px2;  
	  px2=  ( (2*l-1)*x*px1 - l*sqt*px0)/(l-1); 
	  Px[l] = px2; 
	  Pt[l] = -t*px1;
	  P[l] = (x*px2-sqt*px1)/l;
	}
    }

}
	  
template <class T> 
inline void LegendrePolynomialandDiff(int n, double x,  T & P, T & Px) 
{
  /*
  ArrayMem<AutoDiff<1>,10> ad_values(n+1);
  AutoDiff<1> ad_x(x, 0);
  LegendrePolynomial(n, ad_x, ad_values);

  for (int i = 0; i <= n; i++)
    {
      P[i] = ad_values[i].Value();
      Px[i] = ad_values[i].DValue(0);
    }
 
  (*testout) << "P = " << endl << P << endl
	     << "Px = " << endl << Px << endl;
  */
  if(n>=0) 
    {
      P[0] = 1.; 
      Px[0] = 0.;  
      if(n>=1) 
	{
	  P[1] = x;
	  Px[1] = 1.;  
	}       
      double px0 = 0., px1 = 0., px2 =1.;
      for(int l = 2; l<=n; l++) 
	{ 
	  px0=px1; 
	  px1=px2;  

	  px2=  ( (2*l-1)*x*px1 - l*px0)/(l-1); 
	  Px[l] = px2; 
	  P[l] = (x*px2 - px1)/l;
	}
    }
}







/*
  u(0) = 0, u(1) = 1,
  min \int_0^1 (1-x)^{DIM-1} (u')^2 dx

  representation as
  \sum c_i P_i(2x-1)
 */

template <int DIM>
class VertexExtensionOptimal
{
  enum { SIZE = 40 };
  static double coefs[SIZE][SIZE];
  static bool initialized;
public:

  VertexExtensionOptimal ();

  template <typename Tx>
  inline static Tx Calc (int p, Tx x)
  {
    Tx p1 = 1.0, p2 = 0.0, p3;
    Tx sum = 0;

    x = 2*x-1;

    if (p >= 0)
      sum += coefs[0][p];

    for (int j=1; j<=p; j++)
      {
	p3 = p2; p2 = p1;
	p1 = ((2.0*j-1.0)*x*p2 - (j-1.0)*p3) / j;
	sum += coefs[j][p] * p1;
      }

    return sum;
  }

  inline static double CalcDeriv (int p, double x)
  {
    AutoDiff<1> p1 = 1.0, p2 = 0.0, p3;
    AutoDiff<1> sum = 0;
    AutoDiff<1> adx (x, 0);  // \nabla adx = e_0

    adx = 2*adx-1;

    if (p >= 0)
      sum += coefs[0][p];

    for (int j=1; j<=p; j++)
      {
	p3 = p2; p2 = p1;
	p1 = ((2.0*j-1.0)*adx*p2 - (j-1.0)*p3) / j;
	sum += coefs[j][p] * p1;
      }

    return sum.DValue(0);
  }


  /*
    // Based on Jacobi pols, 3D only
  template <typename Tx>
  inline static Tx Calc (int p, Tx x)
  {
    ArrayMem<Tx,20> jacpol(p+1);

    JacobiPolynomial (p, 2*x-1, 1, -1, jacpol);
    
    Tx sum = 0;
    for (int j = 0; j <= p; j++)
      sum += coefs[j][p] * jacpol[j];

    return sum;
  }

  inline static double CalcDeriv (int p, double x)
  {
    ArrayMem<double,20> jacpol(p+1);

    JacobiPolynomial (p, 2*x-1, 2, 0, jacpol);
    
    double sum = 0;
    for (int j = 1; j <= p; j++)
      sum += coefs[j][p] * 0.5 * (j+1) * jacpol[j-1];

    return sum;
  }
  */
};



class VertexStandard
{
public:
  template <typename Tx>
  inline static Tx Calc (int p, Tx x)
  {
    return x;
  }

  inline static double CalcDeriv (int p, double x)
  {
    return 1;
  }
};

/**
   Compute triangle edge-shape functions

   functions vanish on upper two edges

   x,y: coordinates in triangle (-1, 0), (1, 0), (0, 1)

   f_i-2 (x, 0) = IntegratedLegendrePol_i (x)

   f_i-2 ... pol of order i, max order = n

   Monomial extension:
*/

class IntegratedLegendreMonomialExt
{
  enum { SIZE = 1000 };
  static double coefs[SIZE][2];
  
public:

  static void CalcCoeffs ()
  {
    for (int j = 1; j < SIZE; j++)
      {
	coefs[j][0] = double(2*j-3)/j;
	coefs[j][1] = double(j-3)/j;
      }
  }

  template <class Sx, class Sy, class T>
  inline static int CalcTrigExt (int n, Sx x, Sy y, T & values)
  {
    Sy fy = (1-y)*(1-y);
    Sx p3 = 0;
    Sx p2 = -1;
    Sx p1 = x;

    for (int j=2; j<=n; j++)
      {
	p3=p2; p2=p1;
	// p1=( (2*j-3) * x * p2 - (j-3) * fy * p3) / j;
	p1=( double(2*j-3)/j  * x * p2 - double(j-3)/j * fy * p3);
	// p1= coefs[j][0] * x * p2 - coefs[j][1] * fy * p3;
	values[j-2] = p1;
      }     

    return n-1;
  }


  template <class T>
  inline static int CalcTrigExtDeriv (int n, double x, double y, T & values)
  {
    double fy = (1-y)*(1-y);
    double p3 = 0, p3x = 0, p3y = 0;
    double p2 = -1, p2x = 0, p2y = 0;
    double p1 = x, p1x = 1, p1y = 0;

    for (int j=2; j<=n; j++)
      {
	p3=p2; p3x = p2x; p3y = p2y;
	p2=p1; p2x = p1x; p2y = p1y;
	double c1 = (2.0*j-3) / j;
	double c2 = (j-3.0) / j;
	
	p1  = c1 * x * p2 - c2 * fy * p3;
	p1x = c1 * p2 + c1 * x * p2x - c2 * fy * p3x;
	p1y = c1 * x * p2y - (c2 * 2 * (y-1) * p3 + c2 * fy * p3y);
	values (j-2, 0) = p1x;
	values (j-2, 1) = p1y;
      }    
    return n-1;
  }


  template <class Sx, class T>
  inline static int Calc (int n, Sx x, T & values)
  {
    Sx p3 = 0;
    Sx p2 = -1;
    Sx p1 = x;

    for (int j=2; j<=n; j++)
      {
	p3=p2; p2=p1;
	p1=( (2*j-3) * x * p2 - (j-3) * p3) / j;
	values[j-2] = p1;
      }
    return n-1;
  }

  template <class T>
  inline static int CalcDeriv (int n, double x, T & values)
  {
    double p1 = 1.0, p2 = 0.0, p3;

    for (int j=1; j<=n-1; j++)
      {
	p3 = p2; p2 = p1;
	p1 = ((2.0*j-1.0)*x*p2 - (j-1.0)*p3) / j;
	values[j-1] = p1;
      }
    return n-1;
  }
};

#endif