H->HMS conversion HP-15C vs. HP42S vs HP67

[Werner]

(09-06-2018 05:31 PM)Dieter Wrote:  EDIT AND UPDATE:
I now have an idea of what might be going on here: it looks like the usual binary roundoff problem.
This 15C simulator seems to use binary arithmetics as most computer hardware does. This means that 183,38 cannot be represented exactly as a binary value and its decimal equivalent may be 183,37999999999... While decomposing the input into hours, minutes and seconds this in turn is interpreted as 183 hours, 37 minutes and 99,99999... seconds. Et voilà: this indeed equals 183,6444444... hours.

Case closed ?-)
Dieter

Yet Free42 binary (of which I must be the only user?) gets it right ;-)
Werner

[Thomas Okken]

(09-07-2018 09:43 AM)Werner Wrote:  

(09-06-2018 05:31 PM)Dieter Wrote:  it looks like the usual binary roundoff problem.

Yet Free42 binary (of which I must be the only user?) gets it right ;-)

The approach I chose for ->HR for the binary version is robust but not very pretty: multiply x by 1e12, round to nearest integer, and then do everything using integer math. Plus, there is a special case in for |x| < 0.01, where the conversion reduces to a simple division by 0.36.

The idea was to make sure that it will never get anything wrong that the HP-42S gets right, but the simplistic approach I chose does mean a bit of precision is lost. If someone has a better algorithm lying around for IEEE-754 double precision binary, I would be interested.

[Albert Chan]

(09-07-2018 07:07 AM)Pekis Wrote:  And with my funny one-liners:
DEG->DMS: (90*X+INT(60*X)+100*INT(X))/250
DMS->DEG: (250*X-INT(100*X)-60*INT(X))/90

Hi, Pekis

Above work (mostly) in Excel, because Excel does some behind-the-scene rounding.

Example: if angle = 2° 1' , X = 2.01 (IEEE double of 2.0099999...)
DMS->DEG formula:

Python: (502.499999... - 200 - 120) / 90 = 2.0277777... (wrong!)
Excel : (502.5 - 201 - 120) / 90 = 2.0166666... (all terms rounded implicitly)

I would rather do explicit scaling / rounding, before conversion (see post #17):

X = 2.0099999 ... --> Y = 20100.00 (scaled and rounded to hundredth of second)

DMS->DEG: (Y - 40*int(Y/100) - 2400*int(Y/10000)) / 3600 = (20100 - 8040 - 4800) / 3600 ~ 2.0166666...

[AndiGer]

Retro15C has funny results ...
Any positive integer.38 yields in positive integer.6444...
0.38 yields in 0.6333...
-0.30 yields in - no not in -0.5 - but (+)0.1666...

[Gene]

Think I'm not using that one any longer.

[Siegfried]

(09-06-2018 10:37 PM)Thomas Okken Wrote:  I'm guessing there's a typo in his post, because it seems unlikely anyone would mistype the actual test case three times!

Thanks for all the testing and detailed replies.

I am sorry to have added some confusion in my OP which indeed contained a typo. Free42 and Legendary67 returned the same result as Excel, i.e. 146:54:24 (not 44).

So Retro15C is the culprit.

[Albert Chan]

(09-07-2018 02:40 PM)Thomas Okken Wrote:  The approach I chose for ->HR for the binary version is robust but not very pretty: multiply x by 1e12,
round to nearest integer, and then do everything using integer math.

Plus, there is a special case in for |x| < 0.01, where the conversion reduces to a simple division by 0.36.

Hi, Thomas Okken

Just a crazy example, billion degrees 4 minutes.

X = 1000000000.03999996185 ... (binary float)
X * 1e12 = 1000000000039999897600 (binary, exact integer)

With this value, result will be over-counted ~ 100 - 60 = 40 seconds (just like Retro12C)
(note: this is only based on your comment, I have not try it on the actual emulator)

Instead of nearest integer, it might be better with rounded 12 digits value (no scaling needed).
12 digits can safely pushed to 15, see https://www.exploringbinary.com/number-o...nversions/

For accurate dtoa, David Gay's dtoa.c is well tested: http://www.netlib.org/fp/dtoa.c
Or, you can try my code: https://github.com/achan001/dtoa-fast

Edit: some good C compiler can do correct strtod and dtoa functions. With those, above code is not needed

[Thomas Okken]

It doesn't multiply the entire number by 10^12, just the fractional part. But, yes, with a case like that, numerical cancellation will cause the binary round-off to mess up the result. Avoiding that is possible but ugly (I think the 33S does this).

This is why Free42 Decimal exists.

UPDATE: Hadn't read your post carefully enough. Better binary-to-decimal does sound like a good idea. I'll look into it sometime after my vacation...

[Albert Chan]

(09-08-2018 02:35 PM)Thomas Okken Wrote:  It doesn't multiply the entire number by 10^12, just the fractional part. But, yes, with a case like that, numerical cancellation will cause the binary round-off to mess up the result. Avoiding that is possible but ugly (I think the 33S does this).

I thought of that (scale only fractional part) before posting.
But, whether or not the integer part got removed before scaling does not matter:

X - INT(X) is exact = 0.03999996185 ... (DMS->DEG over-counted 40 seconds still here)

In fact, my suggestion of rounded 15 digits required INT(X) not to be removed.

Rounded X (to 15 digits) = "1000000000.04" (character string)

[Albert Chan]

(09-07-2018 03:43 PM)Albert Chan Wrote:  Example: if angle = 2° 1' , X = 2.01 (IEEE double of 2.0099999...)

X = 2.0099999 ... --> Y = 20100.00 (scaled and rounded to hundredth of second)

DMS->DEG: (Y - 40*int(Y/100) - 2400*int(Y/10000)) / 3600 = (20100 - 8040 - 4800) / 3600 ~ 2.0166666...

Uh, I was getting stupid ...
My post #23 already have an elegant method.

Scale by 10000 (to move the troublesome decimal point), then round to 15 digits.
(approximated rounding good enough, just to catch the "40 seconds" bug)

That is all that is needed.
No reason to get precisely rounded 15 digits decimal.

[Albert Chan]

Just think of a great idea: fast, simple, accurate DMS->DEG

The whole point of scaling / rounding is to avoid the "40 seconds" bug.
All we need is to give up 1 ULP binary float accuracy, and we are safe.
(correct strtod conversion guaranteed max error <= 0.5 ULP)

Using this thread original example, X = 183.3799999 ...

X = nextafter(X, X+X) = 183.38000000000002387 ...

Now, safe to apply Pekis DMS->DEG formula, and we are done

(09-07-2018 07:07 AM)Pekis Wrote:  And with my funny one-liners:
DEG->DMS: (90*X+INT(60*X)+100*INT(X))/250
DMS->DEG: (250*X-INT(100*X)-60*INT(X))/90

[Joe Horn]

Unfortunately Pekis' HMS formula suffers from the "60 second bug" which several HP's also suffered from. If you input 9.99999999999 (12 digits total) into any 12-digit HP, and run Pekis' HMS formula, you get exactly 9:59:60 as the result. Should be 9:59:59.9999999 of course.

[Carsen]

(09-09-2018 05:10 AM)Joe Horn Wrote:  Unfortunately Pekis' HMS formula suffers from the "60 second bug" which several HP's also suffered from. If you input 9.99999999999 (12 digits total) into any 12-digit HP, and run Pekis' HMS formula, you get exactly 9:59:60 as the result. Should be 9:59:59.9999999 of course.

The way I see it is that it is more of an "error" of the calculator limited precision. If the calculator had a higher precision, it would certainly return the correct answer with the right amount of 9's trailing at the end of it.

[Albert Chan]

Hi, Joe Horn

I was in the middle of removing Pekins formula ...
DMS->DEG may have problems too, with INT(100*X)
Doing conversion with integer degree and minutes is safer (float seconds OK):

Y = 100*X
Y = nextafter(Y, Y + Y)
DMS->DEG = (2.5*Y - int(Y) - 60*int(Y/100)) / 90

Max error = decimal to binary error + scaling error = 0.5 + 0.5 = 1.0 ULP
--> 1 ULP is enough to get to safety, away from "40 seconds bug"

Edit: conclusion is correct, but ULP math was wrong, see post 55

[Werner]

(09-09-2018 05:10 AM)Joe Horn Wrote:  Unfortunately Pekis' HMS formula suffers from the "60 second bug" which several HP's also suffered from. If you input 9.99999999999 (12 digits total) into any 12-digit HP, and run Pekis' HMS formula, you get exactly 9:59:60 as the result. Should be 9:59:59.9999999 of course.

Hi, Joe!

9.99999999999 ->HMS is 9:59:59.999999964, exactly, and so the rounded 12-digit answer should be 10:00:00, which is what the 48G gets. Unless you meant that Pekis' formula should have returned either 10 or 9:59:59.9999999

Cheers, Werner

[Werner]

A better HMS formula, courtesy of John H. Meyers in comp.sys.hp48, long ago:

H.MS=H.-FP(H.)*0.4-FP((H.-FP(H.)*0.4)*100)*0.004

or

T := H.- FP(H.)*0.4 ;
H.MS := T - FP(T*100)*0.004 ;

translated into RPN that would be

ENTER
FRC
.4
x
-
1E4
%
FRC
4E-3
x
-

At least, it returns 9.59599999999 for 9.99999999999 ..

There's an equivalent HR formula:

H.=H.MS+FP(H.MS*100)/150+FP(H.MS+FP(H.MS*100)/150)/1.5

or

T := H.MS + FP(H.MS*100)/150;
H. := T + FP(T)/1.5;

Cheers, Werner

[teenix]

(09-09-2018 05:10 AM)Joe Horn Wrote:  Unfortunately Pekis' HMS formula suffers from the "60 second bug" which several HP's also suffered from. If you input 9.99999999999 (12 digits total) into any 12-digit HP, and run Pekis' HMS formula, you get exactly 9:59:60 as the result. Should be 9:59:59.9999999 of course.

The HP67 gives

9.5959999999 in DSP 9 mode
9.5960 in DSP 4 mode
9.60 in DSP2 mode
10 in DSP 0 mode

It is just following its normal display rounding code. In this case, it doesn't know that certain digits are hours, minutes, or seconds.

cheers

Tony

[Albert Chan]

(09-09-2018 06:13 AM)Albert Chan Wrote:  Y = 100*X
Y = nextafter(Y, Y + Y)
DMS->DEG = (2.5*Y - int(Y) - 60*int(Y/100)) / 90

Max error = decimal to binary error + scaling error = 0.5 + 0.5 = 1.0 ULP
--> 1 ULP is enough to get to safety, away from "40 seconds bug"

Above can be simplified and maybe more accurate by using non-integer factors

Y = 100.00000000000001 * X
DMS->DEG = (2.5*Y - int(Y) - 60*int(Y/100)) * 0.011111111111111108

first factor (with 15 zeroes) = 100.000 ... = 1 ULP above 100.0, to fix "40 seconds" bug
second factor (with 15 ones) = 0.01111... = 2 ULP below 1/90.0, to compensate over-estimation

For angle=183°38' , above produce error of +0.73 ULP (against exact 183 + 38/60)
For 1000000000°4', above produce error of +0.47 ULP ... not too bad.

[Thomas Okken]

The problem with using nextafter() is that the closest decimal equivalent isn't always going to be higher. In other words, you're always rounding up, but sometimes 0.3999 really is 0.3999, not a mangled 0.4000.

How about this: instead of multiplying by 10^12, I'll multiply by whichever power of 10 makes the integral part of the number have 12 digits (or 15, or whatever choice is safe), then round to the nearest integer, then remove what was originally the integer part, and finally normalize the scale of the remainder. That should preserve the nice properties of the current Free42 Binary implementation, while improving its handling of binary round-off as much as possible (without doing the equivalent of a full-blown dtoa(), anyway).

[Albert Chan]

(09-09-2018 01:54 PM)Thomas Okken Wrote:  The problem with using nextafter() is that the closest decimal equivalent isn't always going to be higher. In other words, you're always rounding up, but sometimes 0.3999 really is 0.3999, not a mangled 0.4000.

Please see my previous post #38. I decided not to do nextafter(), saving a step.
Always rounding-up also have a compensating always rounding-down.

Anyway, we are talking 1 ULP adjustment, not thousands (rounding to 12 digits).
The effect is so tiny (just enough to avoid the bug), you barely see it:

X = 0.3998999 ... (39 minutes, 99 seconds)
Y = X * 100.00000000000001 = 39.9900000000000019 ... ~ 39.99 (still 99 seconds)
DMS->DEG = 0.6774999999999999 (1.08 ULP below exact 39/60 + 99/3600)

X = 0.4699999 ... (47 minutes, but possible "40 seconds" bug)
Y = 47.000000000000007 ...
DMS->DEG = 0.7833333333333333 (0.07 ULP below exact 47/60)

Even if user input is 0.46999 99999 99999 8, above code still intepret it as 46 minutes 100 seconds.
So, no worries about mis-intepreting user input, 0.3999 will not mangled to 0.4000

Last 8 changed to 9 will confuse it (as 47 minutes). But, you are now reaching beyond IEEE double ...