PC Review


Reply
Thread Tools Rating: Thread Rating: 1 votes, 1.00 average.

Excel and the Math Coprocessor for DLLs

 
 
Lynn McGuire
Guest
Posts: n/a
 
      5th Apr 2012
On 4/5/2012 9:48 AM, joeu2004 wrote:
> "Martin Brown" <|||newspam|||@nezumi.demon.co.uk> wrote:
>> If it allows embedded assembly then the following
>> code in VC will do the test explicitly if inserted
>> before "return chptst"

>
> Let's not get side-tracked by methods to make the one example behave consistently and as expected. AFAIK, that is not the issue that
> Lynn is dealing with. (Although he is the final arbiter on that point.)
>
> IIRC, the issue is: Lynn discovered that (some) floating-point computations (initially not the chptst example) did not return exactly
> the same result when called from VBA and when linked into a program. He probably expected an exact zero, which made the difference
> "obvious" to him.
>
> I suspect Lynn interpreted that as a defect, and in Googling for an explanation, he stumbled across the decades-old and long-since
> fixed FDIV defect. When Lynn implemented the FDIV example (incorrectly), he discovered that it did indeed return zero in his program,
> but not exactly zero (albeit infinitesimally close) when called from VBA.
>
> Incorrectly thinking that this explained the original disparity that he encountered, Lynn became fixated on this one example.
>
> But my point is: any example would do. So there is no point in taking draconian steps to make the one example -- the chptst
> computation -- work.
>
> Moreover, although it has been interesting to explore the behavoir of the FPU control word with respect to VBA v. Watcom and VC++, I
> believe that is not the right solution for Lynn.
>
> Instead of setting the FPCW precision mode, Lynn should recognize that such floating-point inconsistencies are common place. For
> example, it is not unusual to experience them when doing "the same" calculation in Excel and in VBA.
>
> Again, the solution is to explicitly round the results of calculation when you expect results to be consistent to a specific precision.
>
> Alternatively, design algorithms to tolerate infinitesimal differences in floating-point calculations. For example, instead of test
> A=B, test ABS(A-B)<e, where "e" is some difference that reflects your tolerance for error, e.g. 0.005.


Nope, I have had the FDIV check in our software
since 1995. We had several customers all of
sudden reporting bad results with their new PCs.
Here is the comment that I put in the code for
that subroutine:

C 03/27/95 Lynn McGuire add pentium math coprocessor test

That test used to work just fine. I suspect that
the Watcom floating point code optimization has
gotten better since 1995. Since I don't have a
bad FPU to check this code with anymore, I
probably need to pull it out of our code.

Lynn
 
Reply With Quote
 
 
 
 
joeu2004
Guest
Posts: n/a
 
      5th Apr 2012
"Lynn McGuire" <(E-Mail Removed)> wrote:
> Nope, I have had the FDIV check in our software
> since 1995.


Okay. My bad! In that case, Martin's assembly-language implementation is
just the thing for you. Your original test design was incorrect; it worked
only by coincidence.

 
Reply With Quote
 
 
 
 
Martin Brown
Guest
Posts: n/a
 
      6th Apr 2012
On 05/04/2012 21:06, joeu2004 wrote:
> "Lynn McGuire" <(E-Mail Removed)> wrote:
>> Nope, I have had the FDIV check in our software
>> since 1995.

>
> Okay. My bad! In that case, Martin's assembly-language implementation is
> just the thing for you. Your original test design was incorrect; it
> worked only by coincidence.


Strictly it worked OK in a standard Fortran REAL*8 aka 53bit mantissa
environment provided that the compiler doesn't get too clever. The store
to an 8 byte real is harmless there.

However, everything breaks if the arithmetic 64bit mantissa but only 53
bits stored the result ends up different in the least significant bit.

This should not matter at all unless your algorithm is extremely
sensitive to roundoff - you still have almost 16 good decimal digits.

All decent commercial compilers are smart enough to optimise divide by a
constant to multiple by a reciprocal these days. So the test is not
actually testing the divide instruction at all any more.

--
Regards,
Martin Brown
 
Reply With Quote
 
joeu2004
Guest
Posts: n/a
 
      7th Apr 2012
"Martin Brown" <|||newspam|||@nezumi.demon.co.uk> wrote:
> provided that the compiler doesn't get too clever


__Of_course__, all of my comments assume that compile-time optimizations are
avoided.

That goes without saying, if our intent is to test the FDIV defect or the
effect of difference precision modes (_PC_64, _PC_53, _PC_24).

(Although it was useful to Lynn that you discovered that his compiler is
optimizing the code.)


"Martin Brown" <|||newspam|||@nezumi.demon.co.uk> wrote:
> On 05/04/2012 21:06, joeu2004 wrote [in response to Lynn]:
>> Your original test design was incorrect;
>> it worked only by coincidence.

>
> Strictly it worked OK in a standard Fortran REAL*8
> aka 53bit mantissa environment [...]. The store to
> an 8 byte real is harmless there.
>
> However, everything breaks if the arithmetic 64bit
> mantissa but only 53 bits stored the result ends up
> different in the least significant bit.


When I said that Lynn's original chptst implementation was incorrect, I was
referring to the rounding of the __intermediate__ result (only the division)
to 64-bit representation (53-bit mantissa)

For the FDIV defect, if x is 4195835 and y is 3145727, then (x/y)*y-x is
exactly zero whether all calculation is done __consistently__ using 64-bit
or using 80-bit representation (the latter with a 64-bit mantissas). Of
course, the final result is stored with 64-bit representation.

(And again, __of_course__ I am only talking about the case when compile-time
optimizations are avoided.)

But instead of a single expression, Lynn[*] implemented this as div = x/y,
then chptst = div*y-x.

That worked (i.e. chptst is exactly zero) for him in the past only because
apparently he was using 64-bit representation __consistently__.

It failed to work (i.e. chptst is not exactly zero) when Excel VBA used
80-bit representation, but only because Lynn is not using 80-bit
representation __consistently__.

Indeed, it is the store into div, with only 53-bit precision, that causes
the problem. It is not "harmless" in this particular case.

That is what I meant when I said it worked (in the past) only by
coincidence. The "coincidence" was that he was using __consistent__
precision in the past, namely 64-bit representation (53-bit precision).



"Martin Brown" <|||newspam|||@nezumi.demon.co.uk> wrote:
> This should not matter at all unless your algorithm
> is extremely sensitive to roundoff - you still have
> almost 16 good decimal digits.


It is not clear to me what your point is.

But in a previous response, I alluded to the fact that _in_general__, we
cannot expect (x/y)*y-x to be exactly zero, even when we use __consistent__
precision for the calculation.

(And again, __of_course__ I am only talking about the case when compile-time
optimizations are avoided.)

Here are some random examples, all using x and y with the same absolute and
relative magnitudes as Nicely's example for the FDIV defect. By that, I
mean: x and y are 7-digit numbers; and 1 < INT(x/y) < 2.

80bit 64bit mixMode x y
=0 =0 <>0 4195835 3145727 (FDIV bug)
=0 =0 =0 1300666 1233646
=0 <>0 <>0 1695563 1538366
=0 <>0 =0 none
<>0 =0 <>0 1923832 1204810
<>0 =0 =0 none
<>0 <>0 <>0 1867447 1462980
<>0 <>0 =0 none


-----[*] I don't know if Lynn invented the div/chptst implementation or he copied
it from some reference. When I asked him for a reference, he pointed me
only to the wiki FDIV bug page,
http://en.wikipedia.org/wiki/Pentium_FDIV_bug.

I do not see Lynn's div/chptst implementation on that page or on any
referenced page.

On the contrary, I do find Nicely's original pentbug.c file, which contains
a __consistent__ implementation using 64-bit representation. Excerpted:

double lfNum1=4195835.0, lfNum2=4195835.0, lfDenom1=3145727.0,
lfDenom2=3145727.0, lfQuot, lfProd, lfZero;
/* The duplicate variables are used in an effort to foil compiler
optimization and compile-time evaluation of numeric literals. */
lfQuot=lfNum1/lfDenom1;
lfProd=lfDenom2*lfQuot;
lfZero=lfNum2 - lfProd;

Aside #1.... Note that Nicely went out of his way to avoid compile-time
optimizations. Or so he hoped ;-).

If Lynn had followed Nicely's implementation correctly, he wouldn't have had
any problem with the FDIV constants despite the change in precision mode
between Watcom languages and Excel VBA and despite the Watcom compile-time
optimizations.

Aside #2.... IMHO, Nicely was wrong to use 64-bit representation to
demonstrate the FDIV bug. Luckily, it worked, but only because the FDIV
error was so large, about 8E-5. He might have overlooked more subtle errors
due to rounding the intermediate 80-bit results to 64-bit representation.

 
Reply With Quote
 
joeu2004
Guest
Posts: n/a
 
      7th Apr 2012
Errata.... I wrote:
> If Lynn had followed Nicely's implementation
> correctly, he wouldn't have had any problem
> with the FDIV constants despite the change in
> precision mode between Watcom languages and
> Excel VBA


But perhaps Nicely's implementation or something nearly identical was not
available when Lynn wrote his(?) algorithm.

The pentbug.c file that I find today has only a 2011 copyright. But it
presumably was written in or before 2003, since that is the date attributed
to the pentbug.zip file on Nicely's website.

If the algorithm had been made public at some earlier date (e.g. 1995), I
would think that Nicely would know that the earlier date(s) of "publication"
should be included in the copyright.


PS.... I wrote:
> IMHO, Nicely was wrong to use 64-bit representation
> to demonstrate the FDIV bug. Luckily, it worked,
> but only because the FDIV error was so large, about
> 8E-5. He might have overlooked more subtle errors
> due to rounding the intermediate 80-bit results to
> 64-bit representation.


I meant to add.... I also think that it was "wrong" to use (x/y)*y-x to
demonstrate the FDIV bug in the first place.

Well, I do agree that it had "marketing appeal" insofar as it might have
made the defect more relevant to some people.

But my point is.... IMHO, that is not the correct way to test for the
presence of the defect.

Actually, it is Tim Coe who discovered the error with 4195835.0/3145727.0
per se. It represents the largest known error among the many defective FDIV
results.

And it does appear that Nicely or Coe did look at the result of the division
using 80-bit arithmetic, at least eventually.

Nicely's published results of the division are (the comma shows 15
significant digits to the left):

1.33382044913624,10025 (Correct value)
1.33373906890203,75894 (Flawed Pentium)

Note that the representation of the "correct value" is more precise than the
64-bit representation can be. The 64-bit representation is exactly
1.33382044913624,10930577198087121360003948211669921875, somewhat greater
than 1.333...,10025. The next closest 64-bit representation is exactly
1.33382044913624,0871013114883680827915668487548828125, which is too low.

Based on the information available today, I would implement the FDIV bug
test as follows.

#include <math.h>

double ident(double x) { return x; } // thwart compiler optimizations

int fdivTest() // 1 = okay; 0 = error
{
const double num = 4195835;
const double denom = 3145727;
const double quot = 1.33382044913624;
return (fabs(ident(num)/ident(denom) - quot) < 5e-15);
}

It does not matter whether that is evaluated consistently in 64-bit (53-bit
mantissa) or 80-bit (64-bit mantissa) representation or a mix.

FYI, we cannot expect the division result to have the same 64-bit binary
representation as the constant 1.33382044913624.

In fact, the two closest representations of the constant are
1.33382044913623,9982834695183555595576763153076171875 and
1.33382044913624,02048793001085869036614894866943359375.




 
Reply With Quote
 
 
 
Reply

Thread Tools
Rate This Thread
Rate This Thread:

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are Off


Similar Threads
Thread Thread Starter Forum Replies Last Post
Excel and the Math Coprocessor for DLLs Lynn McGuire Microsoft Excel Discussion 55 7th Apr 2012 05:26 AM
Excel is messing with the math coprocessor Lynn McGuire Microsoft Excel Programming 6 16th Mar 2012 04:34 PM
Excel is messing with the math coprocessor Lynn McGuire Microsoft Excel Discussion 6 16th Mar 2012 04:34 PM
excel math forumla? (simple math problem inside for math people!) Jason Microsoft Excel Discussion 3 16th Feb 2006 11:54 AM
math coprocessor? =?Utf-8?B?R2F5bGU=?= Microsoft Access 9 18th May 2005 03:19 PM


Features
 

Advertising
 

Newsgroups
 


All times are GMT +1. The time now is 03:39 AM.