From btyner at stat.purdue.edu  Tue Dec 19 19:58:02 2006
From: btyner at stat.purdue.edu (Benjamin Tyner)
Date: Tue, 19 Dec 2006 13:58:02 -0500
Subject: [Rd] precision when calling a C function; presence of Fortran call
Message-ID: <4588363A.5050305@stat.purdue.edu>

I'm trying to figure out why the presence of a Fortran call affects the 
result of a floating-point operation. I have C functions

void test1(int *n, double *f){
  int outC;
  double c0;

  c0 = (double) *n * *f;

  outC = floor(c0);
  printf("when f computed by R, C says %d by itself\n",outC);
}

void test2(int *n, double *f){
  extern int F77_NAME(ifloor)(double *);
  int outC,outFor;
  double c0;

  c0 = (double) *n * *f;

  outFor = F77_CALL(ifloor)(&c0);
  outC = floor(c0);
  printf("when f computed by R, C says %d, Fortran says %d\n",outC,outFor);
}

where the Fortran function ifloor is

      integer function ifloor(x)
      DOUBLE PRECISION x
      ifloor=x
      if(ifloor.gt.x) ifloor=ifloor-1
      end

void test3(){
  int outC;
  double f, c0;
  int n;

  n = 111;
  f = 40. / (double) n;

  c0 = (double) n * f;
  outC = floor(c0);
  printf("when f computed by C, C says %d by itself\n",outC);

}

void test4(){
  extern int F77_NAME(ifloor)(double *);
  int outC,outFor;
  double f, c0;
  int n;

  n = 111;
  f = 40. / (double) n;

  c0 = (double) n * f;
  outFor = F77_CALL(ifloor)(&c0);
  outC = floor(c0);
  printf("when f computed by C, C says %d, Fortran says %d\n",outC,outFor);

}

For convenience, I've put all this in a package at 
http://www.stat.purdue.edu/~btyner/test_0.1-1.tar.gz ; just install, 
load, and run test() to see the results. On my system (linux, i686), 
they are:

 > library(test)
 > test()
when f computed by R, C says 39 by itself
when f computed by R, C says 40, Fortran says 40
when f computed by C, C says 40 by itself
when f computed by C, C says 40, Fortran says 40

That is, with n=111 and f=40/111 passed in from R, test1 gives a value 
of 39. This is not a problem; I am well aware of the limitations of 
floating point. However what concerns me is that test2 gives a value of 
40. It's almost as if C precision is reduced by the presence of calling 
the Fortran function ifloor.  Furthermore, when n and f are computed in 
C, the result is always 40. So my questions are:
1. How to explain the change from 39 to 40 seemingly due to the Fortran 
call being present?
2. When Fortran is not present, why does it matter whether f is computed 
by R or C?

Thanks,
Ben



[Rd] precision when calling a C function; presence of Fortran call
Benjamin Tyner