"Dik T. Winter" <Di********@cwi.nl> wrote in message news:<I5********@cwi.nl>...

In article <62**************************@posting.google.com > da**************@gmail.com (Daniel Vallstrom) writes: > da**************@gmail.com (Daniel Vallstrom) wrote in message news:<62**************************@posting.google. com>... > > [snip question asking how to compute the following loop in a fast manner]

> >

> > for ( ; n != 0; n-- )

> > {

> > x0 = nextafter( x0, DBL_MAX );

> > }

My suggestion for the loop would be breaking up the floating-point number in

its constituent parts, and work with them. In that case you add n (suitably

scaled) to the mantissa and recombine. Be careful with overflow in the

mantissa when you add... See frexp and ldexp (at least on my system).

Thanks again. Here is a simple solution:

int x0Exp;

frexp( x0, &x0Exp );

x0 += ldexp( n * nextafter( 0.0, DBL_MAX ), x0Exp - DBL_MIN_EXP );

I'm still miffed that I didn't find the above solution directly

(eventhough hindsight is always 20/20). I'll post a full program with

the solution at the end.

Daniel Vallstrom

// Tests consecutive nextafter calls and ways to compute the same thing fast

// without iterated nextafter calls.

// Daniel Vallstrom, 041019.

// Compile with e.g: gcc -std=c99 -pedantic -Wall -O3 -lm nextafterNTest.c

/* The program should print something like:

n: 0x1000000 (16777216)

x0: 0x1.001p+8

xn: n nextafter( x0, DBL_MAX ) iterations: 0x1.0010001p+8

The xn loop took 0.8 cpu seconds.

m: n * nextafter( 0.0, DBL_MAX ): 0x0.0000001p-1022

e: exp2( floor( log2( x0 ) ) ): 0x1p+8

z: e | m: 0x1.0000001p+8

z1: z - e: 0x1p-20

z2: x0 + z1: 0x1.0010001p+8 (==xn)

x0Exp: frexp( x0, &x0Exp ): 9

mnorm: DBL_MIN + m: 0x1.0000001p-1022

zb: ldexp( mnorm, x0Exp-DBL_MIN_EXP ): 0x1.0000001p+8 (==z)

zb1: zb - e: 0x1p-20 (==z1)

zb2: x0 + zb1: 0x1.0010001p+8 (==xn)

zc1: ldexp( m, x0Exp-DBL_MIN_EXP ): 0x1p-20 (==z1)

zc2: x0 + zc1: 0x1.0010001p+8 (==xn)

*/

#include <stdio.h>

#include <float.h>

#include <math.h>

#include <assert.h>

#include <time.h>

int main( void )

{

#if FLT_RADIX != 2

#error Case not supported.

#endif

clock_t clockAtLoopStart; // Used to time the start of loops.

clock_t clockAtLoopEnd; // Used to time the end of loops.

// The number of nextafter iterations to use. Can be anything.

unsigned int n = ( 1u << 24 );

printf( "\nn: %#x (%u)\n", n, n );

// The base double to apply nextafter to. Can be anything.

double x0 = 0x1.001p8;

printf( "x0: %a\n", x0 );

// The nth double greater than x0.

double xn = x0;

clockAtLoopStart = clock();

for ( unsigned int k = 0; k != n; k++ )

{

xn = nextafter( xn, DBL_MAX );

}

clockAtLoopEnd = clock();

printf( "xn: n nextafter( x0, DBL_MAX ) iterations: %a\n", xn );

// xn: n nextafter( x0, DBL_MAX ) iterations: 0x1.0010001p+8

printf( " The xn loop took %.1f cpu seconds.\n",

(double) ( clockAtLoopEnd - clockAtLoopStart ) / CLOCKS_PER_SEC );

// The xn loop took 0.8 cpu seconds.

// The goal is to compute xn fast. At least gcc -O3 does *not* optimize

// iterated nextafter calls (the above loop). If the loop (n) is long (or

// many loops have to be done) it takes too long to run.

// Calculate a useful mantissa:

double m = nextafter( 0.0, DBL_MAX );

m = n*m;

printf( "m: n * nextafter( 0.0, DBL_MAX ): %a\n", m );

// m: n * nextafter( 0.0, DBL_MAX ): 0x0.0000001p-1022

// Calculate a useful exponent part.

double e = exp2( floor( log2( x0 ) ) );

printf( "e: exp2( floor( log2( x0 ) ) ): %a\n", e );

// e: exp2( floor( log2( x0 ) ) ): 0x1p+8

// Combine e and m.

double z = e;

for ( unsigned char * zPtr = &z, * zEnd = zPtr + sizeof(z), * mPtr = &m;

zPtr != zEnd; zPtr++, mPtr++ )

{

*zPtr |= *mPtr;

}

printf( "z: e | m: %a\n", z );

// z: e | m: 0x1.0000001p+8

// Remove the faulty first mantissa bit in z, the one coming from e.

double z1 = z - e;

printf( "z1: z - e: %a\n", z1 );

// z1: z - e: 0x1p-20

// Finally we are ready to add the value to x0.

double z2 = x0 + z1;

printf( "z2: x0 + z1: %a (==xn)\n", z2 );

// z2: x0 + z1: 0x1.0010001p+8 (==xn)

assert( z2 == xn );

// Here is an approach using frexp and ldexp instead of bit-oring.

int x0Exp;

frexp( x0, &x0Exp );

printf( "x0Exp: frexp( x0, &x0Exp ): %d\n", x0Exp );

// x0Exp: frexp( x0, &x0Exp ): 9

// Normalize the mantissa.

double mnorm = DBL_MIN + m;

printf( "mnorm: DBL_MIN + m: %a\n", mnorm );

// mnorm: DBL_MIN + m: 0x1.0000001p-1022

// Increase the exponent of mnorm to match the exponent of x0.

double zb = ldexp( mnorm, x0Exp-DBL_MIN_EXP );

printf( "zb: ldexp( mnorm, x0Exp-DBL_MIN_EXP ): %a (==z)\n", zb );

// zb: ldexp( mnorm, x0Exp-DBL_MIN_EXP ): 0x1.0000001p+8 (==z)

assert( zb == z );

// The rest is just like above.

// Remove the faulty first mantissa bit in zb.

assert( e == ldexp( 1.0, x0Exp-1 ) );

double zb1 = zb - e;

printf( "zb1: zb - e: %a (==z1)\n", zb1 );

// zb1: zb1 - e: 0x1p-20

assert( zb1 == z1 );

// We are ready to add the value to x0.

double zb2 = x0 + zb1;

printf( "zb2: x0 + zb1: %a (==xn)\n", zb2 );

// zb2: x0 + zb1: 0x1.0010001p+8 (==xn)

assert( zb2 == xn );

// Here is a better solution.

double zc1 = ldexp( m, x0Exp-DBL_MIN_EXP );

printf( "zc1: ldexp( m, x0Exp-DBL_MIN_EXP ): %a (==z1)\n", zc1 );

// zc1: ldexp( m, x0Exp-DBL_MIN_EXP ): 0x1p-20 (==z1)

assert( zc1 == z1 );

// Add the value to x0.

double zc2 = x0 + zc1;

printf( "zc2: x0 + zc1: %a (==xn)\n", zc2 );

// zc2: x0 + zc1: 0x1.0010001p+8 (==xn)

assert( zc2 == xn );

return 0;

}