Velocity Reviews - Computer Hardware Reviews

Velocity Reviews > Newsgroups > Programming > C Programming > Are multiplication and division of double exactly defined under C99with __STDC_IEC_559__ ?

Reply
Thread Tools

Are multiplication and division of double exactly defined under C99with __STDC_IEC_559__ ?

 
 
Francois Grieu
Guest
Posts: n/a
 
      06-09-2010
[reposted with tweaks, and code that compiles]

Assume a conforming C99 implementation that defines __STDC_IEC_559__.
"double" is a floating point type capable of exactly representing all
integers in the range [0..0x1FFFFFFFFFFFFF] (that's 53 bit). Also assume
FE_TONEAREST is in effect.

Can it be said that multiplication and division of double are exactly
defined? Or/and that when x and y are in range [0..0x1FFFFFFFFFFFFF],
the result of the following functions mym() and myd() is fully specified ?

Note: that would seem to imply over 105 bits of internal precision on
intermediary results before the rounding occurs.

TIA,

Francois Grieu


#include <stdlib.h>
#include <math.h>
#include <fenv.h>
#include <stdio.h>

#ifndef __STDC_IEC_559__
// #error "No claim IEC 60559 FP is available"
#endif

#ifndef FE_TONEAREST
#error "error: FE_TONEAREST is not available"
#endif

typedef unsigned long long tu53;

// question: is that portable?
// it would imply implies a 106-bit wide multiplier output
tu53 mym(tu53 x, tu53 y)
{
return llrint( ( ((double)1/(1ull<<53)) * x ) * y );
}

// question: is that portable?
// it would imply implies a 106-bit wide multiplier output
tu53 myd(tu53 x, tu53 y)
{
return y > x ? llrint( ( (double)(1ull<<53) * x ) / y ) : 0;
}

tu53 myatou53(const char *inStr)
{
long long x;
return inStr!=NULL && (x = atoll(inStr))>=0 && x<1ull<<53 ? x : 0;
}

// for each pair of arguments x y, sghow x y mym(x,y) myd(x,y)
int main(int argc, char *argv[])
{
int j;
if (fesetround(FE_TONEAREST)==0)
for( j=2; j<argc; j+= 2 )
{
tu53 x,y;
x = myatou53(argv[j-1]);
y = myatou53(argv[j]);
// <OT> if that does not work under Win32, change ll to I64 </OT>
printf("#\t%llu\t%llu\t%llu\t%llu\n",
x, y, mym(x,y), myd(x,y) );
}
return EXIT_SUCCESS;
}
 
Reply With Quote
 
 
 
 
Thad Smith
Guest
Posts: n/a
 
      06-10-2010
Francois Grieu wrote:
> [reposted with tweaks, and code that compiles]
>
> Assume a conforming C99 implementation that defines __STDC_IEC_559__.
> "double" is a floating point type capable of exactly representing all
> integers in the range [0..0x1FFFFFFFFFFFFF] (that's 53 bit). Also assume
> FE_TONEAREST is in effect.
>
> Can it be said that multiplication and division of double are exactly
> defined?


After a single multiply or divide, followed by a cast to double or storage to a
double variable, the resulting expression should be determined by the IEC 559
specification.

> Or/and that when x and y are in range [0..0x1FFFFFFFFFFFFF],
> the result of the following functions mym() and myd() is fully specified ?


If the input parameters are constrained to 53 bits or less and the arithmetic
results are an integer of 64 or fewer bits, I believe so. I think calling the
llrint() function, even if internal inline, forces a rounding to strict double
precision prior to integer conversion.

> Note: that would seem to imply over 105 bits of internal precision on
> intermediary results before the rounding occurs.


No. Using a shift and conditional add algorithm on normalized values, unused
low-order bits can be safely discarded. I think two guard bits, updated with
each iteration, are sufficient (with single bit multiply cycles) to distinguish
0, <0.5, 0.5, and >0.5 LSB for correct rounding.

--
Thad
 
Reply With Quote
 
 
 
 
Francois Grieu
Guest
Posts: n/a
 
      06-10-2010
On 10/06/2010 06:53, Thad Smith wrote:
> Francois Grieu wrote:
>> [reposted with tweaks, and code that compiles]
>>
>> Assume a conforming C99 implementation that defines __STDC_IEC_559__.
>> "double" is a floating point type capable of exactly representing all
>> integers in the range [0..0x1FFFFFFFFFFFFF] (that's 53 bit). Also assume
>> FE_TONEAREST is in effect.
>>
>> Can it be said that multiplication and division of double are exactly
>> defined?

>
> After a single multiply or divide, followed by a cast to double or
> storage to a double variable, the resulting expression should be
> determined by the IEC 559 specification.


That's also my understanding of the rehash of the standard that I read
here and there.

>> Or/and that when x and y are in range [0..0x1FFFFFFFFFFFFF],
>> the result of the following functions mym() and myd() is fully
>> specified ?

>
> If the input parameters are constrained to 53 bits or less and the
> arithmetic results are an integer of 64 or fewer bits,

[that's hopefully the case in my test code]
> I believe so. I
> think calling the llrint() function, even if internal inline, forces a
> rounding to strict double precision prior to integer conversion.
>
>> Note: that would seem to imply over 105 bits of internal precision on
>> intermediary results before the rounding occurs.

>
> No. Using a shift and conditional add algorithm on normalized values,
> unused low-order bits can be safely discarded. I think two guard bits,
> updated with each iteration, are sufficient (with single bit multiply
> cycles) to distinguish 0, <0.5, 0.5, and >0.5 LSB for correct rounding.


Indeed, I missed that. More generally if the multiplication is performed
in n successive multiply/add/round steps, we can live with ceil(53/n)+53
bits of intermediary result, or something on that tune.
I guess something similar is possible with division.

Thanks!

Francois Grieu
 
Reply With Quote
 
 
 
Reply

Thread Tools

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
Re: __STDC_IEC_559__ (defined or !defined ?) Keith Thompson C Programming 0 08-17-2010 04:36 PM
Are multiplication and division of double exactly defined under C99with __STDC_IEC_559__ ? Francois Grieu C Programming 0 06-08-2010 11:59 AM
Source of term "multiplication" in matrix multiplication William Hughes C Programming 13 03-15-2010 02:04 PM
cannot convert parameter from 'double (double)' to 'double (__cdecl *)(double)' error Sydex C++ 12 02-17-2005 06:30 PM
need help with large int multiplication and division akickdoe22@hotmail.com C++ 1 01-21-2005 02:46 AM



Advertisments