Fixed range reduction for large multiples of pi/2 on systems with

broken assignment to floats (e.g., i386 with gcc -O, but not amd64 or
ia64; i386 with gcc -O0 worked accidentally).

Use an unnamed volatile temporary variable to trick gcc -O into clipping
extra precision on assignment.  It's surprising that only 1 place needed
to be changed.

For tanf() on i386 with gcc -O, the bug caused errors > 1 ulp with a
density of 2.3% for args larger in magnitude than 128*pi/2, with a
maximum error of 1.624 ulps.

After this fix, exhaustive testing shows that range reduction for
floats works as intended assuming that it is in within a factor of
about 2^16 of working as intended for doubles.  It provides >= 8
extra bits of precision for all ranges.  On i386:

range                       max error in double/single ulps    extra precision
-----                       -------------------------------    ---------------
0 to 3*pi/4                 0x000d3132  /  0.0016              9+ bits
3*pi/4 to 128*pi/2          0x00160445  /  0.0027              8+
128*pi/2 to +Inf            0x00000030  /  0.00000009          23+
128*pi/2 up, -O0 before fix 0x00000030  /  0.00000009          23+
128*pi/2 up, -O1 before fix 0x10000000  /  0.5                 1

The 23+ bits of extra precision for large multiples corresponds to almost
perfect reduction to a pair of floats (24 extra would be perfect).

After this fix, the maximum relative error (relative to the corresponding
fdlibm double precision function) is < 1 ulp for all basic trig functions
on all 2^32 float args on all machines tested:

          amd64     ia64      i386-O0   i386-O1
	  ------    ------    ------    ------
cosf:     0.8681    0.8681    0.7927    0.5650
sinf:     0.8733    0.8610    0.7849    0.5651
tanf:     0.9708    0.9329    0.9329    0.7035

1 files changed, 1 insertions, 0 deletions

diff --git a/lib/msun/src/k_rem_pio2f.c b/lib/msun/src/k_rem_pio2f.c
index e1e9e35..958c863 100644
--- a/lib/msun/src/k_rem_pio2f.c
+++ b/lib/msun/src/k_rem_pio2f.c
@@ -169,6 +169,7 @@ recompute:
 	    case 2:
 		fw = 0.0;
 		for (i=jz;i>=0;i--) fw += fq[i];
+		fw = *(volatile float *)&fw;	/* clip any extra precision */
 		y[0] = (ih==0)? fw: -fw;
 		fw = fq[0]-fw;
 		for (i=1;i<=jz;i++) fw += fq[i];