inconsistent Float16 printing

[stevengj]

In Julia 0.3 and 0.4, I see:

julia> float16(44099)
float16(44100.0)

julia> a = Float16[44099]
1-element Array{Float16,1}:
 44096.0

julia> a[1]
float16(44100.0)

[eschnett]

Note that float16's epsilon for this number is 32. 44096 is divisible by 32, 44100 is not -- that is, 44096 is the correct output. Maybe 44100 is produced by rounding 44096 to 3 significant digits?

Or maybe it should be 4.41e4?

[andreasnoack]

I get

julia> float16(44099)
float16(44096.0)

with

julia> versioninfo()
Julia Version 0.3.3-pre+78
Commit 75b347b* (2014-11-19 10:23 UTC)
Platform Info:
  System: Darwin (x86_64-apple-darwin14.0.0)
  CPU: Intel(R) Core(TM) i7-4870HQ CPU @ 2.50GHz
  WORD_SIZE: 64
  BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY Haswell)
  LAPACK: libopenblas
  LIBM: libopenlibm
  LLVM: libLLVM-3.3

[quinnj]

This is grisu-related. See #6608 and the discussion here.

The short answer is because regular show uses the grisu shortest representation whereas the array show code uses a fixed precision. shortest here is 41100.0 because it's the fewest number of non-zero digits needed to be input back in to the same number, as can be seen by

In  [3]: bits(float16(44099))

Out [3]: "0111100101100010"

In  [4]: bits(float16(44100.0))

Out [4]: "0111100101100010"

Sorry I haven't done more about this, I'm just pretty slammed for the next 3 weeks or so and it's not entirely obvious what should be done here (though #6608 has some excellent recent suggestions to address the array printing).

[eschnett]

The number 44100.0 carries no indication that it has only 3 significant digits. Either it needs to be printed as 4.41e4, or as 44096.0 (or maybe as 44100.0±50.0).

[jiahao]

@eschnett Your comments imply that floating point numbers carry some intrinsic imprecision with them. On the contrary, this is literally the No. 1 item on the very top of Kahan's list of floating point misconceptions (e.g. here). The notion of significant digits simply does not apply here since floating point numbers are exact.

The real issue here is that there can be multiple decimal representations for a given floating point number, analogous to how there are two decimal representations for the real number 1 = 0.999... . In the OP, 44100.0 and 44096.0 are equally valid decimal representations for the Float16 number +1.0101100010₂ * 2 ^ 11110₂, the issue is that we don't pick the same decimal representation consistently.

[eschnett]

@jiahao No, I didn't imply this. When I spoke of significant digits, I spoke of decimal digits in the output, not the floating point number.

Each floating point number is equal to exactly one rational number, in this case 44096. There is only one valid decimal representation for this number, namely 44096.

One can note that 3 decimal digits suffice to represent a Float16 number (as you describe), and thus round 44096 to 3 decimal digits. This gives 44100. However, when one rounds, one needs to indicate this. Thus one needs to write either 4.41e4, or 44100±50.

Both 44096 and 44100, when rounded to the nearest Float16 number, are equal, since both are rounded to 44096. But calling 44100 a valid representation of 44096 is only correct if you state at the same time that you are rounding to 3 decimal digits.

Compare the output of @printf "%f":

julia> @printf("%f", float32(2^40))
1099511627776.000000
julia> @printf("%f", float16(44096))
44096.000000

These numbers (2^40 and 44096) can be represented exactly as Float32 and Float16, respectively, and they are not rounded when output. Float32 has about 7 decimal digits of precision, but Julia (and other languages) output more digits, giving the exact number.

[jiahao]

But calling 44100 a valid representation of 44096 is only correct if you state at the same time that you are rounding to 3 decimal digits.

Agreed, but that is not what I said at all. What I said was that both 44100.0 and 44096.0 are correct decimal representations of +1.0101100010₂ * 2 ^ 11110₂.

[stevengj]

I think that Kahan would say that +1.0101100010₂ * 2 ^ 11110₂ represents the integer 44096, exactly, and only that number. In standard notation, "44100.0" is not read as a representation of 44096, it is read as a representation of the integer 44100 (which is not in the Float16 subset of ℚ).

Of course, in general we cannot print an exact representation of a binary floating-point number x in decimal, and so we are forced to print a value that merely rounds to x: I would say that it reproduces x, not that it represents x. However, when the floating-point value represents a number that can be expressed exactly in decimal with the requested number of digits, it would be nice to print that.

cc: @alanedelman for his gloss on the holy writ.

[stevengj]

In any case, it would be nice to be consistent.

[eschnett]

+1.0101100010₂ * 2 ^ 11110₂ is equal to +10101100010.0₂ * 2 ^ 101000₂ (shifting the point), which is equal to 1376 * 2^5 (converting to decimal), which is equal to 44096. There is no rounding involved. One may need to round when converting to Float16, but not when converting from Float16 to a decimal representation.

I general, we can print an exact representation. Every Float16 has an exact decimal representation. The question is just how many digits this requires. If we want to round to fewer digits, so be it, but that's then a question of rounding a decimal representation.

[eschnett]

@stevengj I consider this a bug. 44100.0 has three zeros that indicate a precision that is not there, since the number was rounded. This is misleading. Other languages don't behave this way:

#include <math.h>
#include <stdio.h>
int main(int argc, char** argv)
{
  printf("%f\n", pow(2.0, 100.0));
  return 0;
}

This outputs 1267650600228229401496703205376.000000, which is the exact result. There is no rounding to 17 digits, which would suffice to reproduce the number.

[quinnj]

@eschnett this is certainly a bug when considering non-Julia input/output; that is, if 44100.0 is written out to a file and read back in another programming language or something, then obviously this presents problems. For simple input/output within julia, float16(4110.0) is the same input/output, which I think is what others are arguing isn't technically wrong.

I think we all agree though that this needs to be improved. I'll have some time in a few weeks, but I'm afraid I can't do much to help until then.

[stevengj]

Right, we can print an exact decimal repr, but not necessarily with the requested number of digits. But in the case of show(x::Float16) where x is an integer. I don't understand the motivation for not printing all the decimal digits. Can someone more familiar with grisu enlighten me?

[JeffBezanson]

I still get the Float16(44100.0) output.

[vtjnash]

yeah, i had to revert part of it to fix some tests. i'll get back to this soon and close this again

[JeffBezanson]

I'm not sure this is fixed. Now I get

julia> Float16(44099)
Float16(4.41e4)

julia> Float16[44099]
1-element Array{Float16,1}:
 44096.0

julia> print_shortest(STDOUT,  Float16(44099))
44100

julia> Int(Float16(44099))
44096

The first output is definitely better, but it's not clear why it needs to be different from the output in the array. The output of print_shortest seems pretty suboptimal.

[vtjnash]

They start down different code paths that can't readily be rejoined at the end. The first attempts to compute the minimal representation, then decide on the format of the digits. The second attempts to print exactly (no more than) N digits. This causes an issue for Float16 (but not more accurate types), since it has less than 6 digits of accuracy (float has 6, double has 15)

[vtjnash]

I'm not sure this is fixed. Now I get

those are all correct outputs now. suboptimal is perhaps just a matter of taste, since 4.41e4 would have required printing more characters, it is not the shortest representation. (showcompact is used in the array output to request printing of 5 digits).

[JeffBezanson]

Should we open a separate issue for print_shortest? It seems silly to print 44100 when the exact value 44096 is the same number of characters.

[vtjnash]

see #9072 (comment)

[JeffBezanson]

Fair enough. We should edit the help text to say that it minimizes nonzeros too.

[vtjnash]

"it prints the minimum number of consecutive non-zero digits in the minimum number of characters"?

that at least is the basic description of the algorithm