<div dir="ltr">Hello everybody,<br>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.<br>
<br><br>#include <Rcpp.h><br>using namespace Rcpp;<br><br>#include <stdio.h><br>#include <math.h><br><br>#define FALSE  0<br>#define TRUE   1<br><br><br>//Test function for Brent method<br>double AFunction(double x) {<br>
  return (x*x-1/2);<br>}<br><br>// TRUE if x1*x2 negative<br>int RootBracketed(double x1,double x2) {<br>  int result;<br>  if ((x1 > 0 && x2 > 0) || (x1 < 0 && x2 < 0)) <br>    result = FALSE;<br>
  else<br>    result = TRUE;<br>  return result;<br>}<br><br>// returns the minimum of two real numbers<br>double Minimum(double x1,double x2) {<br>  double result;<br>  if (x1 < x2) result = x1;<br>  else result = x2;<br>
  return result;<br>}<br><br>/****************************************************<br>*              Brent Method Function                *<br>* ------------------------------------------------- *<br>* The purpose is to find a real root of a real      *<br>
* function f(x) using Brent's method.               *<br>*                                                   *<br>* INPUTS:  x1,x2     : interval of root             *<br>*          Tolerance : desired accuracy for root    *<br>
*          maxIter   : maximum number of iterations *<br>*                                                   *<br>* OUTPUTS: The function returns the root value      *<br>*          ValueAtRoot : value of f(root)           *<br>
*          niter    : number of done iterations     *<br>*          error    : =0, all OK                    *<br>*                   : =1, no root found in interval *<br>*                   : =2, no more iterations !      *<br>
****************************************************/<br>double BrentRoots( double x1, double x2, double Tolerance,<br>                   int maxIterations,<br>               double *valueAtRoot,<br>                   int *niter, double *error )  {<br>
<br>  double FPP = 1e-11, nearzero = 1e-20;<br><br>  double result, AA, BB, CC, DD, EE, FA, FB, FC, Tol1, PP, QQ, RR, SS, xm;<br>  int i, done;<br><br>  i = 0; done = FALSE;   error = 0;<br>  AA = x1;  BB = x2;  FA = AFunction(AA); FB = AFunction(BB);<br>
  if (!(RootBracketed(FA,FB))) <br>    *error = 1;<br>  else {<br>    FC = FB;<br>    do {<br>      if (!(RootBracketed(FC,FB))) {<br>        CC = AA; FC = FA; DD = BB - AA; EE = DD;<br>      }<br>      if (fabs(FC) < fabs(FB)) {<br>
        AA = BB; BB = CC; CC = AA;<br>        FA = FB; FB = FC; FC = FA;<br>      }<br>      Tol1 = 2.0 * FPP * fabs(BB) + 0.5 * Tolerance;<br>      xm = 0.5 * (CC-BB);<br>      if ((fabs(xm) <= Tol1) || (fabs(FA) < nearzero)) {<br>
        result = BB;<br>        done = TRUE;<br>        *valueAtRoot = AFunction(result);<br>      } // A root has been found<br>      else {<br>        if ((fabs(EE) >= Tol1) && (fabs(FA) > fabs(FB))) {<br>
          SS = FB/ FA;<br>          if (fabs(AA - CC) < nearzero) {<br>            PP = 2.0 * xm * SS;<br>            QQ = 1.0 - SS;<br>          }<br>          else {<br>            QQ = FA/FC;<br>            RR = FB /FC;<br>
            PP = SS * (2.0 * xm * QQ * (QQ - RR) - (BB-AA) * (RR - 1.0));<br>            QQ = (QQ - 1.0) * (RR - 1.0) * (SS - 1.0);<br>          }<br>          if (PP > nearzero) QQ = -QQ;<br>          PP = fabs(PP);<br>
          if ((2.0 * PP) < Minimum(3.0*xm *QQ-fabs(Tol1 * QQ), fabs(EE * QQ))) {<br>            EE = DD;  DD = PP/QQ;<br>          }<br>          else {<br>            DD = xm;   EE = DD;<br>          }<br>        }<br>
        else {<br>          DD = xm;<br>          EE = DD;<br>        }<br>        AA = BB;<br>        FA = FB;<br>        if (fabs(DD) > Tol1) <br>          BB = BB + DD;<br>        else {<br>          if (xm > 0) BB = BB + fabs(Tol1);<br>
          else BB = BB - fabs(Tol1);<br>        }<br>        FB = AFunction(BB);<br>        i++;<br>      }<br>  }  while ((!done) && (i < maxIterations));<br>    if (i >= maxIterations) *error = 2;<br>  }<br>
  *niter = i;<br>  return result;<br>} // BrentRoots()<br>// End of file zbrent.cpp<br><br><br><br>// Below is a simple example of exporting a C++ function to R. You can<br>// source this function into an R session using the Rcpp::sourceCpp <br>
// function (or via the Source button on the editor toolbar)<br><br>// For more on using Rcpp click the Help button on the editor toolbar<br><br>// [[Rcpp::export]]<br><br> NumericVector fun560 (NumericVector a){<br>//int n=a.size();<br>
<br>double   e,x1,x2,err,yr,r;<br>int      maxiter,k;<br><br>x1=0.0;<br>x2=1.0;<br>e =1e-10;<br>maxiter = 10;<br><br><br><br>a[1] = BrentRoots(x1,x2,e,maxiter,&yr,&k,&err);<br>a[2] = BrentRoots(x1,x2,e,maxiter,&yr,&k,&err);<br>
<br>//for (int ii = 0; ii <= nrow-1; ii++) {<br>  <br>//          Am(ii,0)= a[ii]*(1.0+B(ii,0));<br>          <br>//     for (int jj = 0; jj<= ncol-2; jj++){  <br>       <br>//        Am(ii,jj+1) = Am(ii,jj)*(1.0+B(ii,jj+1));<br>
<br>return a;<br> }   <br></div>