Just been trying to get the polygon triangulation code from this package: http://www-2.cs.cmu.edu/~quake/triangle.html to dyn.load into R. Uh oh. Floating point exceptions. Track it down to some FPU diddling that the author deems is necessary. Here's my minimal code that breaks: flipme.c: #include <fpu_control.h> void flipme(){ int cword; cword=4210; _FPU_SETCW(cword); printf("%f\n",sqrt(2.0)); } Compile with "R SHLIB flipme.c", then dyn.load and run: > dyn.load("flipme.so") > .C("flipme") Process R floating point exception at Mon Nov 18 17:41:50 2002 The debugger tells me its the printf line that kills it. Anyone have an idea what's going on here? I've tried a few things: * using _FPU_GETCW I've determined that R hasn't changed the control word from the system default. * I wrote a short main program that used dlopen() to call flipme - and it worked fine, so the problem isn't intrinsic to dynamically loaded code. * the author of the triangle package has put this in the code that manipulates the FPU control word: /* Don't change this routine unless you fully understand it. */ - I dont understand it, I changed it, the code broke. Anyone with deeper understanding of floating point, iee754, and assorted other Recipes For Chaos (RFCs) care to comment before we go back to using a file or a pipe to communicate between R and the code? Baz
i386 floating point tweaking
2 messages · Barry Rowlingson, Duncan Murdoch
On Mon, 18 Nov 2002 17:57:47 +0000, you wrote in message <3DD92A1B.1020707@lancaster.ac.uk>:
Just been trying to get the polygon triangulation code from this package: http://www-2.cs.cmu.edu/~quake/triangle.html to dyn.load into R. Uh oh. Floating point exceptions. Track it down to some FPU diddling that the author deems is necessary. Here's my minimal code that breaks: flipme.c: #include <fpu_control.h> void flipme(){ int cword; cword=4210; _FPU_SETCW(cword); printf("%f\n",sqrt(2.0)); } Compile with "R SHLIB flipme.c", then dyn.load and run:
dyn.load("flipme.so")
.C("flipme")
Process R floating point exception at Mon Nov 18 17:41:50 2002
Dynamic code that changes the floating point word should change it
back when it returns. R depends on it being correct.
The value you're using is hex 1072, which means among other things
that calculations should be done to 24 bit precision, and various
exceptions should not be masked. This won't work for R, which assumes
full precision and masking of various things. The safe thing to do
would be to save the CW on entry to each of these routines, set it to
1072 for the duration, then on exit clear any pending exceptions, and
restore the R setting.
I don't know how to do that in C. In assembler it would be
fstcw saved_cw
do some stuff
fclex
fldcw saved_cw
Duncan Murdoch