- if (e == 0) {
- /* Minimum exponent indicates zero or a subnormal number. The
- * subnormal exponent is a sentinel value that shouldn't be taken
- * literally, so we should fix that. If the number is actually zero
- * then the exponent won't matter much so don't bother checking.
- */
-
- e = 1;
- } else {
- /* It's a normal number. In which case there's an implicit bit which
- * we can now set.
- */
-
- t |= 0x00100000;
- }
-
- /* All that remains is to stuff the significant and exponent into a
- * floating point number. We'll have to do this in pieces, and we'll
- * lean on the floating-point machinery to do rounding correctly.
- */
- x = ldexp(t, e - 1043) + ldexp(lo, e - 1075);
- if (s) x = -x;
- }
-
- /* And we're done. */
- *x_out = x; return (0);
-}
-
-/* --- @f64_to_k64@ --- *
- *
- * Arguments: @double x@ = a floating-point number
- *
- * Returns: A 64-bit encoding of @x@.
- *
- * Use: Encodes @x@ as a `binary64' value. See `buf_putf64' for the
- * caveats.
- */
-
-static kludge64 f64_to_k64(double x)
-{
- kludge64 k;
- uint32 lo, hi, t;
- int e; double m;
-
- /* Some machinery before we start. */
-
- if (NANP(x)) {
- /* A NaN. */
- hi = 0x7ff80000; lo = 0;
- } else if (INFP(x)) {
- /* Positive or negative infinity. */
- hi = NEGP(x) ? 0xfff00000 : 0x7ff00000; lo = 0;
- } else if (x == 0) {
- /* Positive or negative zero. */
- hi = NEGP(x) ? 0x80000000 : 0; lo = 0;
- } else {
- /* A normal or subnormal number. Now we have to do some actual work. */
-
- /* Let's get the sign dealt with so we don't have to worry about it any
- * more.
- */
- if (!NEGP(x)) hi = 0;
- else { x = -x; hi = 0x80000000; }
-
- /* Now we start on the value. The first thing to do is to split off the
- * exponent. Our number will be %$m \cdot 2^e$%, with %$1/2 \le m < 1$%.
- */
- m = frexp(x, &e);
-
- /* If our number is too big, we'll round it to infinity. This will
- * happen if %$x \ge 2^{1024}$%, i.e., if %$e > 1024$%.
- */
- if (e > 1024)
- { hi |= 0x7ff00000; lo = 0; }
- else {
- /* Our number is sufficiently small that we can represent it at least
- * approximately (though maybe we'll have to flush it to zero). The
- * next step, then, is to pull the significand bits out.
- */
-
- /* Determine the correct exponent to store. We're not going to bias it
- * yet, but this is where we deal with subnormal numbers. Our number
- * is normal if %$x \ge 2^{-1022}$%, i.e., %$e > -1022$%. In this
- * case, there's an implicit bit which we'll clear. Otherwise, if it's
- * subnormal, we'll scale our floating-point number so that the
- * significand will look right when we extract it, and adjust the
- * exponent so that, when we're finally done, it will have the correct
- * sentinel value.
- */
- if (e > -1022) m -= 0.5;
- else { m = ldexp(m, 1021 + e); e = -1022; }
-
- /* Now we pull out the 53 bits of the significand. This will, in
- * general, leave a tail which we address through rounding. Scale it
- * up so that we end up with %$0 \le m' < 2$%; then we round up if
- * %$m > 1$%, or if %$m = 1$% and the low bit of the significand is
- * set.
- */
- t = ldexp(m, 21); m -= ldexp(t, -21);
- lo = ldexp(m, 53); m -= ldexp(lo, -53);
- m = ldexp(m, 54);
-
- /* Round the number if necessary. */
- if (lo&1 ? m >= 1.0 : m > 1)
- { lo = U32(lo + 1); if (!lo) t++; }
-
- /* Now we just put the pieces together. Note that our %$e$% is one
- * greater than it should be, because our implicit bit should have
- * been the unit bit not the 1/2 bit.
- */
- hi |= ((uint32)(e + 1022) << 20) | t;
- }
- }
-
- /* Convert to external format and go home. */
- SET64(k, hi, lo); return (k);
-}
-
-/*----- External functions ------------------------------------------------*/
-
-/* --- @buf_getf64{,l,b} --- *
~mdw / mLib / blobdiff