code: (svn://svn.kuehne.cn/dstress/run/creal_04.d)
# int main(){
#         real x1 = 2.0;
#         real y1 = 3.0;
#         real x2 = 5.0;
#         real y2 = 7.0;
#         ireal i = 1.0i;
#
#         creal a = x1 + y1*i;
#         creal b = x2 + y2*i;
#         creal c = a / b;
#
#         creal d = ((x1*x2 + y1*y2) + i*(x2*y1 - x1*y2)) / (x2*x2 + y2*y2);
#         assert(c == d);
#
#         d = ((x1 + i*y1)*(x2 - i*y2)) / ((x2 + i*y2)*(x2 - i*y2));
#         assert(c == d);
#         return 0;
# }

Both assertions fail. This seems to be a rounding problem.

Thomas

In article <cm8cvv$1uq8$1@digitaldaemon.com>, Thomas Kuehne says...
>
>code: (svn://svn.kuehne.cn/dstress/run/creal_04.d)
># int main(){
>#         real x1 = 2.0;
>#         real y1 = 3.0;
>#         real x2 = 5.0;
>#         real y2 = 7.0;
>#         ireal i = 1.0i;
>#
>#         creal a = x1 + y1*i;
>#         creal b = x2 + y2*i;
>#         creal c = a / b;
>#
>#         creal d = ((x1*x2 + y1*y2) + i*(x2*y1 - x1*y2)) / (x2*x2 + y2*y2);
>#         assert(c == d);
>#
>#         d = ((x1 + i*y1)*(x2 - i*y2)) / ((x2 + i*y2)*(x2 - i*y2));
>#         assert(c == d);
>#         return 0;
># }
>
>Both assertions fail. This seems to be a rounding problem.

But that's expected.
you have to setup a reasonable rounding diff.
didn't we learn that on the first class of numerical analizes
(or whatever is called in english)?

Ant

Ant schrieb am Dienstag, 2. November 2004 17:54:
> In article <cm8cvv$1uq8$1@digitaldaemon.com>, Thomas Kuehne says...
>>
>>code: (svn://svn.kuehne.cn/dstress/run/creal_04.d)
>># int main(){
>>#         real x1 = 2.0;
>>#         real y1 = 3.0;
>>#         real x2 = 5.0;
>>#         real y2 = 7.0;
>>#         ireal i = 1.0i;
>>#
>>#         creal a = x1 + y1*i;
>>#         creal b = x2 + y2*i;
>>#         creal c = a / b;
>>#
>>#         creal d = ((x1*x2 + y1*y2) + i*(x2*y1 - x1*y2)) / (x2*x2 +
>>#         y2*y2); assert(c == d);
>>#
>>#         d = ((x1 + i*y1)*(x2 - i*y2)) / ((x2 + i*y2)*(x2 - i*y2));
>>#         assert(c == d);
>>#         return 0;
>># }
>>
>>Both assertions fail. This seems to be a rounding problem.
> 
> But that's expected.
> you have to setup a reasonable rounding diff.
> didn't we learn that on the first class of numerical analizes
> (or whatever is called in english)?

Don't know the english term either but yes I expected a rounding difference.

The interesting part is that:
1) doing the same in cfloat(a,b,c,d) doesn't throw an Assertion
2) the results of both fromulars (d 1. and d 2.) are identical

I know, neither on it's own is unusal but the combination smells a bit fishy.


> you have to setup a reasonable rounding diff.
Yeah, I'd like to but http://www.digitalmars.com/d/float.html says
# Rounding Control
# IEEE 754 floating point arithmetic includes the ability to set 4 different
# rounding modes. D adds syntax to access them: [blah, blah, blah] [NOTE:
# this is perhaps better done with a standard library call]

Thomas

"Thomas Kuehne" <thomas-dloop@kuehne.cn> wrote in message > The interesting part is that:
> 1) doing the same in cfloat(a,b,c,d) doesn't throw an Assertion

Floats are only 32 bits long, whereas reals are 80 bits long. Any off bits in the 68th place in the 80 bit real will be lopped off for 32 bit results.

> 2) the results of both fromulars (d 1. and d 2.) are identical

Printing them out by default rounds them to 6 digits of precision.

> I know, neither on it's own is unusal but the combination smells a bit fishy.

It's a complicated topic. If you want to investigate it further, try printing out the various values in hex format (%a).

Walter wrote:

>>1) doing the same in cfloat(a,b,c,d) doesn't throw an Assertion
> 
> Floats are only 32 bits long, whereas reals are 80 bits long. Any off bits
> in the 68th place in the 80 bit real will be lopped off for 32 bit results.

On the X86 arch, that is.

On the PowerPC, "real" is currently also 64 bit. (no "long double" type)
Eventually it will be 128 bit, when that support is less buggy in gcc...

It is enabled in gcc/gdc with the "-mlong-double-128" compile option,
but gives code generation errors on most things besides trivial code :-(

I believe it is handled in the same way as 64-bit integer instructions
are on 32-bit hardware, that is by "cheating" and dividing it into two.

--anders