bisect.c

language: C
license: GPL 2

Code for Snippet:

                
#include <R.h>
#include <Rinternals.h>
 
#define FOREVER for(;;)
 
SEXP            getListElement(SEXP, char *);
SEXP            anyFC(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP);
 
SEXP
BISECT(SEXP f, SEXP A, SEXP B, SEXP ITMAX, SEXP EPS, SEXP VERBOSE, SEXP env)
{
    int             itmax = INTEGER(ITMAX)[0], verbose = LOGICAL(VERBOSE)[0], itel = 0;
    SEXP            R_fcall = R_NilValue, X = R_NilValue;
    double          x, a = REAL(A)[0], b = REAL(B)[0], eps = REAL(EPS)[0], fx, fa, fb;
    PROTECT(R_fcall = lang2(f, R_NilValue));
    PROTECT(X = allocVector(REALSXP, 1));
    SETCADR(R_fcall, A);
    fa = REAL(eval(R_fcall, env))[0];
    SETCADR(R_fcall, B);
    fb = REAL(eval(R_fcall, env))[0];
    FOREVER {
        x = (a + b) / 2.0;
        REAL(X)[0] = x;
        SETCADR(R_fcall, X);
        fx = REAL(eval(R_fcall, env))[0];
        if (fx < 0) {
            if (fa < 0) {
                a = x;
                fa = fx;
            } else {
                b = x;
                fb = fx;
            }
        } else {
            if (fa < 0) {
                b = x;
                fb = fx;
            } else {
                a = x;
                fa = fb;
            }
        }
        if (verbose > 0)
            Rprintf("%3d %20.15f %20.15f %20.15f %20.15f\n", itel, a, b, x, fx);
        if ((fabs(a - b) < eps) || (itel == itmax))
            break;
        itel++;
    }
    UNPROTECT(2);
    return X;
}
 
comments powered by Disqus

Info

Link to this snippet:


Download to Code Collector

To use the direct link to your snippet on CodeCollector.net either copy the html from the above section or drag the Download to Code Collector to where you would like to use it.

More Info:

Times Viewed: 549
Date Added: 2013-03-07 16:45:17
Last Modified: 2013-04-17 22:12:19

Web Analytics