Numerical accuracy on the HP Prime (in Python, Home, and CAS)

[ftneek]

The main purpose of this thread is to try and gain insight about how a cosine value is computed on the Prime, Xcas, etc., discuss limitations of floating point arithmetic, constraints of the HP Prime, etc., and see if it would be possible to increase the accuracy of some operations (like cos, MOD, etc) over floating point arithmetic.

I recently noticed this thread, where it was reveled that the G2 returns an incorrect result of -Inf when evaluating cos(1e50) from the CAS View. This and this thread also had some interesting details.

It made me curious, what should the correct (or expected) result be? I decided to check some different sources, put the results into a spreadsheet, and tried to sort it in a way that made sense (to me). I’ve bolded the results from the HP Prime, results agreeing with other platforms until the last digit are highlighted in green, and results not agreeing with other platforms are highlighted in red. I’ve also attached a python script which simply prints the result of each calculation in case you’d like to try it on another platform and compare the results.

results:
  results_cos(1e50).xlsx (Size: 12.3 KB / Downloads: 20)

python script:
  cos(1e50).txt (Size: 671 bytes / Downloads: 8)

As you can see, the HP Prime’s Python App environment had 8 results that were significantly different from Python or other MicroPython ports (half of them a reported domain error), and 4 results that differed only in the last digit. I won’t complain about the last digit being different, but I’m curious to know a reason for the difference. It seems the Prime might have rounded up while the other Python environments did not.

The Home environment matches the result on Wolfram alpha when in Degrees and Gradians mode, but not Radians. Presumably the G2’s CAS should match the results from the 14425 VC, though it does not seem to be as accurate as python or Wolfram alpha, likely due to the nature of CAS being an exact environment.

Still, is it possible to do better than the current results in Home and CAS?

My first thought is that for an angle less than 0 or greater than 2*pi you can add or subtract multiples of 2*pi until it is within the interval [0,2*pi], or in other words use the congruent angle mod 2*pi. This seems to happen automatically in most other Python implementations.

Home:
13 MOD PI -> 0.43362938564 (ok)
13 MOD 2*PI -> 3.14159265359
The last result should have matched the one above it.
13 -2*(2*pi) -> 0.43362938564 (ok)
13 MOD (2*PI) -> 0.43362938564 (ok)
(1ᴇ50) MOD 2*PI -> 0
1ᴇ50 MOD (2*PI) -> 6.149975477

There is an issue when the dividend must be evaluated and is not enclosed in parentheses. A result of 3.14159265359 does not make sense if the dividend was 2 or PI, let alone 2*PI. It seems MOD is severely limited by floating point roundoff. The last MOD statement shown does not agree with Python (~3.92) or Wolfram alpha (~4.05). The Prime could probably benefit from an improved floating point MOD routine, even for fairly simple inputs. (parentheses bug).

CAS (Exact mode):
13 MOD PI -> -4*π+13 (nice)
13 MOD 2*PI -> π (should have matched above)
13 MOD (2*PI) -> -4*π+13

CAS (Exact or approx does not matter):
1e50 MOD (2*PI) -> 6.64613997892e35

The parentheses issue is present in CAS as well, but at least CAS can return an Exact answer. Unfortunately 1e50 was not reduced to a value less than 2*pi (not even just the significant digits). Maybe another improvement to consider for Home would be to return Exact MOD results in some cases if Intelligent Math is enabled.

So, how would you compute cos(1e50) on a G2 as accurately as possible? Is it possible to beat the Python results? Feel free to share your comments, ideas, programs, results from other calculators, other problems where the numerical accuracy could potentially be improved through one method or another, or explain why it wouldn’t be possible to compute with more precision.

[Nigel (UK)]

One comment: on the Prime, 1e50 is taken by the CAS to be an approximate number. It would be better to ask for cos(10^50), and then use approx() to get a decimal result.

Xcas gives the correct answer for approx(cos(10^50)) (if DIGITS is set to a large enough value) but the Prime gives the wrong answer still. I don’t think you can avoid this on the Prime, where CAS precision is 15 or 16 digits at most. You would surely need more than 50 digits of PI to get any correct digits for 10^50 MOD (2*PI).

An unrelated point: I think that 13 MOD 2*PI = 3.14159… is fine; the calculator is working from left to right, so it’s doing (13 MOD 2) * PI.

Nigel (UK)

[ftneek]

Thank you, that result makes sense now. Maybe it is not a parentheses 'bug' after all. Likely, it is just another instance where the lack of parentheses makes the problem ambiguous or up for interpretation depending on operator precedence.

If I wanted the result of the MOD operation to multiplied by pi I would have written it as (13 MOD 2)*pi, however at the same time you can probably say that if I wanted it the way I was expecting, I should have written it as 13 MOD (2*pi). Is one interpretation generally followed in different programming languages?

And you are right I did not think to increase the precision in Xcas, now I see a much more accurate result. Perhaps I'll try to add the precision used later, though in all my tests I believe I just used the default settings.

My other concern is why the G2's Python results seem to be (very) different sometimes and why they're sometimes different by 1 digit.

edit: if I were to make a suggestion regarding the MOD operator, it would be to add a Textbook Display template, so that it is displayed as ([] MOD []) and you simply need to fill in the boxes, that should make the order of operation more clear.

[parisse]

1e50 mod 2*pi can not be computed with any significant digit using "normal" precision. The reason is that the relative error rounding any floating point number will be say 1e-17 (it's 2^(-53) for double precision inside MicroPython, 2^(-48) for Xcas/CAS normal precision), that's an absolute error of 1e33 for 1e50, and that's much greater than 2*pi.
This is something that everyone doing scientific computation should know, if it's not teached in a course during highschool or first years of University, there is a real problem.

[ftneek]

Thank you parisse. Rounding error is bound to happen when working with limited precision, the results of such a computation are expected to have different error margin in different precision environments.

For clarification, I did not mean to imply the result of limited precision float arithmetic was incorrect. Part of my curiosity was to see if someone would be able to simulate 'extended precision' for a MOD operation through a user program, perhaps using strings even if inefficient, simply for novelty or amusement.

I was confused about the operator precedence at first, and added a suggestion regarding placing parenthesis around MOD when in textbook display mode.

Another reason for the thread was to note some differences between Python and the MicroPython implementation on the Prime, so that users might be able to adjust their scripts accordingly if necessary.