Bad rounding in flat Modelica output

[casella]

Description

Flat Modelica output of constant binding values sometimes leads to bad rounding format

Steps to Reproduce

Flattening this model

model TestBadRounding
  parameter Real p = -0.94;
end TestBadRounding;

produces this flat output

class TestBadRounding
  parameter Real p = -0.9399999999999999;
end TestBadRounding;

Expected Behavior

Flat output should preferably be

class TestBadRounding
  parameter Real p = -0.94
end TestBadRounding;

I'm not an expert here, but I guess this is a classical CS problem 😄

Reported by @laviniaghilardi

[perost]

-0.94 can't be represented exactly as a double precision floating point number, the closest seems to be -0.939999999999999946... So in that regard it's correctly rounded since the next digit would be a 4.

But round trip serialization of floating point numbers is really hard to do correctly, enough so that it's been called C++17's Final Boss. We currently store floating point numbers as strings in the absyn structure in order to avoid these kind of issues with e.g. the diff algorithm, but the frontend converts them to floating point in order to be able to do anything with them. And storing them as strings wouldn't help if you write e.g. -0.84 - 0.1 instead.

One option would be to use Ryu, the library mentioned in the "C++17's Final Boss" video. That's something we would need to discuss though.

[perost]

@adrpo @sjoelund

[sjoelund]

If we would want to keep the representation the same, we would need to use some arbitrary precision library during the flattening, I suppose.

The output is going to parse to the same value if double precision is used. If you remove a single 9 from that, you will get a different result.

Suppose the code actually was -0.939999999999999946 - would you then expect -0.94 to be printed? That's what that Ryu library would print (I believe). If we would change to something like that library, I suppose it could potentially parse to a different number if the precision of the other tool was greater. I'm not sure any solution is preferred :)

[sjoelund]

Smaller output files are nicer, I suppose. So it might be "better"

[perost]

I wonder if it's actually a good idea to always print the same numbers as the user gave though. E.g. -0.94 can't be represented exactly in IEEE 754 double precision, so pretending that it can might not actually be in the best interest of the users. This might be a case of "ignorance is bliss" though.

[casella]

Arbitrary precision seems overkill to me. After all, Modelica is used for physical modelling, there's no way you can get physical parameters right with 16 decimal digits. The best known physical constant is the electron's magnetic momentum, which is known with an error of 7.6e-13, well above machine-epsilon for double precision, i.e. 2.2e-16. In most cases you get far, far less, like 1e-6, 1e-4 or even less.

The case where we'd need more precision is actually the one in the application where we spotted the problem, namely coefficients of polynomials, because of the well-known Wilkinson's effect, which makes the results of high-order polynomials very sensitive to small errors in the coefficients of the higher-order terms.

Using Ryu is one option. Another option, if possible, is to round the number to the 15-th decimal digit, instead of the 16-th, as it currently happens. We would lose a little bit of precision (5e-16), but that's no big deal. On the other hand, we would have the guarantee of getting back round number representations. I guess using "%.15g" as option for printf would do.

[casella]

I wonder if it's actually a good idea to always print the same numbers as the user gave though. E.g. -0.94 can't be represented exactly in IEEE 754 double precision, so pretending that it can might not actually be in the best interest of the users.

Remember that our users are engineers, not mathematicians or CS nerds. If you type 0.94 and you get 0.939999999999 you don't learn anything, you just get annoyed.

[sjoelund]

If you type 0.94 and you get 0.939999999999 you don't learn anything, you just get annoyed.

So what if they type 0.939999999999 and get 0.94?

[casella]

They would never type 0.939999999999 in the first place 😄

[perost]

The case where we'd need more precision is actually the one in the application where we spotted the problem, namely coefficients of polynomials, because of the well-known Wilkinson's effect, which makes the results of high-order polynomials very sensitive to small errors in the coefficients of the higher-order terms.

This is only about how the numbers are printed in the flat model though. If you write e.g. -0.94 it will still become 0.939999999999... in the simulation regardless of what print function we use, so printing -0.94 would hide the issue to the user if the precision actually matters.

[casella]

Well, I remember doing this with a pocket calculator in the early '80s. I did 1 divided by three, then multiplied by three and got 0.99999999. I was obviously annoyed. More recent calculators don't do this. Maybe because they are cheating, but what the hell...

[perost]

Well, I remember doing this with a pocket calculator in the early '80s. I did 1 divided by three, then multiplied by three and got 0.99999999. I was obviously annoyed. More recent calculators don't do this. Maybe because they are cheating, but what the hell...

Decent calculators use arbitrary precision numbers, since speed isn't usually a major concern for a calculator.

[sjoelund]

Using Ryu is one option. Another option, if possible, is to round the number to the 15-th decimal digit, instead of the 16-th, as it currently happens. We would lose a little bit of precision (5e-16), but that's no big deal. On the other hand, we would have the guarantee of getting back round number representations. I guess using "%.15g" as option for printf would do.

%.15g will have the same problem for other numbers. Edit: And for some things like pi you get a different number when parsing it again (like for example in C-code).

[sjoelund]

They would never type 0.939999999999 in the first place

They might put in pi though.

[casella]

The case where we'd need more precision is actually the one in the application where we spotted the problem, namely coefficients of polynomials, because of the well-known Wilkinson's effect, which makes the results of high-order polynomials very sensitive to small errors in the coefficients of the higher-order terms.

This is only about how the numbers are printed in the flat model though. If you write e.g. -0.94 it will still become 0.939999999999... in the simulation regardless of what print function we use, so printing -0.94 would hide the issue to the user if the precision actually matters.

What happens in the simulation is not really relevant, the result with two values of the parameter differing by 1ppq (1 part per quadrillion) will be the same from all practical purposes.

The point is: what is the purpose of the flat output? At this point, it's mostly for humans to inspect the result of flattening. From this point of view, those roundings would be quite annoying, while the additional "precision" they bring is absolutely irrelevant in practice.

[casella]

%.15g will have the same problem for other numbers.

I don't think so. The approximation of IEEE 754 double precision is one part in 2^52, where 52 is the number of bits in the mantissa. This corresponds to 2.2e-16. Hence, if you round to the nearest number with 15 digits of precision, decimal numbers with originally few decimal digits will be rounded to the nearest number with 15 digits, which will be, e.g. 0.9400000000000, and then the trailing zeros will of course be chopped out.

Edit: And for some things like pi you get a different number when parsing it again (like for example in C-code).

True. But in practice, you won't type 3.141592653589793, you'll type Modelica.Constants.pi, and then you'll get whatever rounded approximation of pi you get, and of course it will be fine

[sjoelund]

I don't think so. The approximation of IEEE 754 double precision is one part in 2^52, where 52 is the number of bits in the mantissa. This corresponds to 2.2e-16. Hence, if you round to the nearest number with 15 digits of precision, decimal numbers with originally few decimal digits will be rounded to the nearest number with 15 digits, which will be, e.g. 0.9400000000000, and then the trailing zeros will of course be chopped out.

Then the following should give you 0, right? It won't.

3.141592653589793 - 3.14159265358979;

The reason is that you assumed because the mantissa is stored in 52 bits, that gives you 52 bits of precision. Actually, you get 53 bits. So we need to use %.16g (or some library guaranteeing that the number will parse into the same double precision float). %.15g just won't do.

[sjoelund]

True. But in practice, you won't type 3.141592653589793, you'll type Modelica.Constants.pi, and then you'll get whatever rounded approximation of pi you get, and of course it will be fine

We use Flat Modelica for code generation and the same routines are used to generate C code. It must parse back to exactly the same representation.

[casella]

Then the following should give you 0, right? It won't.

3.141592653589793 - 3.14159265358979;

Sorry @sjoelund, this is not relevant. Let's try to stay focused on the issue.

We have a (small) problem: some decimal literal values, that look nice in the source code (e.g. -0.94), but look very ugly in the flat output (e.g. -0.939999999999999), because they don't have an exact FP representation.

My suggestion for the flat output is the following: display the maximum number of decimal digits that avoid bad rounding of decimal numbers that do not have an exact FP representation, while minimizing the loss of information compared to the binary FP representation.

My understanding is that 15 decimal digits bring a very small loss of information w.r.t. the binary representation (2-3 bits at most), which is not relevant in general because we never have such precision in the phyiscal value of parameters, while guaranteeing that the FP representation of -0.94 will be displayed as such, and not as -0.93999999999999999

In fact, if you are worried about loss of precision, we could make the max number of decimal digits an option of the flat modelica output. We could set it to 16 digits by default, then one could set it to whatever value he/she likes

Do we agree on that?

[casella]

We use Flat Modelica for code generation and the same routines are used to generate C code. It must parse back to exactly the same representation.

I understand we use the same internal representation (AST or something), but we are not using the same textual output that I get when I click on Instantiate Model in OMC directly to generate the C code. And even if we did, there is no reason why we couldn't have different representations of constants for these two outputs, that have different purposes.

Do I miss something?

[sjoelund]

We use Flat Modelica for code generation and the same routines are used to generate C code. It must parse back to exactly the same representation.

I understand we use the same internal representation (AST or something), but we are not using the same textual output that I get when I click on Instantiate Model in OMC directly to generate the C code. And even if we did, there is no reason why we couldn't have different representations of constants for these two outputs, that have different purposes.

Do I miss something?

We do when we use Flat Modelica as an intermediate stage to send between Modelica tools. And yes, we could have different outputs, but why? You'll just get a lot of people asking us why we're not using the full precision.

Using the Ryu library would probably be OK. Always output something that makes all bits of double precision the same. And if it has a shorter representation, output that.

[casella]

We do when we use Flat Modelica as an intermediate stage to send between Modelica tools. And yes, we could have different outputs, but why?

To avoid the ugly -0.93999999999 output 😃. That's a user requirement I am getting right now.

If there are other ways to avoid that, e.g. Ryu, I'm fine with it.

[sjoelund]

Well, that library should print a nice -0.94 which would also parse into exactly the same double precision float. So it should be OK also for the code generation, etc. So no different outputs needed, I hope.

[phannebohm]

Well, I remember doing this with a pocket calculator in the early '80s. I did 1 divided by three, then multiplied by three and got 0.99999999. I was obviously annoyed. More recent calculators don't do this. Maybe because they are cheating, but what the hell...

Calculators do fascinating things these days, like finding rational multiples of pi where they shouldn't.

[casella]

I love the guy, he's so incredibly witty :)

[adrpo]

With PR #11168 we get:

adrpo33@ida-0030 MINGW64 /c/home/adrpo33/dev/OMTesting
# omc TestBadRounding.mo
class TestBadRounding
  parameter Real p = -9.4E-1;
end TestBadRounding;

Seems better than the current -0.9399999999999999 but is not as the original.

[perost]

With PR #11168 we get:

adrpo33@ida-0030 MINGW64 /c/home/adrpo33/dev/OMTesting
# omc TestBadRounding.mo
class TestBadRounding
  parameter Real p = -9.4E-1;
end TestBadRounding;

Seems better than the current -0.9399999999999999 but is not as the original.

It seems like ryu doesn't generate the shortest fixed notation, see #199. std::to_chars in C++17 does and seem to be exactly what we want (it returns -0.94 in this case), but it wasn't implemented for floating point until GCC 11.

[casella]

The whole thing is a bit confusing to me.

• The library's README mentions "Ryu's printf" but there is no such function as printf in the library, only d2s, d2fixed, d2exp, f2s, see ryu.h
• On the other hand, the main claim of the README is that "This project contains routines to convert IEEE-754 floating-point numbers to decimal strings using shortest, fixed %f, and scientific %e formatting", so I guess d2s ad f2s are shortest, d2fixed is %f and d2exp is %e. I would have written that in the ryu.h file for clarity, explaining what each function does, but never mind. Apparently, documenting source code properly is still considered a bit uncool 😅
• The comment in #199 says "This repo doesn't generate shortest fixed notation". What? Isn't that the main claim of the library (".. using shortest formatting").

Apparently, "shortest" means "shortest exponential". I checked the paper, the algorithm for compute_shortest returns a decimal number d_0 and an exponent e_0. I also checked the the d2s_test.cc, it seems from the arguments of the ASSERT_D2S functions that the output format of d2s is the following:

• a minus sign for negative numbers
• one digit
• a decimal point and more digits if required to reproduce the same IEEE 754 bits as the original double-precision value
• the letter E
• the exponent of the power of 10 with a minus sign if negative, without any sign if positive

If that is the case, we can implement the function ourselves (sorry for the ugly Modelica-ish pseudo-code, my C programming is very rusty so I don't dare doing that in C)

str1 = d2s(r)
exp = the integer number following "E" in str1
sign = -1 if the first character in str1 is "-", else 1
digits = str1 stripped of the leading "-" (if present) and of all the characters from "E" onwards
ndec = number of characters in digits after the "." if "." is present, otherwise 0
str2 = ""
if sign < 0
  add "-" to str2
end if
if exp = 0
  add digits to str2
elseif exp > 0
  add digits[1] to str2
  add digits[3:min(ndec+2,exp+2)] to str2
  if exp > ndec
    for i in 1:(exp-ndec)
      add "0" to str2
    end for
  elseif exp < ndec
    add "." to str2
    add digits[exp+3:end] to str2
  end if
elseif exp < 0 
  add "0." to str2
  for i in 1:(-exp-1)
    add "0" to str2
  end for
  add digits[1] to str2
  add digits[3:end] to str2
end if

if len(str2) <= len(str1) // in case of tie, non-exponential wins
  return str2
else
  return str1
end if

[casella]

With this setup, we should get the following results:

1 -> 1
10 -> 10
100 -> 100
110-> 110
1000 -> 1e3
1e1 -> 10
1e2 -> 100
1e3 -> 1e3

In fact, I don't really like 1000 being represented as 1e3, I guess there's nothing wrong at keeping it as 1000 even if it is one character longer; scientific notation is nice to avoid unnecessary long numbers with a lot (-> more than three) of trailing or leading zeros. Three zeros are still ok. So, I would recommend to change the last lines to

if len(str2) <= len(str1) + 1 // slightly favour non-exponential form
  return str2
else
  return str1
end if

Then, we should check that we get the following values in return from MetaModelica

1 -> 1
10 -> 10
100 -> 100
110-> 110
1000 -> 1000
10000 -> 1e4
1e1 -> 10
1e2 -> 100
1e3 -> 1000
1e4 -> 1e4
1100 -> 1100
1110 -> 1110
123.4 -> 123.4
123.45 -> 123.45
123.456 -> 123.456
123.4567 -> 123.4567
123.45678 -> 123.45678
1.34 -> 1.34
0.94 -> 0.94
0.094 -> 0.094
0.0094 -> 0.0094
0.00094 -> 0.00094
0.000094 -> 9.4e-4
9.4e-3 -> 0.0094
9.4e-4 -> 0.00094
9.4e-5 -> 9.4e-5
1e1 -> 10
1e2 -> 100
1e3 -> 1000
1.1e3 -> 1100
1e4 -> 1e4
1.1e4 -> 11000

as well as the corresponding negative numbers.

Of course it would be really nice if we could keep the original form of the literal constants (exponential vs. fixed point), but I guess this is not really possible because that information is lost when parsing literals, and not even available when computing generic expressions. I think what we'll get is anyway a quite good compromise.

[casella]

@perost could you take care of translating my pseudo-code in plain C, so we can try it out? The way I wrote it should be pretty straigtforward, using *ptr++ statements.

[adrpo]

Seems they already have some patched version that might do what we want: https://github.com/ulfjack/ryu/issues/211#issuecomment-1126691741

[perost]

There's also double_to_fd128 in ryu that returns a struct with the mantissa and exponent, which would probably be a better starting point if we want to implement our own formatting.

[bilderbuchi]

std::to_chars in C++17 does and seem to be exactly what we want (it returns -0.94 in this case), but it wasn't implemented for floating point until GCC 11.

It seems like this could be the solution? Adapting the code from this comment it passes Francesco's list of test cases in https://github.com/OpenModelica/OpenModelica/issues/7465#issuecomment-1716790531, with the exception that the switchover to exponential notation happens at 1e5, not 1e4: https://gcc.godbolt.org/z/xTzr9xKEW

[bilderbuchi]

^ also, this seems to be effectively the same thing that std::cout is doing (but the switchover happens at 1e6):

void print_shortest(const double val) {
    std::cout << val << std::endl;
}

[adrpo]

Ok, I will give it a try, let's see what we get in the testsuite. Currently there seems to be some issues with some tests when using ryu.

[perost]

std::to_chars in C++17 does and seem to be exactly what we want (it returns -0.94 in this case), but it wasn't implemented for floating point until GCC 11.

It seems like this could be the solution? Adapting the code from this comment it passes Francesco's list of test cases in #7465 (comment), with the exception that the switchover to exponential notation happens at 1e5, not 1e4: https://gcc.godbolt.org/z/xTzr9xKEW

The issue is that we support Linux distributions with GCC as old as 8.5.0, so we can't really use std::to_chars.

^ also, this seems to be effectively the same thing that std::cout is doing (but the switchover happens at 1e6):

void print_shortest(const double val) {
    std::cout << val << std::endl;
}

I'm not sure what mechanism std::cout uses, but it seems to have its own share of issues. It prints e.g. 0.9399999999 as 0.94 while std::to_chars returns 0.9399999999. As far as I know only std::to_chars have any round-trip guarantees in the C or C++ standard libraries.

[sjoelund]

This sounds like something we could accept to do differently for the older distributions On 14 Sept 2023, at 10:00, "Per Östlund" wrote: std::to_chars in C++17 does and seem to be exactly what we want (it returns -0.94 in this case), but it wasn't implemented for floating point until GCC 11. It seems like this could be the solution? Adapting the code from this comment https://github.com/ulfjack/ryu/issues/199#issuecomment-894598115 it passes Francesco's list of test cases in #7465 (comment) https://github.com/OpenModelica/OpenModelica/issues/7465#issuecomment-1716790531 , with the exception that the switchover to exponential notation happens at 1e5, not 1e4: https://gcc.godbolt.org/z/xTzr9xKEW https://gcc.godbolt.org/z/xTzr9xKEW The issue is that we support Linux distributions with GCC as old as 8.5.0, so we can't really use std::to_chars. ^ also, this seems to be effectively the same thing that std::cout is doing (but the switchover happens at 1e6): void print_shortest(const double val) { std::cout << val << std::endl; } I'm not sure what mechanism std::cout uses, but it seems to have its own share of issues. It prints e.g. 0.9399999999 as 0.94 while std::to_chars returns 0.9399999999. As far as I know only std::to_chars have any round-trip guarantees in the C or C++ standard libraries. — Reply to this email directly, view it on GitHub https://github.com/OpenModelica/OpenModelica/issues/7465#issuecomment-1718954628 , or unsubscribe https://github.com/notifications/unsubscribe-auth/AAXAODDIRD2FRTB2K6VV633X2K2R5ANCNFSM444YLIIA . You are receiving this because you were mentioned.Message ID: @.***>

[adrpo]

Ouch, looks like we need a new OMDev for this, as we only have gcc 10 in there :)

[adrpo]

I can push my changes but they won't compile on Windows :)

[mahge]

I am not sure if there is a reason for not considering it but another alternative is of course libfmt. It is what the C++20 std::format is based on. I am not sure which algorithm it uses internally but there is a possibility to adjust the formatting to some extent.

I have been using it for a long time and it is good for a lot more reasons than formatting floating point numbers. It is a library that we should definitely consider if we are going with C++ based formatting. It is available in every Linux distribution I have come across and is available in MSYS as well. It requires C++17 but that should be okay for us now.

I made some quick tests and you can see the result here https://godbolt.org/z/cPToEaWGf

[perost]

I am not sure if there is a reason for not considering it but another alternative is of course libfmt. It is what the C++20 std::format is based on. I am not sure which algorithm it uses internally but there is a possibility to adjust the formatting to some extent.

I have been using it for a long time and it is good for a lot more reasons than formatting floating point numbers. It is a library that we should definitely consider if we are going with C++ based formatting. It is available in every Linux distribution I have come across and is available in MSYS as well. It requires C++17 but that should be okay for us now.

I made some quick tests and you can see the result here https://godbolt.org/z/cPToEaWGf

libfmt seems to use a similar algorithm to ryu (Dragonbox in newer versions, Grisu in older), so that might also be a good alternative. And as you say it's a very useful library to have in general.

[mahge]

Ahh that's the one: https://github.com/jk-jeon/dragonbox.

I actually checked it when this ticket was first opened (I think I got there from libfmt actually). But then we went with ryu and I forgot about it completely 🤦‍♂️

[casella]

As I understand, ryu does an excellent job at figuring out the shortest possible decimal representation of the mantissa and exponent. That is the conceptually hard part of the problem (as explained in the companion paper). Turning that output into either decimal or exponential is a no-brainer.

Rather than throwing away ryu and trying with other libraries, which most likely have their own issues (don't work with all distributions/OSs, don't provide the output we want in some cases, require a new OMDev, or whatever), so that after some more work we're back to square one, I would recommend to stick to ryu and program the last step ourselves. It will probably require less time, and we'll get a better output.

There's also double_to_fd128 in ryu that returns a struct with the mantissa and exponent, which would probably be a better starting point if we want to implement our own formatting.

This can be a good idea, even though parsing (-)xxxEyy doesn't seem like a big deal 😅

[casella]

Regarding libmft, are we OK with these requirements?

[perost]

Regarding libmft, are we OK with these requirements?

Yes, very ok, it only requires some C++11 features while OM requires C++17. It should also be noted that fmt is a very popular C++ library (std::format in C++20 is based on fmt::format and more will be included in C++23), and as mentioned would also be generally useful.

[casella]

Can we give it a quick try (standalone, not in OMC) and see if we get something like this?

1 -> 1
10 -> 10
100 -> 100
110-> 110
1000 -> 1000
10000 -> 1e4
1e1 -> 10
1e2 -> 100
1e3 -> 1000
1e4 -> 1e4
1100 -> 1100
1110 -> 1110
123.4 -> 123.4
123.45 -> 123.45
123.456 -> 123.456
123.4567 -> 123.4567
123.45678 -> 123.45678
1.34 -> 1.34
0.94 -> 0.94
0.094 -> 0.094
0.0094 -> 0.0094
0.00094 -> 0.00094
0.000094 -> 9.4e-4
9.4e-3 -> 0.0094
9.4e-4 -> 0.00094
9.4e-5 -> 9.4e-5
1e1 -> 10
1e2 -> 100
1e3 -> 1000
1.1e3 -> 1100
1e4 -> 1e4
1.1e4 -> 11000

[mahge]

I made some quick tests and you can see the result here https://godbolt.org/z/cPToEaWGf

[perost]

I made some quick tests and you can see the result here https://godbolt.org/z/cPToEaWGf

I played around with it a bit and unfortunately found a major issue with fmt. With the default formatting it prints e.g. 940234234000000.0 as 940234234000000.0, while with g it instead prints 9.40234e+14 (std::to_chars gives 9.40234234e+14).

I can't find any way to get it to print the shortest representation without rounding the number. Maybe you have some clue?

[mahge]

You can specify the precision using the . specifier. In this case you can use something like

  fmt::print("{:.9g}\n", 987654321000000.0);

to get 9.87654321e+14. This formatting will correctly(?) give you 9.8765432e+14 for 987654320000000.0 (1 replaced by 0)

You can see the effect on the other cases here: https://godbolt.org/z/q6vvE1rn9

However, I am not sure if this solves the actual problem since it is probably not being as "smart" as we want it to be in figuring out the shortest representation. It is just doing as much precision as we asked.

To be honest I never really tried all its capabilities in printing floating point numbers. I just use it to format strings as replacement for std::cout. So there might be more nuances to it than I can explain but I will keep checking.

[casella]

I can't find any way to get it to print the shortest representation without rounding the number. Maybe you have some clue?

The problem with all these functions is that they won't return explicitly the shortest mantissa that preserves the original bitwise representation (which is exactly what ryu's d2s and double_to_fd128 do!), and rather try to second-guess what the end user would like to see, unless a specific number of significant digits is specified (which we don't want to do). This may or may not be what we actually need in OMC, depending on specific cases.

Unless some of these functions have some explicit option to print out the number with the minimal number of mantissa digits that preserve the bitwise representation (which would be handy, of course), I'm afraid we'll have to do this ourselves, as I proposed above, so we maintain full control over what we get.