Message-ID: <320E4A9F.6E8@pobox.oleane.com>
Date: Sun, 11 Aug 1996 23:03:27 +0200
From: Francois Charton <deef AT pobox DOT oleane DOT com>
Organization: CCMSA
MIME-Version: 1.0
To: djgpp AT delorie DOT com
Subject: problem with libm.a "pow" function (and comment on ran())
Content-Type: text/plain; charset=us-ascii
Content-Transfer-Encoding: 7bit

I posted something last week on FPU exceptions compiling povray 3.0 with 
-lm. 
I tried since to isolate the bug (if it is one)
Here is a small program :

MYPROG.C
#include <math.h>

void main(void)
{
int i;
for(i=1;i<1000000;i++) pow(i/1000000.0,90.0);
}

/* EOF */

When I compile it with  :
gcc -o myprog myprog.c
(I link the libc.a pow function), everything works great

but if I say
gcc -o myprog myprog.c -lm
I get (at i==257) a Floating Point exception error

This is pretty strange as my program only calculate power 90 of values 
strictly between 0 and 1, so the result HAVE to be betweeen 0 and 1, and 
there can be no overflow.

As for underflows, as pow(x,90) is an increasing function of x, it should 
rather happen for i=1...

Weirder still, if I start my for(;;) loop at i=256 (instead of 1) it 
bombs at 259, if I start at 257, it bombs at 260... for i over 300 
everything is allright.

I'm afraid there is something wrong in the implementation of libm.a 
(which seems to originate from SUN). Has anyone seen something like this 
before ?


While I'm at it, I looked yesterday at the code for rand(), found the 
"classical" ANSI linear congruential generator with a small change : 
there is some kind of "bit shuffling" of the result:
 next = (next<<15) ^ (next >> 27);

Are you sure that this doesn't degrade the randomness of the generator ? 
If I remember well, higher bits in such generators are more random than 
lower bits, and the shuffling is more or less an "exchange" between 
higher and lower bits.

Regards,
Francois