Opened 12 years ago

Closed 7 years ago

#7099 closed defect (fixed)

serious incomplete gamma function precision bugs

Reported by: was Owned by: jkantor
Priority: critical Milestone: sage-6.4
Component: numerical Keywords:
Cc: kcrisman Merged in:
Authors: Peter Bruin Reviewers: Ralf Stephan
Report Upstream: N/A Work issues:
Branch: 50b629f (Commits, GitHub, GitLab) Commit: 50b629f157f8547661370c39f98137efd46bab67
Dependencies: Stopgaps:

Status badges

Description

sage: C = ComplexField(1000)
sage: C(2+I).gamma_inc(C(3+I))
0.1215156446645086956511068454478419198494520969688892364501953125000000\
000000000000000000000000000000000000000000000000000000000000000000000000\
000000000000000000000000000000000000000000000000000000000000000000000000\
000000000000000000000000000000000000000000000000000000000000000000000000\
00000000000000 +
0.1015339090798260332775427433604775728781532961875200271606445312500000\
000000000000000000000000000000000000000000000000000000000000000000000000\
000000000000000000000000000000000000000000000000000000000000000000000000\
000000000000000000000000000000000000000000000000000000000000000000000000\
00000000000000*I

What's with all the zeros? Same bad behavior for higher precision and any other random input. This is undoubtedly the typical mistake of not setting the pari precision correctly, etc., etc.

Change History (25)

comment:1 Changed 11 years ago by kcrisman

  • Cc kcrisman added
  • Report Upstream set to N/A

comment:2 Changed 10 years ago by kcrisman

Correct. Here's the problem, I think.

sage: C = ComplexField(1000)
sage: c = C(2+I)
sage: c
2.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000 + 1.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*I
sage: d = c.real()
sage: d
2.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
sage: d._pari_()
2.00000000000000

So here is the problem, right in the change to Pari from the RealField.

comment:3 Changed 8 years ago by jdemeyer

  • Milestone changed from sage-5.11 to sage-5.12

comment:4 Changed 8 years ago by vbraun_spam

  • Milestone changed from sage-6.1 to sage-6.2

comment:5 Changed 8 years ago by pbruin

A probably related bug was just reported on sage-support:

sage: numerical_approx(gamma(9, 10^(-3))-gamma(9), digits=40)
0.000000000000000

The correct answer is approximately -1.1101115655 * 10-28.

Last edited 8 years ago by pbruin (previous) (diff)

comment:6 Changed 8 years ago by pbruin

  • Authors set to Peter Bruin
  • Branch set to u/pbruin/7099-incomplete_gamma
  • Commit set to 3504d34ae7027b90c66694aa22c3351d44cf37e5
  • Status changed from new to needs_info

Here is a somewhat rough fix. It doesn't take into account the precisions of both input parameters separately, but uses the precision of the desired parent in Function_gamma_inc, and the precision of the parent of the first argument in ComplexNumber.gamma_inc() et al. Maybe an expert on this can say whether this behaviour is consistent with other special functions in Sage?

(Note that the relatively low precision in the new doctest is because of cancellation.)

comment:7 Changed 8 years ago by git

  • Commit changed from 3504d34ae7027b90c66694aa22c3351d44cf37e5 to f675b0d8c50c5e34882116784161302659920dd8

Branch pushed to git repo; I updated commit sha1. New commits:

f675b0dTrac 7099: add doctest

comment:8 Changed 8 years ago by jdemeyer

I'm wondering if this a bug in the gamma function or in symbolic functions in general...

comment:9 follow-up: Changed 8 years ago by jdemeyer

Isn't this the only and actual bug? That the following doesn't depend on the prec parameter.

sage: numerical_approx(gamma(9,10^-3), prec=1024)                                    
40320.0000000000

comment:10 Changed 8 years ago by jdemeyer

I don't know whether the change to src/sage/functions/other.py is right (I'm not saying it is wrong, just don't know), but the adding of the precision arguments is obviously needed.

comment:11 Changed 8 years ago by git

  • Commit changed from f675b0d8c50c5e34882116784161302659920dd8 to 78ba9638700d8cd9d4b18d83323af8fe9507a7ba

Branch pushed to git repo; I updated commit sha1. New commits:

0117c69Trac 7099: one new doctest, move another one
78ba963Trac 7099: small simplifications

comment:12 follow-up: Changed 8 years ago by pbruin

  • Priority changed from major to critical
  • Status changed from needs_info to needs_review

Setting priority to critical because the bug leads to mathematically incorrect answers.

The original bug is already solved by passing the correct precision to PARI's incgam().

The bug in comment:5 should be solved by the additional changes to Function_gamma_inc._evalf_().

comment:13 in reply to: ↑ 12 Changed 8 years ago by jdemeyer

Replying to pbruin:

The bug in comment:5 should be solved by the additional changes to Function_gamma_inc._evalf_().

I think it is better not to confuse these 2 different bugs.

comment:14 in reply to: ↑ 9 ; follow-up: Changed 8 years ago by dimpase

Replying to jdemeyer:

Isn't this the only and actual bug? That the following doesn't depend on the prec parameter.

sage: numerical_approx(gamma(9,10^-3), prec=1024)                                    
40320.0000000000

is prec parameter mentioned in the doc on gamma? It doesn't, IMHO. Are there other transcendentals for which this parameter does work?

comment:16 in reply to: ↑ 14 Changed 8 years ago by kcrisman

Isn't this the only and actual bug? That the following doesn't depend on the prec parameter.

sage: numerical_approx(gamma(9,10^-3), prec=1024)                                    
40320.0000000000

is prec parameter mentioned in the doc on gamma? It doesn't, IMHO. Are there other transcendentals for which this parameter does work?

I was pretty sure we were deprecating the prec keyword inside such functions. But this is apparently a keyword of numerical_approx, which is different. In which case it should just pass the precision on correctly somewhere else.

comment:17 Changed 8 years ago by vbraun_spam

  • Milestone changed from sage-6.2 to sage-6.3

comment:18 Changed 7 years ago by rws

With #16697 the output is as expected:

sage: C = ComplexField(1000)
sage: gamma_inc(C(2+I),C(3+I))
0.121515644664508695525971545977439666159749344176962379708992904126499444842886620664991650378432544392118359044438541514683402245033018771644222346410367471459456844674335147722343580581945662693850674590491020834632434082710800093315646442975240326569517738365018117780134100101636704042869033248174 + 0.101533909079826033296475736021224621546966200987295663190553587086145836461236284668967411665020429964946098113930918849948956425984499549094441904693395768367238320065064071027383069839637218088862214571990869510941211277488169032567679631037683814516738122300220474252081775895835843619616213883517*I

This is because I changed the default to mpmath for the gamma functions. The ticket is still relevant for C(2+I).gamma_inc(C(3+I)) however.

comment:19 Changed 7 years ago by rws

  • Status changed from needs_review to positive_review

Patchbot was happy at 6.2beta2 so I guess it's still good. The changes look straightforward, anyway.

comment:20 Changed 7 years ago by rws

  • Reviewers set to Ralf Stephan

comment:21 Changed 7 years ago by vbraun_spam

  • Milestone changed from sage-6.3 to sage-6.4

comment:22 Changed 7 years ago by vbraun

  • Status changed from positive_review to needs_work

On the Arando buildbot:

sage -t --long src/sage/rings/complex_mpc.pyx
**********************************************************************
File "src/sage/rings/complex_mpc.pyx", line 2241, in sage.rings.complex_mpc.MPComplexNumber.gamma_inc
Failed example:
    (1+i).gamma_inc(2 + 3*i)
Expected:
    0.0020969149 - 0.059981914*I
Got:
    0.0020969148 - 0.059981914*I
**********************************************************************
1 item had failures:
   1 of   5 in sage.rings.complex_mpc.MPComplexNumber.gamma_inc
    [375 tests, 1 failure, 0.27 s]
sage -t --long src/sage/rings/complex_number.pyx
**********************************************************************
File "src/sage/rings/complex_number.pyx", line 2023, in sage.rings.complex_number.ComplexNumber.gamma_inc
Failed example:
    (1+i).gamma_inc(2 + 3*i)
Expected:
    0.0020969149 - 0.059981914*I
Got:
    0.0020969148 - 0.059981914*I
**********************************************************************
1 item had failures:
   1 of   8 in sage.rings.complex_number.ComplexNumber.gamma_inc
    [378 tests, 1 failure, 3.12 s]

comment:23 Changed 7 years ago by git

  • Commit changed from 78ba9638700d8cd9d4b18d83323af8fe9507a7ba to 50b629f157f8547661370c39f98137efd46bab67

Branch pushed to git repo; I updated commit sha1. New commits:

50b629fTrac 7099: mark two doctests with "abs tol 2e-10"

comment:24 Changed 7 years ago by pbruin

  • Status changed from needs_work to positive_review

I'm assuming this fix does not absolutely require another review round...

comment:25 Changed 7 years ago by vbraun

  • Branch changed from u/pbruin/7099-incomplete_gamma to 50b629f157f8547661370c39f98137efd46bab67
  • Resolution set to fixed
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.