[R] compile error with C code and standalone R math C library

[Faheem Mitha]

Dear People,

I just went back to an old piece of C code. On trying to compile it with
the R math standalone C library I got the following error. Can anyone
enlighten me what I am doing wrong, if anything? C file (rr-sa.c) follows.

I'm on Debian sarge. I'm running R version 1.8.1. Gcc is version
3.3.1.

Thanks in advance.

                                                        Faheem.

**********************************************************************
faheem ~/co/rr/trunk>gcc -o rr rr-sa.c -lRmath -lm
/usr/lib/gcc-lib/i486-linux/3.3.2/../../../libRmath.so: undefined
reference to `Rlog1p'
collect2: ld returned 1 exit status
**********************************************************************
rr-sa.c
**********************************************************************
#include <stdio.h>
#include <assert.h>
#define T 3
#define INITVAL 1
#define MATHLIB_STANDALONE
#define THETA 2
#include <Rmath.h>
#include <time.h>

void advance(double *node, int *current_pos, double newval);
void retreat(double *node, int *current_pos);
double new_val(double *node, double  theta, int current_pos);
double accept_prob(double *node, double theta, int current_pos);
double accept_prob_fn(double a, double y);

int main()
{
  int currentpos = 0;
  int i;
  double newval, node[T+2], theta = THETA, p, u;
  set_seed(time(NULL), clock());
  node[0]=0;
  for(i=1;i<=T;i++)
    node[i]= INITVAL;    /* for simplicity choose all values the same */
  node[T+1]=0;

  for(i=0; currentpos < T; i++)
    {
      u = unif_rand();
      newval = new_val(node, theta, currentpos);
      p = accept_prob(node, theta, currentpos);
      printf("current position is %u\n",currentpos);
      printf("uniform random variable determining acceptance/rejection is %g\n",u);
      printf("value of acceptance probability is %g\n",p);
      printf("value of proposed new variable is %g\n",newval);
      if(u < p)
	advance(node,&currentpos,newval);
      else
	retreat(node, &currentpos);
    }

  for(i=1;i <= T; i++)
    printf("value of node %u is %g\n",i,node[i]);
  return 0;
}

/* function that moves chain one step forward in the event of an
   acceptance */
void advance(double *node, int *current_pos, double newval)
{
  *current_pos = *current_pos + 1;
  node[*current_pos] = newval;
}

/* function that moves the chain backwards in the event of a rejection*/
void retreat(double *node, int *current_pos)
{
  /* need special handling for small values of current_pos */
  if(*current_pos >= 1)
    *current_pos = *current_pos - 1;
  else if(*current_pos == 0)
    ;         /* do nothing, already at beginning */
}

/* generate new candidate value with appropriate distribution */
double new_val(double *node, double theta, int current_pos)
{
  double a, u, y;
  u = unif_rand();
  a = node[current_pos] + node[current_pos+2];
  if (a != 0)
    y = (1/a)*log( ( exp(a*theta) - exp(-a*theta) )*u + exp(-a*theta) );
  else  /* degenerate case (for a=0) is Unif[-theta,theta] */
    y = 2*theta*u - theta;
  return y;
}

/* generate acceptance probabilities for candidates*/
/* NB: This assumes:
(a) That all the initial values of the starting state are set to -1
(b) That theta > 1.
*/
double accept_prob(double *node, double theta, int current_pos)
{
  double a, p, y, num, denom, leftval, rtval;
  a = node[current_pos + 1];
  y = node[current_pos] + node[current_pos+2];
  num = accept_prob_fn(a, y);
 leftval = accept_prob_fn(a, -theta + node[current_pos+2]);
  rtval = accept_prob_fn(a, theta + node[current_pos+2]);
  denom = fmax2(leftval, rtval); /* function in Rmath which returns
				    max of two doubles */
  p = num/denom;
  assert(p>=0 && p <=1);
  return p;
}

/* auxilary function for accept_prob */
double accept_prob_fn(double a, double y)
{
  double f;
  if(y != 0)
  f = (exp(y) - exp(-y))/(y*exp(a*y));
  else
    f = 2;
  return f;
}