Re: [Rd] precision when calling a C function; presence of Fortran call

From: Duncan Murdoch <murdoch_at_stats.uwo.ca>
Date: Tue 19 Dec 2006 - 19:21:25 GMT

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 functions

On 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@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel



R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel Received on Wed Dec 20 14:26:35 2006

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.1.8, at Wed 20 Dec 2006 - 18:31:03 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-devel. Please read the posting guide before posting to the list.