FLT_ROUNDS

FLT_ROUNDS

Post by Fred J. Tydema » Fri, 23 Jan 2004 04:47:37



The C standard has in 5.2.4.2.2p6:

The rounding mode for floating-point addition is characterized by the
implementation-defined value of FLT_ROUNDS:

 -1 indeterminable
  0 toward zero
  1 to nearest
  2 toward positive infinity
  3 toward negative infinity

All other values for FLT_ROUNDS characterize implementation-defined
rounding behavior.

Now, does anyone know of any machines that have different
rounding modes for different floating-point operators?
That is, does * or / round differently than +?
If so, are any of those machines still in use?
---
Fred J. Tydeman        Tydeman Consulting

+1 (775) 287-5904      Vice-chair of J11 (ANSI "C")
Sample C99+FPCE tests: ftp://jump.net/pub/tybor/
Savers sleep well, investors eat well, spenders work forever.

 
 
 

FLT_ROUNDS

Post by pete » Fri, 23 Jan 2004 06:18:02



> The C standard has in 5.2.4.2.2p6:

> The rounding mode for floating-point addition is characterized by the
> implementation-defined value of FLT_ROUNDS:

>  -1 indeterminable
>   0 toward zero
>   1 to nearest
>   2 toward positive infinity
>   3 toward negative infinity

> All other values for FLT_ROUNDS characterize implementation-defined
> rounding behavior.

> Now, does anyone know of any machines that have different
> rounding modes for different floating-point operators?
> That is, does * or / round differently than +?
> If so, are any of those machines still in use?

I have never been able to imagine
what those things could possibley mean.

If FLT_ROUNDS is 2, does one plus zero equal two ?

--
pete

 
 
 

FLT_ROUNDS

Post by pete » Fri, 23 Jan 2004 06:20:28




> > The C standard has in 5.2.4.2.2p6:

> > The rounding mode for floating-point addition is characterized by the
> > implementation-defined value of FLT_ROUNDS:

> >  -1 indeterminable
> >   0 toward zero
> >   1 to nearest
> >   2 toward positive infinity
> >   3 toward negative infinity

> > All other values for FLT_ROUNDS characterize implementation-defined
> > rounding behavior.

> > Now, does anyone know of any machines that have different
> > rounding modes for different floating-point operators?
> > That is, does * or / round differently than +?
> > If so, are any of those machines still in use?

> I have never been able to imagine
> what those things could possibley mean.

> If FLT_ROUNDS is 2, does one plus zero equal two ?

Or rather since it deals with floats,
is zero plus zero, positive when FLT_ROUNDS is 2 ?

--
pete

 
 
 

FLT_ROUNDS

Post by Eric J. Kosteli » Fri, 23 Jan 2004 08:12:28




Quote:>Now, does anyone know of any machines that have different
>rounding modes for different floating-point operators?
>That is, does * or / round differently than +?
>If so, are any of those machines still in use?

Not quite what you're asking, but Cray C on the SV1 and similar systems
has an option of the form

   -h trunc=b

which causes the compiler to generate code to truncate the least
significant b bits of every floating-point result.  This is on
architectures with Cray (not IEEE) arithmetic, where float and double
have identical representations with 46-bit significands.  AFAIK, this
option can be specified differently for every source code file, so
various parts of the same program can run with effectively different
precision.

See www.cray.com/cgi-bin/swpubs/craydoc30/craydoc.cgi, do a search on
SV1, click on the links for Cray Standard C and C++ Reference Manuals,
and look at Sec. 2.15.7.

--Eric

 
 
 

FLT_ROUNDS

Post by Fred J. Tydema » Fri, 23 Jan 2004 09:33:47



> > If FLT_ROUNDS is 2, does one plus zero equal two ?

> Or rather since it deals with floats,
> is zero plus zero, positive when FLT_ROUNDS is 2 ?

If the hardware conforms to IEEE-754 floating-point standard, then
(+0.0) + (+0.0)
(+0.0) + (-0.0)
(-0.0) + (+0.0)
(-0.0) + (-0.0)
are all +0.0 for FLT_ROUNDS of 0, 1, and 2;
are all -0.0 for FLT_ROUNDS of 3.
---
Fred J. Tydeman        Tydeman Consulting

+1 (775) 287-5904      Vice-chair of J11 (ANSI "C")
Sample C99+FPCE tests: ftp://jump.net/pub/tybor/
Savers sleep well, investors eat well, spenders work forever.
 
 
 

FLT_ROUNDS

Post by pete » Fri, 23 Jan 2004 10:51:12




> > > If FLT_ROUNDS is 2, does one plus zero equal two ?

> > Or rather since it deals with floats,
> > is zero plus zero, positive when FLT_ROUNDS is 2 ?

> If the hardware conforms to IEEE-754 floating-point standard, then
> (+0.0) + (+0.0)
> (+0.0) + (-0.0)
> (-0.0) + (+0.0)
> (-0.0) + (-0.0)
> are all +0.0 for FLT_ROUNDS of 0, 1, and 2;
> are all -0.0 for FLT_ROUNDS of 3.

Thank you.

Which values (if any) of FLT_ROUNDS will
prevent new.c from ending it's loop ?
Are there any other implementation defined aspects which
would cause a problem in the code below ?

The code takes a float value and adds a delta to it in a loop.
Each time, the delta float value is halved.
If the value after adding the delta, is the same as the old value,
then the loop stops.

/* BEGIN new.c */

#include <stdio.h>

int main(void)
{
    double delta, old_float, new_float;

    new_float = delta = 1.0;
    do {
        delta /= 2;
        old_float = new_float;
        new_float += delta;
    } while (new_float != old_float);
    printf("%f\n", new_float);
    return 0;

Quote:}

/* END new.c */

--
pete

 
 
 

FLT_ROUNDS

Post by glen herrmannsfeld » Fri, 23 Jan 2004 15:53:23


(snip)

Quote:> Now, does anyone know of any machines that have different
> rounding modes for different floating-point operators?
> That is, does * or / round differently than +?
> If so, are any of those machines still in use?

The hexadecimal floating point from the IBM S/360 still exists in
current machines, though newer machines also support IEEE.

Arithmetic operations, such as +, -, *, /, all truncate (I think
that is what you mean by round toward zero).

Since S/370, there are LOAD ROUNDED which will convert an extended
precision number to long precision, or a long precision to short
precision, while rounding.  It would seem that if one wanted the
rounded result one could generate the longer result and round it.

(Extended, long, and short precision have 112 bit, 56 bit, and 24 bit
mantissa, respectively.)

-- glen

 
 
 

FLT_ROUNDS

Post by Andrew Nester » Tue, 27 Jan 2004 01:40:26



Quote:> The C standard has in 5.2.4.2.2p6:

> The rounding mode for floating-point addition is characterized by the
> implementation-defined value of FLT_ROUNDS:

>  -1 indeterminable
>   0 toward zero
>   1 to nearest
>   2 toward positive infinity
>   3 toward negative infinity

> All other values for FLT_ROUNDS characterize implementation-defined
> rounding behavior.

> Now, does anyone know of any machines that have different
> rounding modes for different floating-point operators?
> That is, does * or / round differently than +?
> If so, are any of those machines still in use?

Hi Fred,

A TI's TMS320C67xx micro has different rounding for
mpy and for add/sub. See http://dspvillage.ti.com,
look for C6000 CPU user guide, lit. num. SPRU189
(well, it's not exactly a *machine* like S370, yet
a 32 bit microprocessor with lots of onchip peripherals)

The processor has two multipliers M1 and M2 and
two adders, L1 and L2. It is possible not only to
set different rounding modes for the multipliers
and adders, but it is also possible to set
different rounding modes separately for each
multiplier and adder.

Thus, any of M1,M2,L1,L2 can have its own
rounding mode. See the description of FADDCR
and FMCR control registers of C6701/C6711/C6713
CPUs.

The encoding is slightly different from yours:

     0  toward nearest even
     1  toward zero (truncation)
     2  toward +Inf
     3  toward -Inf

and no indeterminable mode either.

A few notes on the program Pete's presented here earlier:

It finds the infinite sum in the form

         s = SUM(i=0,...,inf) 2**(-i)

or
         s = 1 + 0.5 + 0.25 + 0.125 + ...

This sum rapidly converges (and its limit is certainly
well known, it is equal to 1/(1-0.5) = 2, see below) and in
any situation wouldn't cause overflow. I think that
any rounding mode is fine for this loop to reach
exit condition and overall result would differ just
a couple macheps.

Hope this helps,

Andrew

P.S. Please email me if you have further questions,
since I do not look frequently on comp.arch.arithmetic.

P.P.S. A geometric progression in the form

        SUM(k=0,...,inf) a**k

converges provided |a| < 1 and has the limit of 1/(1-a).

It took 25 or 54 iterations for float and double respectively
to converge to 2.0 for the Pete's program.

- Show quoted text -

> ---
> Fred J. Tydeman        Tydeman Consulting

> +1 (775) 287-5904      Vice-chair of J11 (ANSI "C")
> Sample C99+FPCE tests: ftp://jump.net/pub/tybor/
> Savers sleep well, investors eat well, spenders work forever.

 
 
 

FLT_ROUNDS

Post by pete » Tue, 27 Jan 2004 10:30:58



>  I think that
> any rounding mode is fine for this loop to reach
> exit condition and overall result would differ just
> a couple macheps.

> Hope this helps,

> Andrew

> P.S. Please email me if you have further questions,
> since I do not look frequently on comp.arch.arithmetic.

> P.P.S. A geometric progression in the form

>         SUM(k=0,...,inf) a**k

> converges provided |a| < 1 and has the limit of 1/(1-a).

> It took 25 or 54 iterations for float and double respectively
> to converge to 2.0 for the Pete's program.

So, my real question, is about ln(), which is supposed
to return the natural log of x, for positive x.
Is the function impervious to all values of FLT_ROUNDS ?
I have confidence in the square root finding while (B > A) loop,
but the while (x != b), has me concerned.

double ln(double x)
{
    int n;
    double a, b, c;
    static double A, B, C;

    if (!A) {
        A = 1.5;
        do {
            B = A;
            A = 1 / B + B / 2;
        } while (B > A);    /* A = sqrt(2)      */
        B /= 2;             /* B = sqrt(2) / 2  */
        C = ln(A);          /* C = ln(sqrt(2))  */
    }
    for (n = 0; x > A; x /= 2) {
        ++n;
    }
    while (B > x) {
        --n;
        x *= 2;
    }
    a = (x - 1) / (x + 1);
    b = C * n + a;
    c = a * a;
    n = 1;
    do {
        x  = b;
        a *= c;
        n += 2;
        b += a / n;
    } while (x != b);
    return x * 2;

Quote:}

--
pete
 
 
 

FLT_ROUNDS

Post by Terje Mathise » Tue, 27 Jan 2004 16:37:19



> So, my real question, is about ln(), which is supposed
> to return the natural log of x, for positive x.
> Is the function impervious to all values of FLT_ROUNDS ?

I'd start by asking:

"Does this code work at all, without depedning upon some particular
compiler's implementation of the C standard?"

Quote:> I have confidence in the square root finding while (B > A) loop,
> but the while (x != b), has me concerned.

> double ln(double x)
> {
>     int n;
>     double a, b, c;
>     static double A, B, C;

>     if (!A) {

So in effect you're depending upon all static variables being
initialized to zero, or in this case, 0.0?

Quote:>         A = 1.5;
>         do {
>             B = A;
>             A = 1 / B + B / 2;
>         } while (B > A);    /* A = sqrt(2)      */
>         B /= 2;             /* B = sqrt(2) / 2  */
>         C = ln(A);          /* C = ln(sqrt(2))  */

And here you're initializing C with a recursive call to the remainder of
the same function, which is is going to use C, except that it is
currently 0.0? Does that even work? :-)

Quote:>     }
>     for (n = 0; x > A; x /= 2) {
>         ++n;
>     }
>     while (B > x) {
>         --n;
>         x *= 2;
>     }
>     a = (x - 1) / (x + 1);
>     b = C * n + a;

So on the first (recursive) call, when you're trying to determine
ln(sqrt(2)), C will be 0.0, right?

What's wrong with simply stating this directly instead?

   static double A = 1.4142135623730951,
                 B = 0.70710678118654757,
                 C = 0.3465735902799727;

Quote:>     c = a * a;
>     n = 1;
>     do {
>         x  = b;
>         a *= c;
>         n += 2;
>         b += a / n;
>     } while (x != b);
>     return x * 2;
> }

Terje

--

"almost all programming can be viewed as an exercise in caching"

 
 
 

FLT_ROUNDS

Post by Andrew Nester » Tue, 27 Jan 2004 21:16:46




> > Which values (if any) of FLT_ROUNDS will
> > prevent new.c from ending it's loop ?
> > Are there any other implementation defined aspects which
> > would cause a problem in the code below ?

> > The code takes a float value and adds a delta to it in a loop.
> > Each time, the delta float value is halved.
> > If the value after adding the delta, is the same as the old value,
> > then the loop stops.

> > /* BEGIN new.c */

> > #include <stdio.h>

> > int main(void)
> > {
> >     double delta, old_float, new_float;

> >     new_float = delta = 1.0;
> >     do {
> >         delta /= 2;
> >         old_float = new_float;
> >         new_float += delta;
> >     } while (new_float != old_float);
> >     printf("%f\n", new_float);
> >     return 0;
> > }

> > /* END new.c */

> Assume IEEE-754 floating-point.
> If FLT_ROUNDS is 2 (towards +INF), then delta will never become zero,
> and new_float will always keep incrementing towards +INF.  After a
> very long time (after creating every floating-point number from 2.0
> to DBL_MAX), new_float will overflow to infinity, at which point
> the loop will exit (after having done about 2**62 iterations). For
> FLT_ROUNDS being 0, 1, or 3, the loop should exit after 53 iterations.
> So, the loop will exit with FLT_ROUNDS being 0, 1, 2, or 3; faithful
> IEEE-754 being done, and the compiler not generating bad code.

Correct, while I was wrong in my previous post saying all modes give
the same result. Round up mode destroys the iteration process.

In fact, the inequality condition *is* unsafe in this and many other
situations as a termination criterion of an iteration process.

A safe condition can be chosen for either precision, single or double
in the following manner:

        while (delta >= EPSILON*new_float);

where EPSILON is either _FLT_ or DBL_ EPSILON defined in <float.h>

This condition is based on the standard assumption that adding
anything below machine epsilon to 1 is 1 and works in any rounding
mode.

I also have to note here that, while the above is true for x87 FPU,
other machines may deviate from the standard. E.g. C67xx FPU
in round up mode will not round up delta to the smallest representable
number when after the current division by 2 all the significand bits
are zeroes, it will return zero as the result. Thus, for double, the
loop would stop after 2**10 iterations.

--
Andrew

- Show quoted text -

> ---
> Fred J. Tydeman        Tydeman Consulting

> +1 (775) 287-5904      Vice-chair of J11 (ANSI "C")
> Sample C99+FPCE tests: ftp://jump.net/pub/tybor/
> Savers sleep well, investors eat well, spenders work forever.

 
 
 

FLT_ROUNDS

Post by pete » Tue, 27 Jan 2004 21:41:44




> > So, my real question, is about ln(), which is supposed
> > to return the natural log of x, for positive x.
> > Is the function impervious to all values of FLT_ROUNDS ?

> I'd start by asking:

> "Does this code work at all, without depedning upon some particular
> compiler's implementation of the C standard?"

> > I have confidence in the square root finding while (B > A) loop,
> > but the while (x != b), has me concerned.

> > double ln(double x)
> > {
> >     int n;
> >     double a, b, c;
> >     static double A, B, C;

> >     if (!A) {

> So in effect you're depending upon all static variables being
> initialized to zero, or in this case, 0.0?

It's C code. That's the rules of C.

- Show quoted text -

Quote:> >         A = 1.5;
> >         do {
> >             B = A;
> >             A = 1 / B + B / 2;
> >         } while (B > A);    /* A = sqrt(2)      */
> >         B /= 2;             /* B = sqrt(2) / 2  */
> >         C = ln(A);          /* C = ln(sqrt(2))  */

> And here you're initializing C with a recursive
> call to the remainder of
> the same function, which is is going to use C, except that it is
> currently 0.0? Does that even work? :-)

> >     }
> >     for (n = 0; x > A; x /= 2) {
> >         ++n;
> >     }
> >     while (B > x) {
> >         --n;
> >         x *= 2;
> >     }
> >     a = (x - 1) / (x + 1);
> >     b = C * n + a;

> So on the first (recursive) call, when you're trying to determine
> ln(sqrt(2)), C will be 0.0, right?

Yes. At that point, x will be equal to A,
which means that n will be equal to zero,
and (b = C * n + a) will be equal to a.

Quote:> What's wrong with simply stating this directly instead?

>    static double A = 1.4142135623730951,
>                  B = 0.70710678118654757,
>                  C = 0.3465735902799727;

The function is intended to use all of the available precision,
regardless of how much is available.

Quote:> >     c = a * a;
> >     n = 1;
> >     do {
> >         x  = b;
> >         a *= c;
> >         n += 2;
> >         b += a / n;
> >     } while (x != b);
> >     return x * 2;
> > }

--
pete
 
 
 

FLT_ROUNDS

Post by pete » Tue, 27 Jan 2004 21:52:23





> > > Which values (if any) of FLT_ROUNDS will
> > > prevent new.c from ending it's loop ?
> > > Are there any other implementation defined aspects which
> > > would cause a problem in the code below ?

> > > The code takes a float value and adds a delta to it in a loop.
> > > Each time, the delta float value is halved.
> > > If the value after adding the delta, is the same as the old value,
> > > then the loop stops.

> > > /* BEGIN new.c */

> > > #include <stdio.h>

> > > int main(void)
> > > {
> > >     double delta, old_float, new_float;

> > >     new_float = delta = 1.0;
> > >     do {
> > >         delta /= 2;
> > >         old_float = new_float;
> > >         new_float += delta;
> > >     } while (new_float != old_float);
> > >     printf("%f\n", new_float);
> > >     return 0;
> > > }

> > > /* END new.c */

> > Assume IEEE-754 floating-point.
> > If FLT_ROUNDS is 2 (towards +INF), then delta will never become zero,
> > and new_float will always keep incrementing towards +INF.  After a
> > very long time (after creating every floating-point number from 2.0
> > to DBL_MAX), new_float will overflow to infinity, at which point
> > the loop will exit (after having done about 2**62 iterations). For
> > FLT_ROUNDS being 0, 1, or 3, the loop should exit after 53 iterations.
> > So, the loop will exit with FLT_ROUNDS being 0, 1, 2, or 3; faithful
> > IEEE-754 being done, and the compiler not generating bad code.

> Correct, while I was wrong in my previous post saying all modes give
> the same result. Round up mode destroys the iteration process.

> In fact, the inequality condition *is* unsafe in this and many other
> situations as a termination criterion of an iteration process.

I know. That's the trick. I'm interested in examples
where the direct comparison of two floats is appropriate.
So far, the only interesting example that I have, is
with the greater than operator, in a square root finding loop,
an example of which, is in the ln() function definition.

- Show quoted text -

Quote:> A safe condition can be chosen for either precision, single or double
> in the following manner:

>         while (delta >= EPSILON*new_float);

> where EPSILON is either _FLT_ or DBL_ EPSILON defined in <float.h>

> This condition is based on the standard assumption that adding
> anything below machine epsilon to 1 is 1 and works in any rounding
> mode.

> I also have to note here that, while the above is true for x87 FPU,
> other machines may deviate from the standard. E.g. C67xx FPU
> in round up mode will not round up delta to the smallest representable
> number when after the current division by 2 all the significand bits
> are zeroes, it will return zero as the result. Thus, for double, the
> loop would stop after 2**10 iterations.

Thank you.

--
pete

 
 
 

FLT_ROUNDS

Post by Terje Mathise » Wed, 28 Jan 2004 01:36:41




>>So on the first (recursive) call, when you're trying to determine
>>ln(sqrt(2)), C will be 0.0, right?

> Yes. At that point, x will be equal to A,
> which means that n will be equal to zero,
> and (b = C * n + a) will be equal to a.

>>What's wrong with simply stating this directly instead?

>>   static double A = 1.4142135623730951,
>>                 B = 0.70710678118654757,
>>                 C = 0.3465735902799727;

> The function is intended to use all of the available precision,
> regardless of how much is available.

I realized that that had to be the case, but a _very_ long decimal
expansion would still suffice, or wouldn't it? (Those I used was
intended to be  good for about 64-bit mantissas.)

I.e. you are programming with 'double' here, not arbitrary precision. ;-)

Terje

--

"almost all programming can be viewed as an exercise in caching"

 
 
 

FLT_ROUNDS

Post by James Kuyp » Wed, 28 Jan 2004 02:13:04




...
> > I have confidence in the square root finding while (B > A) loop,
> > but the while (x != b), has me concerned.

> > double ln(double x)
> > {
> >     int n;
> >     double a, b, c;
> >     static double A, B, C;

> >     if (!A) {

> So in effect you're depending upon all static variables being
> initialized to zero, or in this case, 0.0?

6.7.8p10: "If an object that has static storage duration is not
initialized explicitly, then: ... if it has arithmetic type, it is
initialized to (positive or unsigned) zero;"

Quote:> >         A = 1.5;
> >         do {
> >             B = A;
> >             A = 1 / B + B / 2;
> >         } while (B > A);    /* A = sqrt(2)      */
> >         B /= 2;             /* B = sqrt(2) / 2  */
> >         C = ln(A);          /* C = ln(sqrt(2))  */

> And here you're initializing C with a recursive call to the remainder of
> the same function, which is is going to use C, except that it is
> currently 0.0? Does that even work? :-)

Of course it does. The call to ln() is recursive only if A is 0.0. By
the time the recursive call is made, A has already been set to a
non-zero value, so the recursion terminates with that second call. C
will have a value of 0.0 during that recursive call, which apparantly
doesn't cause a problem. It acquires the correct value of ln(sqrt(2))
as a result of evaluating the call.

I can't vouch for the validity of the algorithm being used (actually,
I could, but I don't want to spend the time), but it does indeed give
correct results for the values I tried (1.0, 2.0, and 4.0).

- Show quoted text -

Quote:> >     }
> >     for (n = 0; x > A; x /= 2) {
> >         ++n;
> >     }
> >     while (B > x) {
> >         --n;
> >         x *= 2;
> >     }
> >     a = (x - 1) / (x + 1);
> >     b = C * n + a;

> So on the first (recursive) call, when you're trying to determine
> ln(sqrt(2)), C will be 0.0, right?

> What's wrong with simply stating this directly instead?

>    static double A = 1.4142135623730951,
>                  B = 0.70710678118654757,
>                  C = 0.3465735902799727;

Nothing. These are two different valid ways to do it. The method he
gave only produces the wrong result if the code itself is wrong; your
approach also gives a wrong result if the initialization of the static
data is wrong; it's debateable which form is more reliable. Yours is
certainly faster, but is less portable, because you'll need to add
additional digits when porting to machines with smaller values for
DBL_EPSILON. For C99 code, hexadecimal constants would be preferred;
this is precisely the kind of situation for which hexadecimal floating
point constants were invented.