Benjamin Tyner
2006-Dec-19 18:58 UTC
[Rd] precision when calling a C function; presence of Fortran call
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
Duncan Murdoch
2006-Dec-19 19:21 UTC
[Rd] precision when calling a C function; presence of Fortran call
On 12/19/2006 1:58 PM, Benjamin Tyner wrote:> I'm trying to figure out why the presence of a Fortran call affects the > result of a floating-point operation. I have C functionsOn Windows, it's fairly common for runtime libraries to switch the precision from 80 bit to 64 bit. R on Windows tries to guard against this, and on your test package, I get when f computed by R, C says 40 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 It's possible that the floating point change is happening on your platform, or happened before you even got to this call. Here's code that works in Windows to detect such a change across a LoadLibrary call: rcw = _controlfp(0,0) & ~_MCW_IC; /* Infinity control is ignored */ _clearfp(); tdlh = LoadLibrary(path); dllcw = _controlfp(0,0) & ~_MCW_IC; No idea if _controlfp() exists on your system. Duncan Murdoch> > 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 > > ______________________________________________ > R-devel at r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-devel