Re: Irix6 long doubles implemented wrong? (27_io/ostream_inserter_arith)

[Alexandre Oliva <aoliva@redhat.com>]

On Jan 26, 2003, Alexandre Oliva <aoliva@redhat.com> wrote:

> On Jan 19, 2003, Alexandre Oliva <aoliva@redhat.com> wrote:
>> the minimum normal exponent should be set such that all 107
>> significant bits are available (and it's 107, not 106, as I
>> implemented it) for normals.

> It can't possibly be 107.

On second (or rather nth) thought, it possibly can, if you take
advantage of the sign bit of the second double.  But then you have to
allow it to be different from the sign of the first double, which we
don't, and the conversion routines would have to be tweaked to take
advantage of it without losing the additional bit of precision.

Furthermore, a way to implement most of the emulation without having
to use integers has occurred to me.  Here's the idea.

For addition of two long-doubles composed of a pair of doubles

- add the most significant double of each number.  If the result is
  infinity, NaN, copy it to the second double of the result and
  you're done

- otherwise, figure out which of the numbers you're adding has the
  biggest exponent.  Take a copy of the result obtained so far into a
  temporary say tmp and, from it, subtract the most significant double
  of each of the input operands.  Now this temporary contains the
  negated value of whatever bits there were in the high doubles that
  had to be thrown away due to rounding or limited precision from the
  high part of the result.

- assign to the low part of the result the result of adding the
  negated temporary, the low doubles of the input operands, highest
  exponent first

in pseudo-code (untested):

long double tfadd (long double in1, long double in2) {
  long double result;
  double tmp;

  result.high = in1.high + in2.high;

  if (isinf (result.high) || isnan (result.high))
    { result.low = result.high; return result; }

  result.low = in1.low + in2.low;

  if (in1.exponent < in2.exponent)
    result.low -= result.high - in1.high - in2.high;
  else
    result.low -= result.high - in2.high - in1.high;

  return result;
}


For negation, just negate each double individually.

For subtraction, negate the second operand and add.

For widening multiplication (double x double -> long double):

- multiply the doubles into the high part of the result

- if the result is infinity, a nan or zero, copy the high part of the
  low part and stop.

- divide the result by one of the input doubles into tmp1

- subtract the other input double from tmp1 into tmp2

- multiply tmp2 by the other input double into the low part of the
  result

long double tfdfmult (double in1, double in2) {
  long double result;

  result.high = in1 * in2;

  if (isinf (result.high) || isnan (result.high) || result.high == 0.0)
    result.low = result.high;
  else
    result.low = (result.high / in1 - in2) * in1;

  return result;
}


For multiplication of two long doubles:

- widen-multiply each pair of doubles, and add the results from least
  to most significant.

long double tfmult (long double in1, long double in2) {
  if (exponent (in1) < exponent (in2))
    swap (in1, in2);

  return tfdfmult (in1.low, in2.low)  + tfdfmult (in1.low, in2.high)
       + tfdfmult (in1.high, in2.low) + tfdfmult (in1.high, in2.high);
}


To invert: ???  this is the missing bit :-(

To divide, invert the divisor and multiply.


For the numerical experts, are there any flaws in the procedures
above?  Any ideas on how to implement inversion without losing
(much) precision?

-- 
Alexandre Oliva   Enjoy Guarana', see http://www.ic.unicamp.br/~oliva/
Red Hat GCC Developer                 aoliva@{redhat.com, gcc.gnu.org}
CS PhD student at IC-Unicamp        oliva@{lsd.ic.unicamp.br, gnu.org}
Free Software Evangelist                Professional serial bug killer