From murdoch at stats.uwo.ca  Tue Dec 19 20:21:25 2006
From: murdoch at stats.uwo.ca (Duncan Murdoch)
Date: Tue, 19 Dec 2006 14:21:25 -0500
Subject: [Rd] precision when calling a C function;
 presence of Fortran call
In-Reply-To: <4588363A.5050305@stat.purdue.edu>
References: <4588363A.5050305@stat.purdue.edu>
Message-ID: <45883BB5.5050806@stats.uwo.ca>

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



[Rd] precision when calling a C function; presence of Fortran call
Duncan Murdoch