www.digitalmars.com         C & C++   DMDScript  

digitalmars.D.bugs - creal & opDiv

reply Thomas Kuehne <thomas-dloop kuehne.cn> writes:
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
Nov 02 2004
parent reply Ant <Ant_member pathlink.com> writes:
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
Nov 02 2004
parent reply Thomas Kuehne <thomas-dloop kuehne.cn> writes:
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.

# 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
Nov 02 2004
parent reply "Walter" <newshound digitalmars.com> writes:
"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).
Nov 02 2004
parent =?ISO-8859-1?Q?Anders_F_Bj=F6rklund?= <afb algonet.se> writes:
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
Nov 05 2004