Glaring FDIV bug in Pentium!
Terje Mathisen
>Terje Mathisen (Terje.M...@hda.hydro.com) wrote:
>: After my post yesterday about how Intel improved the precision of
>: trancendental functions on the Pentium, I got an interesting email from
>: Alex Wolfe, who works for EETimes.
>
>: It seems someone reported on Compuserve that the Windows calculator gave
>: wrong results when used on a Pentium machine, but correct on a 486!
>
>: I wrote a tiny test program, p87.pas, in Borland Pascal, using Extended real
>: precision (80 bit) for all fp operations to verify this.
>
>: It looks to me like the Pentium, both the 60 and the 90MHz models, only
>: carry fp divisions to single precision.
>[deleted for brevity]
>
>This is not a feature of the compiler either!!
No, it cannot be. I hand-checked the code generated by BP very carefully,
it was using just a single FDIVP and FMULP beside the loads and stores.
>the results below are from Borland C++ under NT3.5
>32bit compile. Code is also included. Sureley this
>must be a difference between the 486 and 586, i.e.
>different FPDIV instructions what I mean this is
>intentional and not a blunder ?? and will be remedied
>by a 586 specific EFPDIV insturction or the like.
The p5 cannot require a different (and undocumented!) opcode to execute the
same FDIV as all the previous generations of x87 (co-)processors.
>
>t0 = 8.24633702441000000000e+11 :4026BFFFFFB829000000
>t1 = 1.00000000000000000000e+00 :3FFF8000000000000000
>t2 = 1.21265962489115780000e-12 :3FD7AAAAAADFDB8E4CCB
>t3 = 9.99999996274709702000e-01 :3FFEFFFFFFF000000001
>
>#include <stdio.h>
>
>void Dump(unsigned char* p)
>{
>int i;
>
>for (i=sizeof(long double)-1; i!=-1; i--)
> printf("%2.2X",(unsigned int)(*(p+i)) );
>printf("\n");
>}
>
>void main(void)
>{
>long double t0, t1, t2, t3;
>
>t0 = 824633702441.0;
>t1 = 1.0;
>t2 = t1 / t0;
>t3 = t2 * t0;
>printf("t0 = %30.20Le :",t0); Dump((unsigned char*)&t0);
>printf("t1 = %30.20Le :",t1); Dump((unsigned char*)&t1);
>printf("t2 = %30.20Le :",t2); Dump((unsigned char*)&t2);
>printf("t3 = %30.20Le :",t3); Dump((unsigned char*)&t3);
>}
I have done some more testing, it seems the FDIV is carried to about 29
bits resolution in the mantissa, which would be a very reasonable value
for single-precision microcode.
Since we seem to get the exact same (but wrong) result, independently of
the Pentium clock frequency, I suspect this could be a stupid micro-code
bug, where someone goofed and entered too few iterations for the division
algorithm.
I will try to hand-code an emulated FDIV next week, and check if any of the
intermediate bit patterns match what the Pentium returns.
If this is the problem, Intel could fix it very rapidly. A hardware problem
would take much longer, I suspect.
-Terje Mathisen (include std disclaimer) <Terje.M...@hda.hydro.com>
"almost all programming can be viewed as an exercise in caching"
Jon Jenkins
: After my post yesterday about how Intel improved the precision of
: trancendental functions on the Pentium, I got an interesting email from
: Alex Wolfe, who works for EETimes.
: It seems someone reported on Compuserve that the Windows calculator gave
: wrong results when used on a Pentium machine, but correct on a 486!
: I wrote a tiny test program, p87.pas, in Borland Pascal, using Extended real
: precision (80 bit) for all fp operations to verify this.
: It looks to me like the Pentium, both the 60 and the 90MHz models, only
: carry fp divisions to single precision.
: The tested expression is (1/824633702441)*824633702441:
: Pentium (60 & 90)
: 8.24633702441000E+0011 = 4026BFFFFFB829000000 824633702441
: 1.00000000000000E+0000 = 3FFF8000000000000000 1
: 1.21265962489116E-0012 = 3FD7AAAAAADFDB8E4CCB 1/824...
: 9.99999996274710E-0001 = 3FFEFFFFFFF000000001 (1/824..)*824...
: 486DX
: 8.24633702441000E+0011 = 4026BFFFFFB829000000
: 1.00000000000000E+0000 = 3FFF8000000000000000
: 1.21265962940867E-0012 = 3FD7AAAAAAEA8638FB73
: 1.00000000000000E+0000 = 3FFF8000000000000000
This is not a feature of the compiler either!!
the results below are from Borland C++ under NT3.5
32bit compile. Code is also included. Sureley this
must be a difference between the 486 and 586, i.e.
different FPDIV instructions what I mean this is
intentional and not a blunder ?? and will be remedied
by a 586 specific EFPDIV insturction or the like.
t0 = 8.24633702441000000000e+11 :4026BFFFFFB829000000
t1 = 1.00000000000000000000e+00 :3FFF8000000000000000
t2 = 1.21265962489115780000e-12 :3FD7AAAAAADFDB8E4CCB
t3 = 9.99999996274709702000e-01 :3FFEFFFFFFF000000001
#include <stdio.h>
void Dump(unsigned char* p)
{
int i;
for (i=sizeof(long double)-1; i!=-1; i--)
printf("%2.2X",(unsigned int)(*(p+i)) );
printf("\n");
}
void main(void)
{
long double t0, t1, t2, t3;
t0 = 824633702441.0;
t1 = 1.0;
t2 = t1 / t0;
t3 = t2 * t0;
printf("t0 = %30.20Le :",t0); Dump((unsigned char*)&t0);
printf("t1 = %30.20Le :",t1); Dump((unsigned char*)&t1);
printf("t2 = %30.20Le :",t2); Dump((unsigned char*)&t2);
printf("t3 = %30.20Le :",t3); Dump((unsigned char*)&t3);
}
----------------------------------------------------------------------
Name: Dr Jon Jenkins
Location: Digital Equipment Corp, NaC,
Burnett Place, Research Park,
Bond University, Gold Coast
QLD, AUSTRALIA 4229
Phone: 61-75-75-0151
Fax: 61-75-75-0100
Internet: jenk...@ozy.dec.com
The opinions expressed above are entirely personal and do not
reflect the corporate policy of DEC or the opinions of DEC management.
-----------------------------------------------------------------------
Andreas Kaiser
> Sureley this must be a difference between the 486 and 586, i.e.
> different FPDIV instructions what I mean this is intentional
> and not a blunder ??
Usually, the division is correct (what did you expect?). Just a few
operands are divided wrong. My results (P90) with ~25.000.000.000
random arguments (within 1..2^46), with even results divided by two
until odd, to assure unique mantissa patterns (the binary exponent
doesn't care, of course).
3221224323
12884897291
206158356633
824633702441
1443107810341
6597069619549
9895574626641
13194134824767
13194134826115
13194134827143
13194134827457
13194138356107
13194139238995
26388269649885
26388269650425
26388269651561
26388276711601
26388276712811
52776539295213
52776539301125
52776539301653
52776539307823
52776553426399
Gruss, Andreas
--------------------
-- Andreas Kaiser -- internet: a...@ananke.s.bawue.de
-------------------- fidonet: 2:246/8506.9
Jon Jenkins
: Jon Jenkins <jenk...@oils.ozy.dec.com> wrote:
: >This is not a feature of the compiler either!!
: >the results below are from Borland C++ under NT3.5
: >32bit compile. Code is also included. Sureley this
: >must be a difference between the 486 and 586, i.e.
: >different FPDIV instructions what I mean this is
: >intentional and not a blunder ?? and will be remedied
: >by a 586 specific EFPDIV insturction or the like.
: "Will be" in what sense?
What I meant was is some RISC chips have different
versions of FPDIV instructions for different FP types.
Not being initmate with Intel x87 architecture I was
wondering whether on other x87 cpus that the compiler
was converting long double to some other fp format
compatible with FPDIV on Intel x87 to be compatible
with older x87 architectures and that this problem
(which appears to be related to the bit pattern of
some operands) would be solved by 586 specific compilers
using a different (extended FPDIV) on the 586 chips ?
: As far as I know, the 586Pentium has no special EFPDIV
: instruction, only a small number of new instructions with very
: specialized uses (e.g., 8-byte compare and exchange).
As I said I cant recall ever having looked at x87 instruction
set but some cpus have both IDIV, EIDIV, FPDIV and EFPDIV
instruction types (as well as the MUL equivalents). I questioned
whether this might be the case.
: By "will be remedied" do you mean "Intel will eventually come out with a
: new version of the Pentium with such an instruction"?
Or a fix for this problem.
Dennis O'Connor~
Intel has announced that they already knew of this bug, that it
only matters to a few (probably all USENET-litterate :-) people,
and that they have changed the masks to fix the problem.
No, I don't know which stepping will have the fix, or how
to interpret the built-in version ID to see if it's a fixed
one or not. I'm an i960(R) architecture, not a Pentium(tm)
processor, guy. Sorry, really. Now you all know what I
know about the problem.
--
Dennis O'Connor doco...@sedona.intel.com
Intel i960(R) Microprocessor Division Solely responsible for what I do.
Guy Harris
>This is not a feature of the compiler either!!
>the results below are from Borland C++ under NT3.5
>32bit compile. Code is also included. Sureley this
>must be a difference between the 486 and 586, i.e.
>different FPDIV instructions what I mean this is
>intentional and not a blunder ?? and will be remedied
>by a 586 specific EFPDIV insturction or the like.
"Will be" in what sense?
As far as I know, the 586^H^H^HPentium has no special EFPDIV
instruction, only a small number of new instructions with very
specialized uses (e.g., 8-byte compare and exchange).
By "will be remedied" do you mean "Intel will eventually come out with a
Gordon Burditt
>: carry fp divisions to single precision.
So what was the Precision Control set to?
For some reason I don't understand, the default control word
on Pentium FPU uses a precision of 24 bits, and the default control
word on a 486 uses extended precision (64 bits). This is documented
but it doesn't seem to warrant mention in the "compatability" section.
They also changed around masked exceptions, too. Why? I wonder also
if they modified the Pentium to respect Precision Control more, as
it doesn't seem to make a lot of difference in the result on a 486.
>: The tested expression is (1/824633702441)*824633702441:
Chances are very good that (1/N)*N is not SUPPOSED to come
out exactly 1 when done at finite precision. Consider (1/3)*3
when done to D decimal places.
>: Pentium (60 & 90)
>: 8.24633702441000E+0011 = 4026BFFFFFB829000000 824633702441
>: 1.00000000000000E+0000 = 3FFF8000000000000000 1
>: 1.21265962489116E-0012 = 3FD7AAAAAADFDB8E4CCB 1/824...
>: 9.99999996274710E-0001 = 3FFEFFFFFFF000000001 (1/824..)*824...
Looks like it's set at single precision and it means it.
>: 486DX
>: 8.24633702441000E+0011 = 4026BFFFFFB829000000
>: 1.00000000000000E+0000 = 3FFF8000000000000000
>: 1.21265962940867E-0012 = 3FD7AAAAAAEA8638FB73
>: 1.00000000000000E+0000 = 3FFF8000000000000000
Please print the contents of the FPU control word, so we can
see what the precision was set to.
My FreeBSD 486DX system doesn't support 10-byte long doubles,
so all the printouts are in the 8-byte double format, and you
can't compare the hex to the numbers above. I used the same
program, modified a little to do the same calculation in each
of 3 precisions. Differences in precision seem to have more
effect on printf than on the calculation itself. (Note: if
you use -O, the program doesn't use the FPU at all except for
the fpsetprec() calls, as it does everything at compile time.)
Single precision
t0 = 8.24633719715391653965e+11 :4267FFFFF7052000
t1 = 1.00000000000000000000e+00 :3FF0000000000000
t2 = 1.21265963850397673340e-12 :3D75555560000000
t3 = 1.00000000000000000000e+00 :3FF0000000000000
Double precision
t0 = 8.24633702441000000000e+11 :4267FFFFF7052000
t1 = 1.00000000000000000000e+00 :3FF0000000000000
t2 = 1.21265962940866686969e-12 :3D7555555D50C71F
t3 = 9.99999999999999888978e-01 :3FEFFFFFFFFFFFFF
Extended precision
t0 = 8.24633702441000000000e+11 :4267FFFFF7052000
t1 = 1.00000000000000000000e+00 :3FF0000000000000
t2 = 1.21265962940866686969e-12 :3D7555555D50C71F
t3 = 9.99999999999999888978e-01 :3FEFFFFFFFFFFFFF
Gordon L. Burditt
sneaky.lonestar.org!gordon
#include <stdio.h>
#include <machine/floatingpoint.h>
void
Dump(unsigned char *p)
{
int i;
for (i = sizeof(long double) - 1; i != -1; i--)
printf("%2.2X", (unsigned int) (*(p + i)));
printf("\n");
}
/* main returns int, not void! */
int
main(void)
{
long double t0, t1, t2, t3;
printf("Single precision\n");
fpsetprec(FP_PS);
t0 = 824633702441.0;
t1 = 1.0;
t2 = t1 / t0;
t3 = t2 * t0;
printf("t0 = %30.20Le :", t0);
Dump((unsigned char *) &t0);
printf("t1 = %30.20Le :", t1);
Dump((unsigned char *) &t1);
printf("t2 = %30.20Le :", t2);
Dump((unsigned char *) &t2);
printf("t3 = %30.20Le :", t3);
Dump((unsigned char *) &t3);
printf("Double precision\n");
fpsetprec(FP_PD);
t0 = 824633702441.0;
t1 = 1.0;
t2 = t1 / t0;
t3 = t2 * t0;
printf("t0 = %30.20Le :", t0);
Dump((unsigned char *) &t0);
printf("t1 = %30.20Le :", t1);
Dump((unsigned char *) &t1);
printf("t2 = %30.20Le :", t2);
Dump((unsigned char *) &t2);
printf("t3 = %30.20Le :", t3);
Dump((unsigned char *) &t3);
printf("Extended precision\n");
fpsetprec(FP_PE);
t0 = 824633702441.0;
t1 = 1.0;
t2 = t1 / t0;
t3 = t2 * t0;
printf("t0 = %30.20Le :", t0);
Dump((unsigned char *) &t0);
printf("t1 = %30.20Le :", t1);
Dump((unsigned char *) &t1);
printf("t2 = %30.20Le :", t2);
Dump((unsigned char *) &t2);
printf("t3 = %30.20Le :", t3);
Dump((unsigned char *) &t3);
return 0;
}
Guy Harris
>Guy Harris (g...@nova.netapp.com) wrote:
>: Jon Jenkins <jenk...@oils.ozy.dec.com> wrote:
>: >This is not a feature of the compiler either!!
>: >the results below are from Borland C++ under NT3.5
>: >32bit compile. Code is also included. Sureley this
>: >must be a difference between the 486 and 586, i.e.
>: >different FPDIV instructions what I mean this is
>: >intentional and not a blunder ?? and will be remedied
>: >by a 586 specific EFPDIV insturction or the like.
>
>: "Will be" in what sense?
>What I meant was is some RISC chips have different
>versions of FPDIV instructions for different FP types.
>Not being initmate with Intel x87 architecture I was
>wondering whether on other x87 cpus that the compiler
>was converting long double to some other fp format
>compatible with FPDIV on Intel x87 to be compatible
>with older x87 architectures and that this problem
>(which appears to be related to the bit pattern of
>some operands) would be solved by 586 specific compilers
>using a different (extended FPDIV) on the 586 chips?
The 32-bit x86 architecture (which *I* consider to include the
floating-point instructions implemented by the 80387, as well as by the
on-chip FPUs on other x86 implementations, just as I consider the MIPS
and SPARC architectures, say, to include floating-point instructions,
even though the original implementations of those architectures put
floating-point functions on a separate chip) has a "precision control"
field in the floating-point unit status word, specifying whether the
add, subtrace, multiply, divide, and square-root instructions should run
with less than the 64-bit-significand "extended precision" default
value.
This was done, according to the 80387 manual, "to provide compatibility
with the earlier-generation arithmetic processors having less precision
than the 80387"; it permits 24-bit single-precision and 53-bit
double-precision to be used.
The floating-point add, subtract, multiply, and divide instructions
convert in-memory operands to extended precision; there are no separate
single-precision or double-precision versions of those instructions.
I.e., no, the x86 floating-point divide instructions don't have
"extended precision" variants, they all take extended-precision
operands, and do the division in a precision that depends on the
"precision control" field in the FPU status word, normally
extended-precision. (The various divide instructions differ in the way
they use the FP register stack, and whether they take an integer divisor
which is converted to a floating-point number as part of the operation,
not in the actual division they perform with the operand values.)
>: As far as I know, the 586Pentium has no special EFPDIV
>: instruction, only a small number of new instructions with very
>: specialized uses (e.g., 8-byte compare and exchange).
>As I said I cant recall ever having looked at x87 instruction
>set
Perhaps you should have done so before speculating, then...
>but some cpus have both IDIV, EIDIV, FPDIV and EFPDIV
>instruction types (as well as the MUL equivalents). I questioned
>whether this might be the case.
...because, had you done so, you would have known not to question that,
as you would have then known that it *wasn't* the case.
>: By "will be remedied" do you mean "Intel will eventually come out with a
>: new version of the Pentium with such an instruction"?
>Or a fix for this problem.
What you said was
Sureley this
must be a difference between the 486 and 586, i.e.
different FPDIV instructions what I mean this is
intentional and not a blunder ?? and will be remedied
by a 586 specific EFPDIV insturction or the like.
and, while there apparently will be "a fix for this problem", the "fix"
apparently won't be *anything* "like" of a "586 specific EFPDIV
instruction", but will be a fix to the regular boring old FP division
instructions to make them do what the manuals for the 80387, Pentium,
and probably 486 all say they should do, i.e. a correct division whose
precision depends only on the "precision control" value in the FPU
status word.
I.e., the only way in which it's "a difference between the 486 and 586"
is that the 586 had a bug in it causing it to *incorrectly* implement
floating-poit division, and was, most definitely, a blunder, not an
intentional choice on the part of Intel.
Tim Coe
>In article <CzE2...@cwi.nl>, d...@cwi.nl (Dik T. Winter) says:
>>
>>Not seen the original, but:
>> > In article <1994Nov16.1...@vitsemi.com>,
>> > Tim Coe <c...@vitsemi.com> wrote:
>> > >(4195835 / 3145727) * 4195835
^^^^^^^
Oops. The number I actually entered on the Pentium was 3145727.
I didn't check what I said as I was writing the post.
>> > >
>> > >The result was 4195579.
>> > >This represents an error of 256 or one part in ~16000.
>> >
>>I would say the error is much larger as the correct result is about 5596491.
>>I think you intended: (4195835 / 3145727) * 3145727. But in that case the
>>result tends to show 11 digits intermediate precision.
>
>On my Pentium 4195835 / 3145727 is 1.333739068902
>on a 486 4195835 / 3145727 is 1.333820449136
>
>It looks like a lot less than 11 digits of intermediate precision to me.
Since then I performed the following calculations in Microsoft
Windows Desk Calculator on a Pentium machine with the following
results:
(41.999999/35.9999999)*35.9999999 - 41.999999 ==> (-0.75)*(2^-13)
(48.999999/41.9999999)*41.9999999 - 48.999999 ==> (-1.0)*(2^-13)
(55.999999/47.9999999)*47.9999999 - 55.999999 ==> (-1.0)*(2^-13)
(62.999999/53.9999999)*53.9999999 - 62.999999 ==> (-1.0)*(2^-13)
(54.999999/59.9999999)*59.9999999 - 54.999999 ==> (-1.0)*(2^-13)
(5244795/3932159)*3932159 - 5244795 ==> (-1.0)*(2^8)
I chose these calculations in anticipation of them exposing further
Pentium FDIV failure modes. They did. The size of the erroneous results
are exactly consistant with the final version of tentive Pentium
divider model included below and in no way can be attributed to
a Desk Calculator bug. The existance of these results pins
most of the digit selection thresholds included in the model.
I also performed the following calculations that did NOT produce erroneous
results:
(38.499999/32.9999999)*32.9999999 - 38.499999 ==> 0
(45.499999/38.9999999)*38.9999999 - 45.499999 ==> 0
I have been following this thread with great interest. One misperception
that needs clearing is that this is an extended precision problem. This
bug hits between 50 and 2000 single precision dividend divisor pairs (out
of a total of 64 trillion.) Another misperception is related to the magnitude
of the relative error. I would propose the following table of probabilities
of getting the following relative errors when performing random double
extended precision divides:
relerror = (correct_result - Pentium_result)/correct_result
Error Range | Probability
-------------------------------------------
1e-4 < relerror | 0
1e-5 < relerror < 1e-4 | 0.3e-11
1e-6 < relerror < 1e-5 | 0.6e-11
1e-7 < relerror < 1e-6 | 0.6e-11
1e-8 < relerror < 1e-7 | 0.6e-11
.
.
1e-18 < relerror < 1e-17 | 0.6e-11
1e-19 < relerror < 1e-18 | 0.6e-11
Examination of the above divide failures reveals that both the dividend
and divisor are integers minus small deltas. Also notable is the induced
error is roughly delta^(2/3). The integers in the divisors are actually
restricted to those listed and their binary scalings. The integers in
the dividends may be much more freely chosen. This type of dividend
divisor pair actually occurs quite often when forward integrating
trajectories off metastable points. This is because metastable points
in systems often have certain exactly integral characteristics and as
a path diverges from the metastable point these characteristics slowly diverge
from their integral values. If the forward integration algorithm
happens to divide these characteristics, and they happen to be for
example 7 and 3, it will get nailed.
The divider model includes support for up to 60 bits of divisor and
up to 64 bits of dividend. The last four bits of dividend are kludged
in.
Here is a list of failing dividend divisor mantissas in hex. A dash
between two numbers indicates an inclusive failing range. Compile
the program and run these numbers through it and watch the bits dance:
800bf6 bffffc
a00ef6 effffc
a808d2 8fffe
e00bd2 bfffe
a7ffd2 8fffe
c3ffd2 a7ffe
dfffd2 bfffe
fbffd2 d7ffe
f9ffdc7 efffe
b9feab7-b9feabf 8fff
b9ffab0e-b9ffab7f 8fffc
-the following double extended pair fails 3 times!!!
c3ffd2eb0d2eb0d2 a7ffe
e00bd229315 bfffe
9fffef5-9fffeff effff4
9ffff21-9ffff3f effff8
9ffff4d-9ffff7f effffc
f008e35-f008e3f 8ffff4
f008e6d-f008e7f 8ffff6
f008ea1-f008ebf 8ffff8
f008ed9-f008eff 8ffffa
f008f0d-f008f3f 8ffffc
f008f45-f008f7f 8ffffe
f008f7e 8ffffff1
f0023e 8fffff8
effff0d 8ffffc
a808d1b-a808d3f 8fffe
a808d67-a808d7f 8fffe4
a808db3-a808dbf 8fffe8
a808dff 8fffec
-Tim Coe c...@vitsemi.com
#include <stdio.h>
main()
{
unsigned r0, r1, r2, r3, r4, r5, r6, s0, s1;
unsigned t0, t1, t2, t3, cycle, f, incorrect, spup;
unsigned thr_m2_m1, thr_m1_0, thr_0_1, thr_1_2, positive, errornum;
char line[30], *linepoint;
r0 = 0x0bffffc0;
r1 = 0;
r2 = 0x0800bf60;
r3 = 0;
printf("First digit of mantissas must be between 8 and f\n");
printf("Enter dividend mantissa in hex: ");
*(line+15) = '0';
scanf("%s", line);
linepoint = line;
while (*linepoint != '\0') linepoint++;
while (linepoint < line + 15) *linepoint++ = '0';
*(line+16) = '\0';
sscanf(line+15, "%x", &spup);
spup = (spup >> 2) | (12 & (spup << 2));
*(line+15) = '\0';
sscanf(line+7, "%x", &r3);
*(line+7) = '\0';
sscanf(line, "%x", &r2);
printf("Enter divisor mantissa in hex: ");
scanf("%s", line);
linepoint = line;
while (*linepoint != '\0') linepoint++;
while (linepoint < line + 15) *linepoint++ = '0';
*(line+15) = '\0';
sscanf(line+7, "%x", &r1);
*(line+7) = '\0';
sscanf(line, "%x", &r0);
r4 = 0;
r5 = 0;
t0 = r2;
while (!(t0 & 1)) t0 = t0 >> 1;
printf("%d\n", t0);
t0 = r0;
while (!(t0 & 1)) t0 = t0 >> 1;
printf("%d\n", t0);
/* These thresholds are VERY tentative. */
/* There may be bugs in them. */
t0 = r0 >> 22;
/* Next threshold is strongly indicated */
/* by the failure of 1/9895574626641 */
if (t0 < 36) thr_0_1 = 3;
/* Next threshold is strongly indicated */
/* by the failure of 1/824633702441 */
else if (t0 < 48) thr_0_1 = 4;
/* Next threshold is strongly indicated */
/* by the failure of 5244795/3932159 */
else if (t0 < 60) thr_0_1 = 5;
else thr_0_1 = 6;
thr_m1_0 = 254 - thr_0_1;
if (t0 < 33) thr_1_2 = 11;
else if (t0 < 34) {
printf("This model does not correctly handle\n");
printf("this divisor. The Pentium divider\n");
printf("undoubtly handles this divisor correctly\n");
printf("by some means that I have no evidence\n");
printf("upon which speculate.\n");
exit();
}
/* Next threshold is strongly indicated */
/* by the failure of 41.999999/35.9999999 */
else if (t0 < 36) thr_1_2 = 12;
else if (t0 < 39) thr_1_2 = 13;
/* Next threshold is strongly indicated */
/* by the failure of 1/1443107810341 and */
/* by the failure of 48.999999/41.9999999 */
else if (t0 < 42) thr_1_2 = 14;
else if (t0 < 44) thr_1_2 = 15;
/* Next threshold is strongly indicated */
/* by the failure of 55.999999/47.9999999 */
else if (t0 < 48) thr_1_2 = 16;
/* Next threshold is strongly indicated */
/* by the failure of 62.999999/53.9999999 */
else if (t0 < 54) thr_1_2 = 18;
/* Next threshold is strongly indicated */
/* by the failure of 54.999999/59.9999999 */
else if (t0 < 60) thr_1_2 = 20;
else thr_1_2 = 23;
thr_m2_m1 = 254 - thr_1_2;
if (t0 == 35) errornum = 22;
else if (t0 == 41) errornum = 26;
else if (t0 == 47) errornum = 30;
else if (t0 == 53) errornum = 34;
else if (t0 == 59) errornum = 38;
else errornum = 128;
incorrect = 0;
cycle = 1;
/* The cycle limit would be ~34 instead of */
/* 12 for double extended precision. */
while (cycle < 12) {
t0 = 255 & ((r2 >> 24) + (r4 >> 24));
if ((t0 > thr_m1_0) || (t0 < thr_0_1)) {
s0 = 0;
s1 = 0;
positive = 0;
printf("next digit 0\n");
}
else if (t0 > thr_m2_m1) {
s0 = r0;
s1 = r1;
positive = 0;
printf("next digit -1\n");
}
else if (t0 < thr_1_2) {
s0 = ~r0;
s1 = ~r1;
positive = 4;
printf("next digit 1\n");
}
else if (t0 & 128) {
s0 = (r0 << 1) | (r1 >> 31);
s1 = r1 << 1;
positive = 0;
printf("next digit -2\n");
}
else {
s0 = ~((r0 << 1) | (r1 >> 31));
s1 = ~(r1 << 1);
positive = 4;
printf("next digit 2\n");
if ((t0 == errornum) && (((r2 >> 21) & 7) == 7) && (((r4 >> 21) & 7) == 7)) {
printf("A bug condition has been detected.\n");
printf("Enter 0 for correct result or 1 for incorrect result: ");
scanf("%d", &incorrect);
if (incorrect) {
/* These amounts that are subtracted from the */
/* remainder have NOT been extensively verified. */
if (errornum == 22) s0 = s0 - (3 << 25);
else s0 = s0 - (4 << 25);
}
}
}
t0 = s0 ^ r2 ^ r4;
t1 = s1 ^ r3 ^ r5;
t2 = (s0 & r2) | (s0 & r4) | (r2 & r4);
t3 = (s1 & r3) | (s1 & r5) | (r3 & r5);
r2 = (t0 << 2) | (t1 >> 30);
r3 = t1 << 2;
r4 = (t2 << 3) | (t3 >> 29);
r5 = (t3 << 3) | positive | (spup & 3);
spup = spup >> 2;
t0 = r2;
f = 32;
while (f--) {
if (t0 & (1 << 31)) putchar('1');
else putchar('0');
t0 = t0 << 1;
}
t0 = r3;
f = 32;
while (f--) {
if (t0 & (1 << 31)) putchar('1');
else putchar('0');
t0 = t0 << 1;
}
putchar('\n');
t0 = r4;
f = 32;
while (f--) {
if (t0 & (1 << 31)) putchar('1');
else putchar('0');
t0 = t0 << 1;
}
t0 = r5;
f = 32;
while (f--) {
if (t0 & (1 << 31)) putchar('1');
else putchar('0');
t0 = t0 << 1;
}
putchar('\n');
t0 = r2 + r4;
f = 32;
while (f--) {
if (t0 & (1 << 31)) putchar('1');
else putchar('0');
t0 = t0 << 1;
}
printf(" iteration number %d\n", cycle++);
}
}
GasBoy
writes:
>>Since then I performed the following calculations in Microsoft
Windows Desk Calculator on a Pentium machine with the following
results:
(41.999999/35.9999999)*35.9999999 - 41.999999 ==> (-0.75)*(2^-13)
(48.999999/41.9999999)*41.9999999 - 48.999999 ==> (-1.0)*(2^-13)
(55.999999/47.9999999)*47.9999999 - 55.999999 ==> (-1.0)*(2^-13)
(62.999999/53.9999999)*53.9999999 - 62.999999 ==> (-1.0)*(2^-13)
(54.999999/59.9999999)*59.9999999 - 54.999999 ==> (-1.0)*(2^-13)
(5244795/3932159)*3932159 - 5244795 ==> (-1.0)*(2^8)
I chose these calculations in anticipation of them exposing further
Pentium FDIV failure modes. They did. The size of the erroneous results
are exactly consistant with the final version of tentive Pentium
divider model included below and in no way can be attributed to
a Desk Calculator bug. The existance of these results pins
most of the digit selection thresholds included in the model.<<
Hello!!???!! Haven't you read anything about the Microsoft Windows
calculator? It's just a tad buggy, in case you didn't know. If you are
trying to find fpu bugs running that program then you would accuse any
processor you ran it on. This has already been established here and
elswhere.
Scott W.
Thomas Koenig
--
Thomas Koenig, Thomas...@ciw.uni-karlsruhe.de, ig...@dkauni2.bitnet.
The joy of engineering is to find a straight line on a double
logarithmic diagram.
Christopher J. Wiggins
: calculator? It's just a tad buggy, in case you didn't know. If you are
: trying to find fpu bugs running that program then you would accuse any
: processor you ran it on. This has already been established here and
: elswhere.
Well personally I took the step of trying the Windows calculator, and found
the bug (as I already had done in a C program under Linux, as well as the
X windows calculator that came with my Slackware distribution). I was pretty
much convinced by then, but decided to take the next obvious step-use the
Windows calculator on one of the secretaries' PC (Both PC are GW2K and were
delivered within a week of each other with, as far as I can tell, the same
software distribution). The calculation this time worked, with the correct
answer produced. Hers is a 486/66, mine a P5-60. This finally convinced me
that something was indeed wrong!
Chris
Thomas Koenig
on the bug that I can see, but whose readers may be affected by it.
Followups directed to comp.sys.intel.]
There is a serious bug in the floating point performance of Intel
Pentiums, being hotly discussed on comp.sys.intel at the moment.
Readers who are interested in this topic are asked to read the
discussion there.
A short synopsis:
Intel Pentiums have a bug in their floating point arithmetic, which
reduces the accuracy of results when dividing in double precision (64
bits) from 16 to as low as 4. This error happens for certain values of
dividend and divisor, which are fairly rare. According to a model by
c...@vitsemi.com (Tim Coe), described in <1994Nov21.0...@vitsemi.com>
on comp.sys.intel, a reduction below 8 significant figures happens about
once every 5*10^10 divisions.
Some claims have been made that Intel is now shipping a fixed version
of Pentiums, which have corrected this problem. However, nobody
has reported a Pentium system for which the bug could not be
reproduced.
Here is an example C program which reproduces the bug, on a Linux
Pentium P60 with gcc 2.5.8:
#include <stdio.h>
int main()
{
double x,y,z;
x = 4195835.0;
y = 3145727.0;
z = x - (x / y) * y;
printf("%f\n",z);
return 0;
}
On a defect Pentium, this program prints "256.000000"; on a 486,
this prints "0.000000". Doing the same calculation on any other
program which relies on Pentium floating point hardware (such as Excel,
the Windows pocket calculator, or a special DOS program written to
test for the presence of the bug) will yield the same, incorrect
result.