[Rcpp-devel] Looping over the parameters of a function

[Petre Caraiani]

Hello everybody,
I have a code in C which computes the root of a function using the Brent
algorithm. The code is attached below. I provided the full code but you can
focus on the function definition and the main program. I can call this
program from R. What I am interested in is looping over the parameters of
the defined function. For example, passing an array a to the C program and
computing the root of the function: (x*x-a[i]/2) for each entry in vector a.


#include <Rcpp.h>
using namespace Rcpp;

#include <stdio.h>
#include <math.h>

#define FALSE  0
#define TRUE   1


//Test function for Brent method
double AFunction(double x) {
  return (x*x-1/2);
}

// TRUE if x1*x2 negative
int RootBracketed(double x1,double x2) {
  int result;
  if ((x1 > 0 && x2 > 0) || (x1 < 0 && x2 < 0))
    result = FALSE;
  else
    result = TRUE;
  return result;
}

// returns the minimum of two real numbers
double Minimum(double x1,double x2) {
  double result;
  if (x1 < x2) result = x1;
  else result = x2;
  return result;
}

/****************************************************
*              Brent Method Function                *
* ------------------------------------------------- *
* The purpose is to find a real root of a real      *
* function f(x) using Brent's method.               *
*                                                   *
* INPUTS:  x1,x2     : interval of root             *
*          Tolerance : desired accuracy for root    *
*          maxIter   : maximum number of iterations *
*                                                   *
* OUTPUTS: The function returns the root value      *
*          ValueAtRoot : value of f(root)           *
*          niter    : number of done iterations     *
*          error    : =0, all OK                    *
*                   : =1, no root found in interval *
*                   : =2, no more iterations !      *
****************************************************/
double BrentRoots( double x1, double x2, double Tolerance,
                   int maxIterations,
               double *valueAtRoot,
                   int *niter, double *error )  {

  double FPP = 1e-11, nearzero = 1e-20;

  double result, AA, BB, CC, DD, EE, FA, FB, FC, Tol1, PP, QQ, RR, SS, xm;
  int i, done;

  i = 0; done = FALSE;   error = 0;
  AA = x1;  BB = x2;  FA = AFunction(AA); FB = AFunction(BB);
  if (!(RootBracketed(FA,FB)))
    *error = 1;
  else {
    FC = FB;
    do {
      if (!(RootBracketed(FC,FB))) {
        CC = AA; FC = FA; DD = BB - AA; EE = DD;
      }
      if (fabs(FC) < fabs(FB)) {
        AA = BB; BB = CC; CC = AA;
        FA = FB; FB = FC; FC = FA;
      }
      Tol1 = 2.0 * FPP * fabs(BB) + 0.5 * Tolerance;
      xm = 0.5 * (CC-BB);
      if ((fabs(xm) <= Tol1) || (fabs(FA) < nearzero)) {
        result = BB;
        done = TRUE;
        *valueAtRoot = AFunction(result);
      } // A root has been found
      else {
        if ((fabs(EE) >= Tol1) && (fabs(FA) > fabs(FB))) {
          SS = FB/ FA;
          if (fabs(AA - CC) < nearzero) {
            PP = 2.0 * xm * SS;
            QQ = 1.0 - SS;
          }
          else {
            QQ = FA/FC;
            RR = FB /FC;
            PP = SS * (2.0 * xm * QQ * (QQ - RR) - (BB-AA) * (RR - 1.0));
            QQ = (QQ - 1.0) * (RR - 1.0) * (SS - 1.0);
          }
          if (PP > nearzero) QQ = -QQ;
          PP = fabs(PP);
          if ((2.0 * PP) < Minimum(3.0*xm *QQ-fabs(Tol1 * QQ), fabs(EE *
QQ))) {
            EE = DD;  DD = PP/QQ;
          }
          else {
            DD = xm;   EE = DD;
          }
        }
        else {
          DD = xm;
          EE = DD;
        }
        AA = BB;
        FA = FB;
        if (fabs(DD) > Tol1)
          BB = BB + DD;
        else {
          if (xm > 0) BB = BB + fabs(Tol1);
          else BB = BB - fabs(Tol1);
        }
        FB = AFunction(BB);
        i++;
      }
  }  while ((!done) && (i < maxIterations));
    if (i >= maxIterations) *error = 2;
  }
  *niter = i;
  return result;
} // BrentRoots()
// End of file zbrent.cpp



// Below is a simple example of exporting a C++ function to R. You can
// source this function into an R session using the Rcpp::sourceCpp
// function (or via the Source button on the editor toolbar)

// For more on using Rcpp click the Help button on the editor toolbar

// [[Rcpp::export]]

 NumericVector fun560 (NumericVector a){
//int n=a.size();

double   e,x1,x2,err,yr,r;
int      maxiter,k;

x1=0.0;
x2=1.0;
e =1e-10;
maxiter = 10;



a[1] = BrentRoots(x1,x2,e,maxiter,&yr,&k,&err);
a[2] = BrentRoots(x1,x2,e,maxiter,&yr,&k,&err);

//for (int ii = 0; ii <= nrow-1; ii++) {

//          Am(ii,0)= a[ii]*(1.0+B(ii,0));

//     for (int jj = 0; jj<= ncol-2; jj++){

//        Am(ii,jj+1) = Am(ii,jj)*(1.0+B(ii,jj+1));

return a;
 }
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20140408/09262f8b/attachment.html>