/*
 *  nde_bp.c by A. Pepelyshev  Copyright (C) 2011
 *  This version distributed under GPL (version 2 or later)
 */

/*
====    nonparametric density estimation using Bernstein polynomials    ====
*/

#include <math.h>

#ifndef max
#  define max(a,b) ((a) > (b) ? (a) : (b))
#  define min(a,b) ((a) < (b) ? (a) : (b))
#endif

#define abso(a) ((a) > 0 ? (a) : -(a) )


#if !defined(PI)
#  define PI 3.14159265
#endif


#include <R.h>
#include <Rinternals.h>
#include <Rdefines.h>
  
SEXP
bernstein_poly_vec(SEXP n, SEXP x )
{
  int i, j, nn;
  double xx, xx1, *pbern, *pn, *px;
  SEXP bern;

    PROTECT(n = AS_NUMERIC(n));
    PROTECT(x = AS_NUMERIC(x));
    pn = NUMERIC_POINTER(n); 
    px = NUMERIC_POINTER(x);
    xx=px[0];
    xx1 = 1-xx;
    nn=pn[0];
    nn=abso(nn);
    PROTECT(bern = NEW_NUMERIC(nn+1));
    pbern = NUMERIC_POINTER(bern);

  if(nn == 0)
  {
    pbern[0] = 1.0;
  }
  else
  {  
   pbern[0] = xx1;
   pbern[1] = xx;
   if(nn>1)
   {
    for (i = 2; i <= nn; i++) 
    {
      pbern[i] = xx * pbern[i-1];
      for( j=(i-1); j>=1; j-- )
      {
        pbern[j] = xx * pbern[j-1] + xx1 * pbern[j];
      }
      pbern[0] = xx1 * pbern[0];
    }
   }
  }
 UNPROTECT(3);
 return(bern);
}