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

[Pekis]

Hello,

Please note the proposed patch in my DEG->DMS formula to fix the "60 seconds" bug by naturally rounding to internal precision:

INT(X*60) is transformed into INT(X*120/2)
INT(X) is transformed into INT(X*2/2)

(patched): DEG->DMS: (90*X+INT(X*120/2)+100*INT(X*2/2))/250
(unchanged): DMS->DEG: (250*X-INT(100*X)-60*INT(X))/90

What do you think of it ?

[Albert Chan]

(09-10-2018 07:31 AM)Pekis Wrote:  (patched): DEG->DMS: (90*X+INT(X*120/2)+100*INT(X*2/2))/250
(unchanged): DMS->DEG: (250*X-INT(100*X)-60*INT(X))/90

What do you think of it ?

I propose patch DMS->DEG, to fix the "40 seconds" bug

DMS->DEG: (250*X-INT(100.00000000000001*X)-60*INT(X))/90

formula is basically unchanged, until X in "40 seconds bug" region.

In other words, X is at most 2 ULP away from "buggy" 100 seconds.
The patch remove the 100 seconds, round-up an extra minute instead.
(Not precisely 1 ULP fix with nextafter(), but good enough)

Your patch DEG->DMS have no effect for binary float, *2 and /2 is exact.

For calculator that used BCD, *2 /2 trick only work if X is "filled" with digits.
But, not always. On HP-12C, X=9.999999995, DEG->DMS = 9.5960 (exact)

Since 9°59'60" = 10°, "60 seconds bug" is not serious (easily fixed in the head)
At worse, this is a display bug, wrongly assumed number is decimal (see post 37).

[Pekis]

(09-10-2018 12:46 PM)Albert Chan Wrote:  Your patch DEG->DMS have no effect for binary float, *2 and /2 is exact.

For calculator that used BCD, *2 /2 trick only work if X is "filled" with digits.
But, not always. On HP-12C, X=9.999999995, DEG->DMS = 9.5960 (exact)

OK, but at least for simple and usual calculations (in BCD), it removes the annoying "60 seconds bug". For example, DEG->DMS(DMS->DEG(0.08)) goes back to 0.08, and not 0.0760.

[Werner]

FWIW, the formula

HR=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;
HR := T + FP(T)/1.5;

that I posted earlier has no problems with the billion degrees 4 minutes problem, and doesn't have to resort to introducing errors of a few ulps.

Cheers, Werner

[Joe Horn]

(09-10-2018 12:46 PM)Albert Chan Wrote:  Since 9°59'60" = 10°, "60 seconds bug" is not serious (easily fixed in the head)
At worse, this is a display bug, wrongly assumed number is decimal (see post 37).

(1) "Not serious"? Would anybody tolerate a digital clock that said the time is 1:29:60? No. It would be considered a buggy clock. Nor should we tolerate an HMS algorithm that returns 60 seconds or 60 minutes.

(2) "Only a display bug" is not true. Pekis' formulas (both his original and his "improved" one) return EXACTLY 1.29600000000 on the HP 50g for an input of 1.49999999999. That's not an artifact of display-mode rounding. It's a bug in the algorithm.

[Albert Chan]

(09-10-2018 01:31 PM)Pekis Wrote:  for simple and usual calculations (in BCD), it removes the annoying "60 seconds bug".
For example, DEG->DMS(DMS->DEG(0.08)) goes back to 0.08, and not 0.0760.

I mis-understood your patch, thinking it cover all cases ...
If simple cases round-trip are what you are after, *2/2 can be removed, to this:

(patched): DEG->DMS: (90*X+INT(X*120/2)+100*INT(X))/250

(09-10-2018 02:10 PM)Werner Wrote:  T := H.MS + FP(H.MS*100)/150;
HR := T + FP(T)/1.5;

that I posted earlier has no problems with the billion degrees 4 minutes problem,
and doesn't have to resort to introducing errors of a few ulps.

Without error analysis (ULPs), you cannot get this right.
There are too many cases above formula can fail: 29', 57', 58', 1°13', 1°14', 1°15', 1°16', ...

Example: 57 minutes should convert to 57/60 = 0.95 degrees

T = 0.56999 ... + 0.00666 ... = 0.57666 ... (I assumed FP meant fractional part)
HR = 0.57666 ... + 0.38444 ... = 0.96111 ... (again, off by 40 seconds)

[Albert Chan]

(09-10-2018 03:32 PM)Joe Horn Wrote:  

(09-10-2018 12:46 PM)Albert Chan Wrote:  Since 9°59'60" = 10°, "60 seconds bug" is not serious (easily fixed in the head)
At worse, this is a display bug, wrongly assumed number is decimal (see post 37).

(1) "Not serious"? Would anybody tolerate a digital clock that said the time is 1:29:60? No. It would be considered a buggy clock. Nor should we tolerate an HMS algorithm that returns 60 seconds or 60 minutes.

(2) "Only a display bug" is not true. Pekis' formulas (both his original and his "improved" one) return EXACTLY 1.29600000000 on the HP 50g for an input of 1.49999999999. That's not an artifact of display-mode rounding. It's a bug in the algorithm.

Calculator really need to know the fractions is minutes and seconds for correct round-up.
In this regard, my Casio get it right, dedicated DMS key for degree/minute/second.
So it knows 60 seconds to a minute, 60 minutes to a degree.
Best of all, user know it too, since it is displayed differently, as in 1°2°3°

Unfortunately, HP pick DD.MMSS decimal representation (a normal number !)
Even the calculator itself got confused. You can press DEG->DMS repeatedly, getting garbage.
It does not know ... (see post 37)

As long as calculator treat the "number" as decimal, it cannot be fix by patching formula.
Raw result of exactly 1.296, or will rounded as 1.296 doesn't matter.
Either way, the result need some cleaning up.

When I say, "not serious", "only a display bug", I meant cleaned-up result comes out the same.

Yes, I agreed this is not good. Did HP Prime still use this DD.MMSS form ?

[Joe Horn]

(09-10-2018 06:18 PM)Albert Chan Wrote:  Did HP Prime still use this DD.MMSS form ?

No, it uses a special display format, similar to Casio. For example, →HMS(pi) yields 3º08'29.73355".

Note: I just looked up the HP48 →HMS function in HP's source code, and it mentions that REV:A of the HP 48 had the 60 second bug! It even gives 1.49999999999 --> 1.2960 as the example! Golly. But then it says that the bugfix "sacrifices a small amount of accuracy in final result due to intermediate rounding." Apparently they multiplied the 15-form number by 0.6 and then rounded the last 3 digits (not sure what that's all about). They said, "The rounding is not essential, but avoids undesirable rounding effects associated with packing the final h.ms result into a 12-digit destination." So they struggled with this problem too, and opted for slightly less accuracy to avoid 60 seconds from ever appearing in the result. That's similar to the Prime, which only shows a maximum of 5 decimal places for the seconds in HMS notation.

[Joe Horn]

(09-10-2018 06:18 PM)Albert Chan Wrote:  As long as calculator treat the "number" as decimal, it cannot be fix by patching formula.

Ho ho! The gauntlet has been thrown! <scampers off to Happy Programming Place to ponder this challenge>

[Dieter]

(09-10-2018 06:53 PM)Joe Horn Wrote:  

(09-10-2018 06:18 PM)Albert Chan Wrote:  Did HP Prime still use this DD.MMSS form ?

No, it uses a special display format, similar to Casio. For example, →HMS(pi) yields 3º08'29.73355".

For the record: the WP31s offers both ways: you can temporarily display an angle in d°mm'ss,ss" format or you can convert it as usual to d.mmssss.

(09-10-2018 06:53 PM)Joe Horn Wrote:  Note: I just looked up the HP48 →HMS function in HP's source code, and it mentions that REV:A of the HP 48 had the 60 second bug! It even gives 1.49999999999 --> 1.2960 as the example! Golly. But then it says that the bugfix "sacrifices a small amount of accuracy in final result due to intermediate rounding."

Simply do a→H and →H.MS sequence and everything is fine again.

On the other hand the HP35s gets it right: it returns exactly 1,3 (1°30'). I assume this is also true for later HP48 ROM revisions.
The 16-digit WP31s displays 1,2960 but this is only due to display rounding. The full 16-digit result is 1,295999999996400.

BTW, both the '41 and the '67 return 1,295999999 for an input of 1,499999999. The true result is 1,29599999964 so this should have been rounded to 1,2960 or 1,3000.

Dieter

[Werner]

(09-10-2018 12:46 PM)Albert Chan Wrote:  I propose patch DMS->DEG, to fix the "40 seconds" bug

DMS->DEG: (250*X-INT(100.00000000000001*X)-60*INT(X))/90

Why not

DMS->DEG: (250*X-INT(100*X+0.25)-60*INT(X))/90

- either X is 'well behaved' and it makes no difference
- or SS=0 and H.MM is not exactly representable (the 40-seconds), and then it rounds up. Independent of the precision used.

Cheers, Werner

[Werner]

(09-10-2018 04:28 PM)Albert Chan Wrote:  Without error analysis (ULPs), you cannot get this right.
There are too many cases above formula can fail: 29', 57', 58', 1°13', 1°14', 1°15', 1°16', ...

Example: 57 minutes should convert to 57/60 = 0.95 degrees

T = 0.56999 ... + 0.00666 ... = 0.57666 ... (I assumed FP meant fractional part)
HR = 0.57666 ... + 0.38444 ... = 0.96111 ... (again, off by 40 seconds)

Oops, of course. Uncanny how you can come up with counterexamples so easily.

But then, instead of FP(100*H.MS) - which is the culprit- , all you need to do is
to take H.MS*100 - IP(H.MS*100+0.25) instead.

So,

T := HMS*100;
T := T - IP(T+0.25);
T := HMS + T/150;
HR := T + FP(T)/1.5;

or, in RPN:

Code:

>LBL "DMS>"
 1E4
 %
 4
 1/X
 RCL+ ST Y
 IP
 -
 150
 /
 +
 ENTER
 FP
 1.5
 /
 +
 END

Seems to work in Free42 Binary.

Cheers, Werner

[Albert Chan]

(09-12-2018 07:06 AM)Werner Wrote:  Why not

DMS->DEG: (250*X-INT(100*X+0.25)-60*INT(X))/90

(09-12-2018 07:14 AM)Werner Wrote:  T := HMS*100;
T := T - IP(T+0.25);
T := HMS + T/150;
HR := T + FP(T)/1.5;

Correctly implemented strtod guaranteed maximum error 0.5 ULP
Scaling by constant also guaranteed maximum error of 0.5 ULP
So, scaled value can at most be 1 ULP below User Decimal Input (where 40 second bug might hit)
Ideally, we want 1 ULP correction: nextafter(100*X, 200*X)

X*100.00000000000001 might correct 1 to 2 ULP (due to rounding).
But, (possible) 1 ULP over-correction is good enough.
The patched version is also fast, same as un-patched version.

Too big over-correction, however, may mis-interpret what the user entered.
Both Pekis and Meyers "+25 seconds" patches have similar counter-examples:

For negative angles, patches make it worse, say, X = -0.10, under-counted 40 seconds.
For positive angles, 75+ seconds, it "fixed" a phantom bug, under-counted 40 seconds.

Example: 85 seconds, both patches returned 0.0125 degree (= 45 seconds).

The patches are not wrong, but limited (non-negative angles and no "bad" seconds).

[Albert Chan]

(09-12-2018 02:08 PM)Albert Chan Wrote:  Correctly implemented strtod guaranteed maximum error 0.5 ULP
Scaling by constant also guaranteed maximum error of 0.5 ULP
So, scaled value can at most be 1 ULP below User Decimal Input (where 40 second bug might hit)

I made an error adding up the ULPs.
What should be added is relative errors, not ULPs

Let binary mantissa of X = m1, of scaled Y = 100*X = m2
So, 2^52 <= m1 < 2^53, and same for m2 (units in ULP)

Max relative (under-estimated) error = 0.5/m1 + 0.5/m2

Max absolute error (ULP) = Max relative error * m2
= 0.5 + 0.5 m2/m1
= 0.5 + 0.5 * max(m2/m1)
= 0.5 + 0.5 * (2^53-1) / 2^52
< 1.5 ULP

Technically +1 ULP still work, but it is cutting it very close.
(If Y < 0.5 ULP below an integer, it will round-up to integer)

100.00000000000001*X patch is safer:
Let m3 = binary mantissa of 100.0 = 0.78125 * 2^53

Max relative (under-estimated) error = 0.5/m1 + 0.5/m2 - 1/m3

Max absolute error (ULP) = Max relative error * m2
= 0.5 + 0.5 m2/m1 - m2/m3
= 0.5 + m2 * (1/(2 m1) - 1/m3)

Since (2 m1) > m3, the factor is negative, we need to minimize both terms

Max absolute (under-estimated) error
= 0.5 + 2^52 * (1/2^53 - 1/(0.78125*2^53))
= 0.36 ULP < 0.5 ULP

In other words, DD.MM * 100.00000000000001 >= DDMM (in binary float)

[Albert Chan]

Above error bound can be tighten a bit, since we know scaling factor ~ 100.0

m1 and m2 are related by the factor k = m3/2^53 = 0.78125

If m1/2^53 >= 0.64, m2/m1 = k, otherwise m2/m1 = 2*k

With factor 100.0:

Max absolute error = 0.5 + 0.5 m2/m1 = 0.5 + k ~ 1.28 ULP
+1 ULP fix reduce this to k - 0.5 ~ 0.28 ULP

With factor 100.00000000000001:

Max absolute error = 0.5 + 0.5 m2/m1 - min(m2/m3)

If m1/2^53 >= 0.64, error = 0.5 + 0.5 k - 0.5/k ~ 0.25 ULP
Otherwise, error = 0.5 + 0.5 (2 k) - k/k = k - 0.5 ~ 0.28 ULP

Both approach have the same lower bound !

[Pekis]

Hello,

Could this formula be tested on various machines ?

I'm trying to help rounding up to next integer when it's very close.
So, I transformed
INT(X) into INT(X+0.75*(1-FP(X)))
and INT(60*X) into INT(60*X+0.75*(1-FP(60*X)))

the DEG->DMS(X) would then be
DEG->DMS(X)=(90*X+INT(60*X+0.75*(1-FP(60*X)))+100*INT(X+0.75*(1-FP(X))))/250

Or if using INT only:
DEG->DMS(X)=(90*X+INT(60*X+0.75*(1-60*X+INT(60*X)))+100*INT(X+0.75*(1-X+INT(X))))/250

I think it's worth trying and I'm ready to withstand the expected bombing

Thanks

[Albert Chan]

(09-13-2018 08:47 AM)Pekis Wrote:  I'm trying to help rounding up to next integer when it's very close.
So, I transformed
INT(X) into INT(X+0.75*(1-FP(X)))
and INT(60*X) into INT(60*X+0.75*(1-FP(60*X)))

My advise is give up fixing the 60 seconds bug (sorry for the put down ...)
The fix may be worse for the "numbers". (and you don't know when that hit)

Before the fix, we may see 60 seconds, sometimes 60 minutes, but it is accurate.
A little cleaning-up, and we are good ...

After the fix, potentially you get into fixing phantom bug, similar to +25 seconds patch:
Just an example, on HP-12C, try X = 1000.999999:

Pekis DEG->DMS formula: (90*X + INT(60*X) + 100*INT(X)) / 250

Old DMS = (90*X + 60059 + 100*1000)/250 = 1000.5960 (= 1001 degree)

Using above "fuzzy" INT for INT, we get:

New DMS = (90*X + 60059 + 100*1001)/250 = 1000.9960 (= 1001 degree 40 minutes)

Similar search should discover "40-seconds bug" too ...

[Pekis]

It was fun ... Thanks for your answers ... Waiting for next time

[Werner]

(09-12-2018 02:08 PM)Albert Chan Wrote:  The patches are not wrong, but limited (non-negative angles and no "bad" seconds).

You're absolutely right - I specifically assumed correct input, and forgot about negative angles.

Cheers, Werner

[AndiGer]

(09-09-2018 08:59 AM)teenix Wrote:  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

Interesting observation ...
Both of my HP67 (1610S, 1705S) give 9.595999999 in DSP 9 mode. But 0.596000000 for both 0.999999999 (9 times 9 input) and .9999999999 (10 times 9 input).
On the other hand n.595999999 for any integer n between 1 and 9.

Legendary 67 gives n.596000000 for any value of the form n.595999999 for n between 0 and 9. But 1.000000000 for .9999999999 (10 times 9).
Andi