fp-bit bugs

[Torbjorn Granlund <tege at matematik dot su dot se>]

I decided to test fp-bit.c with random bit vectors as input.  Not
surprisingly, there are a number of errors remaining, given how much errors
have been found in this code over the years.

1. 1/(-Inf) and (-1)/Inf should yield -0 but instead yield +0.
2. 7e-310/2.22 is improperly rounded.  Lots of numbers with denormalized
   input or output yield improperly rounded quotients for both
   multiplication and division.  The result is 1 ulp too small.
3. Incidentally, 7e-310*2.22 is also improperly rounded.
4. And so is 2.2805499343075318e-159/2.7653578282142482e+154.  In this
   example, the output is a denormalized number.
5. An interesting example is 5.5568968730708003e-163*8.8888643361340676e-162
   where fp-bit computes the product as 0.0, and other IEEE implementations
   4.9406564584124654e-324.  (That last number is incidentally the smallest
   positive number an IEEE double can represent.)

Here is a patch for the first problem.  I am not working on fixing the other
problems.  It probably has to do with improper handling of the sticky bit
when working on denorms.  It is not terribly serious to get 1 ulp wrong
for such numbers since denorms is a bug anyway.

1998-09-18  Torbjorn Granlund  <tege@matematik.su.se>

	* fp-bit.c (pack_d): Do not clear SIGN when fraction is 0.
	(_fpadd_parts): Get sign right for 0.

*** /home/tege/egcs-1.1a/gcc/config/fp-bit.c	Wed Jul  8 15:44:21 1998
--- /home/tege/fp-bit-new.c	Fri Sep 18 19:51:21 1998
*************** pack_d ( fp_number_type *  src)
*** 460,466 ****
    else if (fraction == 0)
      {
        exp = 0;
-       sign = 0;
      }
    else
      {
--- 460,465 ----
*************** _fpadd_parts (fp_number_type * a,
*** 735,741 ****
  	{
  	  tfraction = a_fraction - b_fraction;
  	}
!       if (tfraction > 0)
  	{
  	  tmp->sign = 0;
  	  tmp->normal_exp = a_normal_exp;
--- 734,740 ----
  	{
  	  tfraction = a_fraction - b_fraction;
  	}
!       if (tfraction >= 0)
  	{
  	  tmp->sign = 0;
  	  tmp->normal_exp = a_normal_exp;