chapel-lang / chapel

a Productive Parallel Programming Language
https://chapel-lang.org
Other
1.8k stars 422 forks source link

[Bug]: Hexadecimal Floating Point Output Format #25014

Open damianmoz opened 6 months ago

damianmoz commented 6 months ago

It might be a feature but not a nice one.

What does 0x1.0p-1074, a real(64), print out under %30.18xer (or ....xr)

0.000000000000100000p-1022

I might have got those leading zeros wrong. That format is not very useful.

Or am I doing something really stupid. Thanks.

lydia-duncan commented 6 months ago

Hi Damian! Thanks for reporting. I certainly wouldn't be surprised if we failed to properly test the overlap of all of those format specifiers. But I'm having trouble determining exactly what the correct result should be (in part because I don't think I've ever written exponential notation for anything other than base 10, so I'm having trouble convincing WolframAlpha to tell me what the conversion is).

I think there's a couple of things you could be saying here (where some combination is more than reasonable):

Is one or more of those an accurate representation of your issue? Or have I missed the point entirely?

damianmoz commented 6 months ago

Somewhat. But that is because I wrote this in a rush before I went to bed.

It is not urgent.I will elaborate on the details more clearly later this week.

You do not need to know how to do the conversion. That part of the existing code is correct. However, that code does not correctly handle a subnormal floating point number. I think the fix is straightforward but at 11pm when I posted the issue, I was simply interested in noting the issue By the way, I was trying with '%12.4xr originally and that works flawlessly if it is fed floating point numbers which are what IEEE 754 standard calls normal numbers, i.e. any real(64) >= 2^-1022 or a real(32) >= 2^-126. Anyway, it is before 4am here and I need some more sleep.

lydia-duncan commented 6 months ago

(No need to respond until more normal hours) After a brief Wikipedia dive, I think what you're saying is that it should display as 1.0p-1035 in decimal, and it's not right now because of how the floating point is stored in memory?

damianmoz commented 6 months ago

Note that currently, the code works perfectly for the case where the datum to be written is a normal IEEE 754 number. Here normal means the datum lies in the range 2^-1022..<inf. It does not mean abnormal but rather subnormal.

It should display the tiniest real(64) as 2^-1074 as 1..0p-1074., or something like 2^-1035 as 1.0p-1035, or the general case of M * 2^N where 1 <= M < 2 as M in hexadecimal and N as the exponent of 2 in decimal (which it currently does correctly when N>=-1022). For example, it should display

1.5 * 2^-1030 as 1.8p-1030
1.65625 * 2^-1040 as 1.a8p-1040 // (Note that 1.65625 = 1.625 + 0.03125)
1.5 * 2^-1022 as 1.8p-1022 // this currently works perfectly

How IEEE 754 stores things in memory is what you start with. It is the baseline. It is not the problem.

If you point me at the formatted I/O code, I can probably figure out the mistake.

I am guessing the current code simply takes the floating point number in question, looks at the biased exponent and the significand, removes the bias from the exponent and starts formatting. Certainly that would give you the result you now see.

What was omitted is to upscale the significand and downscale by the exponent if the biased exponent is zero, the former by a multiplication and the latter (because it is an exponent) by the corresponding subtraction. At the moment this is not done when that exponent is zero, which results in a value of M which lies in the range 2^-52<=M<1. The current code correctly displays that value of M correctly but that is not what '%x' format is all about. The leading digit should be a 1 followed by a decimal place. Try using '%a' format in C and look at the result.

I could not find the offending code in a quick perusal but I was probably looking in the wrong place.

damianmoz commented 6 months ago

Just looking at the documentation and it says

%xer
    hexadecimal number using p to mark exponent e.g. 6c.3f7p-2a

Luckily this is wrong. This outputs a garden variety Chapel hexadecimal real but without the leading '0x' that you would see in a Chapel literal, i.e.

const x = 0x1.6c3f7p-19;

This will output as 1.6c3f7p-19.

This is an exponential notation with 'e' replaced by 'p', a base of 2 mandated for the exponent which itself is expressed in decimal and the digits in the coefficient (or significand or more loosely mantissa) expressed in hexadecimal. The output will always have a leading 1 (as far as I know but I am no guru on Chapel's formatted I/O.

As another example, the variable x following

const x = 0x0.6c3f7p-2;

writef("%xr\n", x);
assert(x == 0x1.f8p-5);

will output under '%xr' control as '1.fb8p-5 (which is correct).

The '%xr' format currently prints what C would print as '%a' format. That equivalence should go in your table.

damianmoz commented 6 months ago

I also think that mentioning %xer is confusing when %xr works fine. The e seems redundant to me.

damianmoz commented 6 months ago

This hexadecimal significand and decimal exponent format for '%x' comes from IEEE 754 clause 5.12.3. Technically you should also allow a '%X' format to use a 'P'. But I have never done so in my life.

lydia-duncan commented 6 months ago

If you point me at the formatted I/O code, I can probably figure out the mistake.

I'm not as familiar with this part of our implementation, so this is more of a guess and @mppf (who I think wrote it, albeit a decade ago) may have better pointers. But I believe it's handled by the code that starts in https://github.com/chapel-lang/chapel/blob/main/runtime/src/qio/qio_formatted.c#L4704? That code is pretty extensive, so it may be handled in one of the calls it makes or maybe in something else entirely in the call chain from https://github.com/chapel-lang/chapel/blob/main/modules/standard/IO.chpl#L11385

Just looking at the documentation and it says

%xer
    hexadecimal number using p to mark exponent e.g. 6c.3f7p-2a

Luckily this is wrong. This outputs a garden variety Chapel hexadecimal real but without the leading '0x' that you would see in a Chapel literal, i.e.

const x = 0x1.6c3f7p-19;

This will output as 1.6c3f7p-19.

Huh. Are you saying that we output in hexadecimal format until we get to subnormal numbers? That definitely sounds like something we shouldn't be doing imo - it seems to me like it should at the very least be consistent for all numbers in a particular base, and I'd like it to match the base that was specified by the format specifier.

The '%xr' format currently prints what C would print as '%a' format. That equivalence should go in your table.

That's good to know, thanks! It may be a good idea to make a separate issue reporting that so it doesn't get lost in the details of this one.

I also think that mentioning %xer is confusing when %xr works fine. The e seems redundant to me.

This made me curious, so I wrote a quick little code to check our behavior:

use IO;

var a = 0x34;
var b = 0x1.f8p-5;
writef("%r\n", a);
writef("%xr\n", a);
writef("%xer\n", a);
writef("%r\n", b);
writef("%xr\n", b);
writef("%xer\n", b);

This didn't print how I was expecting. My thinking was that 0x34 should get printed as 0x34 or 3.4p+2 in hexadecimal. But instead this prints as:

52
1.ap+5
1.ap+5
0.0615234
1.f8p-5
1.f8p-5

Am I showing my ignorance here? The number is already using exponential notation printed the same as it was written, why wouldn't it display the same? Does exponential notation use a different base in different base systems? I suppose that makes sense when I think about it, it just isn't something I've had to know before.

It seems like the exponential notation is what you wanted, is that always true or are there real use cases where it would be important to print the hexadecimal number exactly as it was originally written?

mppf commented 6 months ago

I'm not as familiar with this part of our implementation, so this is more of a guess and @mppf (who I think wrote it, albeit a decade ago) may have better pointers. But I believe it's handled by the code that starts in https://github.com/chapel-lang/chapel/blob/main/runtime/src/qio/qio_formatted.c#L4704? That code is pretty extensive, so it may be handled in one of the calls it makes or maybe in something else entirely in the call chain from https://github.com/chapel-lang/chapel/blob/main/modules/standard/IO.chpl#L11385

The C code you pointed to is just parsing the format strings. This section actually does hexadecimal float output:

https://github.com/chapel-lang/chapel/blob/e6e0f2fecfc667eed4d9ea020eb5bdc6fe35b32a/runtime/src/qio/qio_formatted.c#L3298-L3314

damianmoz commented 6 months ago

I (re)learnt something today.

It seems that printf behaviour for a subnormal double and '%a' (or accurate) format is to output such subnormal numbers with an exponent of -1022 and a fractional hexadecimal exponent where the leading digit is zero. I had forgotten this.

However, because qio works internally with a real(64) and not a real(32), when it is given a real(32), its subnormal behaviour is what I would expect, i.e. it does use a non-zero leading digit.

Someone should note that 'writef' is not generic.

I will submit an issue to correct the documentation with '%xer' and suggest it change to '%xr' which also works and in there will put some extra words about the details of the format being used (which somebody can cannibalise if they want).

My apologies for an issue which (in retrospect) is non-issue although it did show up the error in the documentation.

That said, at some stage, I would like re-implement my Simula67-like library which handles these such things almost properly (or at least my definition of the word). Its uses routines which are all called out? which are maybe not as convenient as printf/writef style but avoid all the mess with variable argument lists and are fully generic. I need to find the source and fix the rounding bugs which I never fixed 20 years ago. Maybe something to do to keep me off the streets in Christmas 2026!!

damianmoz commented 6 months ago

Before I write anything, is the syntax for what

%[@][+~][<^>]{xX}[e]r

where things in square brackets are optional, things in squiggly brackets are mandatory, and @ means prepend 0x (or 0X), + means use a plus sign if positive, ~ (tilde) is my poor attempt to imply a space and says use a space if positive, the <^> mean left, center and right justified, x says (use 0x if requested by @ and) use p before the exponent, X says (use 0X if requested by @ and) use P before the exponent, and I think the 'e' is redundant but you can insert it if you like.

For what it is worth, how often do people use I centered formatting? And when they do, for what do people use it. I would like to introduce them to 'tbl'!

I suggest

Floating point values may also be specified using a hexadecimal format using the 0x prefix. In such cases, p is used to indicate the exponent (avoiding ambiguity with the hexadecimal digit e) and indicates a power of 2 rather than 10:

This is what IEEE calls a hexsequence. Strictly it is a hexadecimal floating constant optionally scaled by a decimal binary exponent. Actually C calls the entirety of that literal the hexadecimal floating constant and mandates the decimal binary exponent (and some people call that form, i.e. one which mandates the exponent begun with P or p, as P notation)

Anyway, what do people think about a description in

https://chapel-lang.org/docs/users-guide/base/basicValues.html

replacing those earlier words.

Floating point values may also be specified using a combination hexadecimal/decimal exponential format. It is begun with (optional sign) and '0x', is then followed by a hexadecimal digit sequence, optionally followed by a hexadecimal fraction (a hexadecimal point and a hexadecimal digit sequence), optionally followed by a 'p' (or 'P') and binary exponent (an (optionally signed) decimal digit sequence which indicates the power of 2 which scales the hexadecimal floating constant seen between '0x' and the 'p' (or the 0X and 'P')).

What a mouthful. Maybe somebody can abbreviate my waffle.

Anyway, it is this format (or should it be hexi[deci] notation) that writef outputs when fed a '%x' format control string (or readf inputs).

lydia-duncan commented 6 months ago

Before I write anything, is the syntax for what

%[@][+~][<^>]{xX}[e]r

where things in square brackets are optional, things in squiggly brackets are mandatory, and @ means prepend 0x (or 0X), + means use a plus sign if positive, ~ (tilde) is my poor attempt to imply a space and says use a space if positive, the <^> mean left, center and right justified, x says (use 0x if requested by @ and) use p before the exponent, X says (use 0X if requested by @ and) use P before the exponent, and I think the 'e' is redundant but you can insert it if you like.

Your understanding is correct, though {xX} is only necessary if you want it to print in hexadecimal, you can just write %r for other situations.

For what it is worth, how often do people use I centered formatting? And when they do, for what do people use it. I would like to introduce them to 'tbl'!

It was something we recently added, so I don't know that we have much information for usage in Chapel, but we got the idea from Python iirc.

Anyway, what do people think about a description in

https://chapel-lang.org/docs/users-guide/base/basicValues.html

replacing those earlier words.

Seems reasonable to me, though I do think we can probably iterate on it a bit more in the interest of clarity. And we'll still want to do the update you suggested about mentioning %xr is similar to %a in C.

damianmoz commented 6 months ago

I noticed in passing that what IEEE calls a Hex-Sequence format is what Fortran calls a Hex-Significand format which is arguably more accurate because the latter spells out that it is only the significand that is hexadecimal Sadly, neither format's name gives a hint as to what base the exponent is. An interesting observation.

I will think about and document the words we can use to describe '%xr' on the weekend and others can improve on them next weekl

damianmoz commented 6 months ago

Sorry. I bumped UPDATE COMMENT too early. There was a lot of cut and pasted rubbish in the first attempt.

It was mentioned above:

Your understanding is correct, though {xX} is only
necessary if you want it to print in ```hexadecimal,
you can just write %r for other situations.

I was only intending to tweak '%xr' I was not proposing to change anything else.

That said I find the description of '%6.4r' as

as with %r but padded on the left to
6 columns and with 4 significant digits

a big odd grammatically. Somebody else when they first read it thought it meant there were 6 blanks on the left but maybe they forgot their glasses.

It displays in a numeric field which is 6 columns wide
including an explicit or implicit sign (which by default
is blank if positive), with a total of 4 significant digits.
It is left padded with blanks.

Well, in this case 6-4(digits)-1(sign) = 1 blank if using a non-exponential format

Too verbose. But the current description is too brief.

It does need to be pointed out that Chapel's %w.nr format is unrelated to C's %w.pf or Fortran's Fw.d for that matter, or the 'g' format of either, or even Chapel's own %w.ndr. I am still not sure which I prefer. I really like Chapel's paradigm because it is easier to describe mathematically.