Opened 9 years ago

Closed 7 years ago

#10164 closed defect (fixed)

Few digits of precision in N().

Reported by: gerbicz Owned by: jason
Priority: major Milestone: sage-5.1
Component: misc Keywords: N, digits, numerical approximation
Cc: kcrisman Merged in: sage-5.1.beta3
Authors: Robert Gerbicz, Douglas McNeil Reviewers: Karl-Dieter Crisman, Benjamin Jones
Report Upstream: N/A Work issues:
Branch: Commit:
Dependencies: Stopgaps:

Description (last modified by benjaminfjones)

For N(1+10**-400000,digits=400001) the displayed number of digits is only 400000, it is printing 1. followed by 399999 zeroes. The reason is that in functional.py: prec = int((digits+1) * 3.32192) + 1. However log(10)/log(2)~3.3219280948874>3.32192, so if digits is large the used precision will be smaller by some digits than the requested number of digits.

The suggestion is to use 3.32193 instead of 3.32192, see the trivial patch.


Apply: trac_10164_folder_sans_whitespace.patch to Sage library

Attachments (6)

10164.patch (644 bytes) - added by gerbicz 9 years ago.
trac_10164_round_up_log10_2.patch (8.0 KB) - added by dsm 7 years ago.
roud up log(10,2)
trac_10164_round_up_log10_2.2.patch (8.0 KB) - added by dsm 7 years ago.
round up log(10,2)
trac_10164_round_up_log10_2_v3.patch (8.3 KB) - added by dsm 7 years ago.
v3: use a more precise approx
trac_10164_reviewer.patch (10.6 KB) - added by benjaminfjones 7 years ago.
name the magic 3.23... number, remove trailing whitespace
trac_10164_folder_sans_whitespace.patch (12.9 KB) - added by benjaminfjones 7 years ago.
last two patches folded togethers with whitespace removed

Download all attachments as: .zip

Change History (33)

Changed 9 years ago by gerbicz

comment:1 Changed 7 years ago by kcrisman

  • Cc kcrisman added
  • Keywords beginner added

See also this sage-devel thread. Doug reports the following evil lines:

misc/functional.py:            prec = int((digits+1) * 3.32192) + 1 
rings/complex_interval.pyx:    bits = 
max(int(3.32192*len(s_real)),int(3.32192*len(s_imag))) 
rings/complex_number.pyx:    bits = 
max(int(3.32192*len(s_real)),int(3.32192*len(s_imag))) 
rings/real_mpfi.pyx:        bits = int(3.32192*len(s)) 
rings/real_mpfr.pyx:            bits = int(3.32192*sigfigs)+1 
symbolic/expression.pyx:                prec = int((digits+1) * 3.32192) + 1 

So a full patch would fix all of these, though we definitely want to keep the original author on. Also need doctests for these large cases, of course.

I can't believe I never saw this ticket.

comment:2 Changed 7 years ago by dsm

I'm on this one.

comment:3 Changed 7 years ago by dsm

Okay, I have a patch. Testing now. Turns out it's a little trickier to doctest than I thought, because some of the bugs can't be triggered yet because of other limitations. :-/

Changed 7 years ago by dsm

roud up log(10,2)

Changed 7 years ago by dsm

round up log(10,2)

comment:4 follow-up: Changed 7 years ago by dsm

  • Authors changed from Robert Gerbicz to Robert Gerbicz, Douglas McNeil
  • Status changed from new to needs_review

Okay, so I've changed all six instances of "3.32192" (which seems to be the most common approximation to log(10,2) in use). There are three other references matching 3.32* (in rings/arith.py, interfaces/gp.py, and libs/pari/gen.pyx) I left alone. Those all have at least 10 digits, so the problem would be pushed to higher levels even if it occurs, and there was in every case something I wasn't confident enough of to change. I don't know much about how pari's precision works, for example.

I think it's a step in the right direction, anyhow.

comment:5 in reply to: ↑ 4 Changed 7 years ago by kcrisman

  • Milestone set to sage-4.8
  • Reviewers set to Karl-Dieter Crisman
  • Status changed from needs_review to needs_work

I think it's a step in the right direction, anyhow.

Very clever testing.

Very minor 'needs work' issues:

  • I think you could also add as a doctest the example at #12163 (a dup, more or less). We do like doctesting problems people officially put on tickets. That sort of example has the added benefit of being a little easier to read :) Of course, if more than one person disagrees, it's not a big deal.
  • Needs to be done - your doc lines are really long. Compare
    s_real -- a string that defines a real number (or something whose
    
    to
    Make sure we've rounded up log(10,2) enough to guarantee sufficient precision (trac #10164).
    
    Try to shorten those puppies up to the 80-character limit (if I recall correctly) for fine terminal viewing.

Because this has the potential to do weird things to doctests, I'm not going to give positive review in any case until I run long tests. But looks good other than those tiny issues.

Also, a question: I assume it's okay (and already the case, in fact) that we have MORE digits than "officially requested" by the user.

comment:6 Changed 7 years ago by schilly

I like this so far, but I have a maybe dumb question. Why not use more digits? in "raw" python, 3.3219280948873626 is the answer on my 64bit machine. changing the last digit from 6 to 7 should be just perfect.

comment:7 Changed 7 years ago by kcrisman

Long doctests are a-ok, so just the patch itself would be treated differently based on comments.


Harald's comment seems quite reasonable. One could go either way, though, unless we are really trying to get exactly the 'right' number of digits, which apparently hasn't been the case (we treat as lower bound).

comment:8 Changed 7 years ago by eviatarbach

Why not get exactly the correct number of digits? Certainly the added precision will fix the problem for much larger numbers.

The question is how many digits we can feasibly expect people to want to calculate. 3.32193 works for 1000000 digits, 3.3219281 works for 10000000, but I don't see why we can't use 3.3219280948873627 to future-proof it.

To calculate the precision error with different approximations of log(10, 2):

(floor((n + 1) * a) + 1) - (floor((n + 1) * log(10,2).n(100000)) + 1),

where n is the number of digits you want, and a is the approximation. log(10,2).n(100000) is just used for comparison.

comment:9 Changed 7 years ago by dsm

I simply followed the OP's suggestion, which works just fine: adding more digits would only change things by ~1e-7 or so.

"Why not get exactly the correct number of digits?"

Well, since I don't know if the last digits are guaranteed to be correct anyway, I'm not sure how important it is to be precise when we're not accurate, esp. given that we don't seem to have been in the past. Of the two, not using as much precision as requested is a far more serious problem than using a (very tiny) bit more. But it's probably easier for us simply to change it than to waste time arguing about the use case. :-)

As for: "3.32193 works for 1000000 digits, 3.3219281 works for 10000000, but I don't see why we can't use 3.3219280948873627 to future-proof it."

I think there's some confusion here. 3.32193 works for _any_ number of digits, because it's an over-estimate. It's already completely future-proofed.

comment:10 Changed 7 years ago by eviatarbach

I see your point. I guess the way to do it would be to use the overestimate and slice to the correct number of digits, but this may be difficult since there doesn't seem to be an index method for RealNumber?.

comment:11 Changed 7 years ago by ncarter

  • Keywords beginner removed

This has gone beyond beginner capacities, due to the extended discussion above.

comment:12 Changed 7 years ago by eviatarbach

  • Status changed from needs_work to needs_review

The patch needs review.

comment:13 Changed 7 years ago by eviatarbach

dsm, would you object to having a more precise overestimation, such as 3.32192809488736235?

comment:14 Changed 7 years ago by dsm

I don't really object, although putting the overestimate right at the end might cause problems for people with more digits than they can possibly have memory for, as it's hard to trust arithmetic down at the lsb level unless you think harder about things than this problem warrants. For example, float("3.32192809488736235") < math.log(10,2), even though 3.32192809488736235 > log(10,2), which might not be what you'd guess.

Anyway, I'm fine with adding more digits, as long as it's still an overestimate w.r.t. math.log. If 3.321928094887363 is okay, then I'll update the patch, fix the docs, and we can put this to bed.

comment:15 Changed 7 years ago by eviatarbach

Sounds good!

Changed 7 years ago by dsm

v3: use a more precise approx

comment:16 Changed 7 years ago by dsm

Okay; anyone want to have a look and make sure I didn't break anything while changing it?

comment:17 Changed 7 years ago by dsm

Oh, yeah. And unfortunately putting in the test from trac #12163 doubles the runtime of testing functional.py from ~9s to ~22s for me. I asked on IRC and the response was that we could probably get away without marking it up as a long test, but that's easily changed.

comment:18 Changed 7 years ago by kcrisman

I don't see why you couldn't just mark the second (huge) test as # long time. That seems like a decent compromise.

comment:19 Changed 7 years ago by davidloeffler

Apply trac_10164_round_up_log10_2_v3.patch

(for patchbot)

comment:20 Changed 7 years ago by davidloeffler

  • Description modified (diff)

comment:21 Changed 7 years ago by benjaminfjones

  • Reviewers changed from Karl-Dieter Crisman to Karl-Dieter Crisman, Benjamin Jones

The v3 patch looks very good. I added a reviewer patch that defined a constant LOG_TEN_TWO_PLUS_EPSILON in order to name the magic number used to overestimate log(10, 2). I'm running tests now.

comment:22 Changed 7 years ago by eviatarbach

Hello,

Great that this patch is moving forward!

Do you think you could just add a comment on the tests explaining that the decimal point adds an extra 1 to the length? Otherwise it seems inaccurate until further inspection.

Thanks.

Changed 7 years ago by benjaminfjones

name the magic 3.23... number, remove trailing whitespace

comment:23 Changed 7 years ago by benjaminfjones

Thanks eviatar, I added the comment to the docstring.

I ran doc tests on all files that were touched and those passed, no problems. I'll wait on the patchbot's approval before giving positive review.

Patchbot, apply:

  1. trac_10164_round_up_log10_2_v3.patch
  2. trac_10164_reviewer.patch
Last edited 7 years ago by benjaminfjones (previous) (diff)

Changed 7 years ago by benjaminfjones

last two patches folded togethers with whitespace removed

comment:24 Changed 7 years ago by benjaminfjones

  • Description modified (diff)

comment:25 Changed 7 years ago by benjaminfjones

Patchbot is refusing to cooperate. I'm running make ptestlong manually.

comment:26 Changed 7 years ago by benjaminfjones

  • Status changed from needs_review to positive_review

All tests pass, positive review.

comment:27 Changed 7 years ago by jdemeyer

  • Merged in set to sage-5.1.beta3
  • Resolution set to fixed
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.