Opened 5 years ago

Closed 3 years ago

#18210 closed defect (invalid)

numerical bug in incomplete gamma function

Reported by: VivianePons Owned by:
Priority: major Milestone: sage-duplicate/invalid/wontfix
Component: symbolics Keywords: sd67
Cc: Merged in:
Authors: Reviewers:
Report Upstream: Fixed upstream, in a later stable release. Work issues:
Branch: Commit:
Dependencies: Stopgaps:

Description (last modified by rws)

Wrong values from incomplete gamma:

$ ./sage -v
SageMath Version 6.7.beta1, Release Date: 2015-04-15

$ ./sage
sage: gamma(60, 30).numerical_approx()
-1.28306738270893e74

Reason is a bug in Pari:

? \p 50
   realprecision = 57 significant digits (50 digits displayed)
? incgam(60,30)
%2 = 1.3868299023788799504133735635863795825413935892945 E80
? \p 30
   realprecision = 38 significant digits (30 digits displayed)
? incgam(60,30)
%3 = -1.23084064495468737276106696496 E74

Previous description:

The following code:

sage: var('R k')
(R, k)
sage: integrated = -(gamma(1/2*k, 1/2*R*k) - gamma(1/2*k))/gamma(1/2*k)
sage: plot(integrated.subs(k=41), (R,0,6))

gives a Seg fault error and crashes Sage entirely:

line 134: 4216 Segmentation fault (core dumped)

I believe it comes from the symbolic computation because the following code is working:

sage: var('R k')
(R, k)
sage: integrated = -(gamma(1/2*k, 1/2*R*k) - gamma(1/2*k))/gamma(1/2*k)
sage: plot(integrated.subs(k=41.), (R,0,6))

(Note the 41. instead of 41)

I have no idea where this comes from but it's really bad, because whatever is going wrong, crashing Sage is bad!

Attachments (6)

working-graph-2012.png (285.4 KB) - added by buck 5 years ago.
broken-graph-2015.png (51.1 KB) - added by buck 5 years ago.
unstable incomplete gamma.ipynb.json (360.8 KB) - added by buck 5 years ago.
demo of the gamma function's instability
unstable incomplete gamma.pdf (111.9 KB) - added by buck 4 years ago.
unstable incomplete gamma.ipynb (360.8 KB) - added by buck 4 years ago.
broken-graph-2016.png (56.5 KB) - added by buck 4 years ago.

Download all attachments as: .zip

Change History (47)

comment:1 Changed 5 years ago by VivianePons

More funny stuff:

sage: a = 319830986772877770815625*sqrt(pi)
sage: b = 1048576*gamma(41/2, 41/2*x)
sage: f = a - b
sage: plot(f, (R,0,1))  # works fine
sage: plot(f / sqrt(pi), (R,0,1))  # BOOM

comment:2 Changed 5 years ago by nbruin

This is recent because on 6.5beta1 the error doesn't occur.

It looks like some infinite recursion in ginac. Excerpt from traceback obtained from sage -gdb after segv has occurred:

...
#12573 0x00007fffc8185e3d in GiNaC::mul::eval(int) const ()
   from /usr/local/sage/sage-git/local/lib/libpynac.so.1
#12574 0x00007fffc80f170b in GiNaC::ex::construct_from_basic(GiNaC::basic const&) () from /usr/local/sage/sage-git/local/lib/libpynac.so.1
#12575 0x00007fffc80acb59 in GiNaC::ex::ex(GiNaC::basic const&) ()
   from /usr/local/sage/sage-git/local/lib/libpynac.so.1
#12576 0x00007fffc8185e3d in GiNaC::mul::eval(int) const ()
   from /usr/local/sage/sage-git/local/lib/libpynac.so.1
#12577 0x00007fffc80f170b in GiNaC::ex::construct_from_basic(GiNaC::basic const&) () from /usr/local/sage/sage-git/local/lib/libpynac.so.1
#12578 0x00007fffc80acb59 in GiNaC::ex::ex(GiNaC::basic const&) ()
   from /usr/local/sage/sage-git/local/lib/libpynac.so.1
...

Changed 5 years ago by buck

Changed 5 years ago by buck

comment:4 Changed 5 years ago by buck

I think there's more to it than that even. Using the limited-precision numerics, I get this plot:



In 2012, when the exact numerics still worked, I got this graph:



comment:5 Changed 5 years ago by kcrisman

  • Type changed from PLEASE CHANGE to defect

comment:6 follow-up: Changed 5 years ago by rws

Interestingly, neither the first snippet in the description nor the BOOM snippet in comment:1 will crash on 6.7beta1/OpenSuSE. What machine/system is that, Viviane/buck/nbruin?

comment:7 Changed 5 years ago by buck

I'm on OS X 10.10.3

comment:8 Changed 5 years ago by VivianePons

Right, I just tried on 6.7 and indeed it is not crashing!

Buck I'd be happy to know if you obtain plot that "make sense" on 6.7. Sorry for the extra compiling, to get 6.7, just checkout the develop ranch and do a git pull. Let me know if you have problems.

comment:9 in reply to: ↑ 6 Changed 5 years ago by nbruin

Replying to rws:

Interestingly, neither the first snippet in the description nor the BOOM snippet in comment:1 will crash on 6.7beta1/OpenSuSE. What machine/system is that, Viviane/buck/nbruin?

(on fedora) 6.7beta1 doesn't segfault anymore either. I haven't checked if the precision issues remain.

comment:10 Changed 5 years ago by rws

No crash on OSX/6.7.beta1 either. Any crashes on 6.7beta1 at all?

comment:11 follow-up: Changed 5 years ago by buck

I also get no segfault in 6.7, but the numerical issue persists. Separate issue I suppose.

comment:12 in reply to: ↑ 11 Changed 5 years ago by rws

Replying to buck:

I also get no segfault in 6.7, but the numerical issue persists. Separate issue I suppose.

Since no one can reproduce the crash on a recent Sage we will use this ticket for the numeric issue. I tried the code in the sage-devel posting but had no glitches in the graph. Can you please give a fully contained code snipped that shows the glitches on a recent Sage system?

comment:13 Changed 5 years ago by buck

@rws: This is a reposting of my sage-support thread here where I gave some small reproductions of the problem.


The incomplete gamma function is resulting in negative outputs with positive real input, which shouldn't happen.

sage: gamma(60, 30).numerical_approx()
-1.28306738270893e74

Attached is a notebook showing the same in graphical form.

If I feed in numbers of a higher precision, I get the correct result:

sage: R10 = RealField(200)
sage: gamma(R10(60), R10(30)).numerical_approx()
1.38682990237888e80

My final problem is that there doesn't seem to be any way to get higher precision numbers to feed all the way through the plot function:

sage: show(plot3d(log(gamma(k, k*x)), (x, R10(0), R10(5/2)), (k, R10(1), R10(101))))
Launched jmol viewer for Graphics3d Object

(you would see a 3d graph with a bite out of it where the precision issue exists)


I'm able to reproduce it today on the latest sage:

$ ./sage -v
SageMath Version 6.7.beta1, Release Date: 2015-04-15

$ ./sage
sage: gamma(60, 30).numerical_approx()
-1.28306738270893e74

Changed 5 years ago by buck

demo of the gamma function's instability

comment:14 Changed 5 years ago by rws

  • Description modified (diff)
  • Summary changed from Symbolic computation of Gamma related function crashes Sage to numerical bug in incomplete gamma function

comment:15 follow-up: Changed 5 years ago by rws

  • Report Upstream changed from N/A to Not yet reported upstream; Will do shortly.

BTW this will be fixed with #16697 (needs review) because there we switched to mpmath for numerics. You can already use this in Sage via:

sage: from sage.libs.mpmath import utils as mpmath_utils
sage: import mpmath
sage: mpmath_utils.call(mpmath.gammainc, 60, 0, 30)
1.28307801824120e74

This ticket may be closed therefore (after the Pari devs have ack'ed the bug) as duplicate of #16697 (augmenting its description with the bug).

comment:16 in reply to: ↑ 15 Changed 5 years ago by rws

Replying to rws:

sage: from sage.libs.mpmath import utils as mpmath_utils
sage: import mpmath
sage: mpmath_utils.call(mpmath.gammainc, 60, 0, 30)
1.28307801824120e74

Ah no, that is wrong as well. I'll report to Johansson too.

comment:17 Changed 5 years ago by rws

  • Report Upstream changed from Not yet reported upstream; Will do shortly. to Reported upstream. No feedback yet.

comment:18 Changed 5 years ago by rws

No, mpmath is correct. I confused lower/upper. So, as said, use the #16697 branch or better, review it:

sage: gamma(60,30).n()
1.38682990237888e80

comment:19 Changed 4 years ago by rws

  • Report Upstream changed from Reported upstream. No feedback yet. to Reported upstream. Developers acknowledge bug.

I was wrong referring to #16697 because this is a different bug. This issue is still open upstream.

comment:20 Changed 4 years ago by buck

rws, could you please copy the relevant bits from #16697?

I'm going to work to close this in coming weeks.

comment:21 Changed 4 years ago by buck

Discusion migrated from #16697


by buck:

... this patch doesn't fix the numerical instability in incomplete gamma as promised elsewhere.

sage: gamma_inc_lower(25, 14.5).n()
-6.63538851954338e22
sage: gamma_inc_lower(25, 14.5).n(algorithm='mpmath')
-6.63538851954338e22
sage: gamma_inc_lower(25, 14.5).n(algorithm='pari')
-6.63538851954338e22

by buck:

Of note, gamma (even incomplete) is strictly positive in the reals; it's an integral of a strictly positive function.

https://en.wikipedia.org/wiki/Incomplete_gamma_function#Definition

by buck:

Wolfram gives the value of the same as 4.7e21, which seems correct. ​http://www.wolframalpha.com/input/?i=Gamma%5B25%2C+0%2C+14.5%5D

by buck:

Actually, if I ask pari directly, it gets the correct value, so maybe the algorithm argument is broken?

sage: pari('gamma(25) - incgam(25, 14.5)')
4.71197173824555 E21
sage: pari('incgamc(25, 14.5)')
4.71197173824555 E21

by buck:

Same story if I ask mpmath directly. I can't imagine what's wrong.

sage: call(mpmath.gammainc, 25, 0, 14.5)
4.71197173824555e21

by rws:

Replying to buck:

Same story if I ask mpmath directly. I can't imagine what's wrong.

'algorithm' isn't propagated to evalf so we always get the buggy pari values. As we would have to switch the default anyway as soon as a later pari version fixes the bug, I'll now simply disable pari. We'll have to live with mpmath's 53 bits for now.

by rws:

Replying to buck:

Actually, if I ask pari directly, it gets the correct value, so maybe the algorithm argument is broken?

sage: pari('gamma(25) - incgam(25, 14.5)')
4.71197173824555 E21
sage: pari('incgamc(25, 14.5)')
4.71197173824555 E21

This works because it uses the gp interface (the libpari might be broken on our side). Ah, that's two different bugs. However, even if it is in our interface to libpari and we fix this, or if we use the gp interface, there is still the unfixed pari bug at ​http://pari.math.u-bordeaux.fr/cgi-bin/bugreport.cgi?bug=1689PARI/GP

by rws:

Actually the 53-bit restriction applies only to larger values, so all is not lost.

comment:22 follow-up: Changed 4 years ago by buck

@rws: Could you expand on the "two different bugs" and link them (and/or create tickets for them)?

Also, could you explain why "because it uses the gp interface" would obscure the negative-number return value?

comment:23 in reply to: ↑ 22 Changed 4 years ago by rws

Replying to buck:

@rws: Could you expand on the "two different bugs" and link them (and/or create tickets for them)?

Okay, this ticket is the (60,30) bug from comment:13, i.e. gamma(60, 30).numerical_approx(). The other one is the fact that gamma(25,14.5) and CC(25).gamma_inc(14.5), which both use libpari, give

sage: CC(25).gamma_inc(14.5)
6.86802286928673e23
sage: gamma(25,14.5)
6.86802286928673e23

while

sage: pari('incgam(25, 14.5)')
6.15736429994994 E23

But then, we see also

sage: pari('incgam(60,30)')
-1.23084064495469 E74
sage: gamma(60,30).n()
-1.28306738270893e74

so there is a mismatch between the two interfaces, and the gp interface is sometimes right where libpari is wrong.

Also, could you explain why "because it uses the gp interface" would obscure the negative-number return value?

Because the gp interface is sometimes right, apparently.

However, I don't think it useful to open two tickets as long as we don't know if there really are two different bugs.

comment:24 Changed 4 years ago by buck

The bug reported in the OP seems gone; I can't reproduce the segfault. Can we change the topic of the bug to strictly about the numerical issue? That seems to be what most of the discussion was about.

comment:25 Changed 4 years ago by buck

@rws: I'll spend some time digging into this myself, and produce a pari and/or gp patch if I can, but I'm a bit lost as to where to start. Could you give some clues as to how I might pointpoint the issue, what code I should look at patching?

comment:26 Changed 4 years ago by rws

In my opinion the Pari bug (which shows itself when using the GP/Pari program without Sage) should be fixed first. This means looking into the Pari code and reporting to its developers on that ticket already given in comment:17. After that it remains to be seen what's still wrong. For completeness Sage's libpari interface is in libs/pari.

comment:27 Changed 4 years ago by rws

Karim Belabas of Pari just sent:

I believe the problem is now fixed in master by the following commit
(rewritten from Henri Cohen's proposed patch)

commit 4276cc667123a92373c4aff7195eb728dd39f6e1
Author: Karim Belabas <Karim.Belabas@math.u-bordeaux1.fr>
Date:   Fri Jan 15 15:58:17 2016 +0100

    HC140- incgam(30,60) < 0. More generally, wrong results for s >> 1 [#1689]

Sorry it took so long. Please test and report problems !

comment:28 Changed 4 years ago by buck

Thanks! Building master sagemath branch now...

comment:29 Changed 4 years ago by buck

I still get the same broken graph in sage master.

Sage is at a1b60f2c3e06171caa4aa044a73ea58236686d30 and pari is at 'GP/PARI CALCULATOR Version 2.8.0 (development git-6157df4)'.

I can't tell immediately if that's before or after the referenced pari commit above.

comment:30 Changed 4 years ago by buck

Ah after cloning pari, it's apparent that the installed sage version is four months behind the fixing commit.

I'll try to find in the docs how to get sage to install this newer version.

comment:31 Changed 4 years ago by buck

Procedure:

cd $PARI_ROOT
TAG=$(git describe)     
git archive HEAD -o $TAG.tar.gz --prefix=$TAG/
mv pari-2.8-2230-g450ce38.tar.gz $SAGE_ROOT/upstream/ 
cd $SAGE_ROOT
./sage --fix-pkg-checksums
./sage -i -f -c pari

Changes to sage: (I had to touch up one of the patches to get it to apply.)

  • build/pkgs/pari/checksums.ini

    diff --git a/build/pkgs/pari/checksums.ini b/build/pkgs/pari/checksums.ini
    index c62c530..aabe08d 100644
    a b  
    11tarball=pari-VERSION.tar.gz
    2 sha1=fa23e0c8b6e38a356048d19224dc9b9658d77724
    3 md5=c753faaa4780de5ad8d461f0ffd70ecf
    4 cksum=1220765312
     2sha1=55299bfe042da491648897e830030083809d9967
     3md5=97738f8e92bba498f7dfd723c8e9d910
     4cksum=1221580104
  • build/pkgs/pari/package-version.txt

    diff --git a/build/pkgs/pari/package-version.txt b/build/pkgs/pari/package-version.txt
    index 2b25bd1..5fdd71e 100644
    a b  
    1 2.8-1813-g6157df4.p0
     12.8-2230-g450ce38
  • build/pkgs/pari/patches/public_memory_functions.patch

    diff --git a/build/pkgs/pari/patches/public_memory_functions.patch b/build/pkgs/pari/patches/public_memory_functions.patch
    index b3726d7..ee14fa4 100644
    a b index 7067183..4ede6ed 100644 
    99+void *  pari_mainstack_malloc(size_t size);
    1010+void    pari_mainstack_mfree(void *s, size_t size);
    1111+void    pari_mainstack_free(struct pari_mainstack *st);
    12  void    paristack_alloc(size_t rsize, size_t vsize);
    1312 void    paristack_newrsize(ulong newsize);
    1413 void    paristack_resize(ulong newsize);
     14 void    paristack_setsize(size_t rsize, size_t vsize);
    1515diff --git a/src/language/init.c b/src/language/init.c
    1616index 7b5922d..2a578d7 100644
    1717--- a/src/language/init.c

I also had to install bison to get the build to run, which I hadn't installed before, but I *had* built pari before (I think!) so I believe that's a new build dependency there?

Result: Successfully installed pari-2.8-2230-g450ce38; The PARI self-tests all passed

The sage console fails to start however, with ImportError: sage/local/lib/python2.7/site-packages/sage/libs/pari/gen.so: undefined symbol: gp_algcenter

Last edited 4 years ago by buck (previous) (diff)

comment:32 Changed 4 years ago by buck

Subsequent make fails during pari autoregen:

python -c "from sage_setup.autogen.pari import rebuild; rebuild()"
Generating PARI functions: (!_) (#_) (%) (%#) (+_) (-_) Catalan Col Colrev (DEBUGLEVEL) Euler I List Map Mat Mod (O) 
...
galoispermtopol galoissubcyclo galoissubfields galoissubgroups gamma gammah gammamellininvTraceback (most recent call last):
  File "<string>", line 1, in <module>
  File "sage_setup/autogen/pari/__init__.py", line 5, in rebuild
    G()
  File "sage_setup/autogen/pari/generator.py", line 247, in __call__
    self.handle_pari_function(**v)
  File "sage_setup/autogen/pari/generator.py", line 173, in handle_pari_function
    args, ret = parse_prototype(prototype, help)
  File "sage_setup/autogen/pari/parser.py", line 204, in parse_prototype
    raise ValueError('unknown prototype character %r' % c)
ValueError: unknown prototype character 'b'

I note that the build gets into a weird state afterward because the file being generated (src/sage/libs/pari/auto_gen.pxi) is written to even though its generation failed. We could avoid this state by making a patch like so:

  • src/sage_setup/autogen/pari/generator.py

    diff --git a/src/sage_setup/autogen/pari/generator.py b/src/sage_setup/autogen/pari/generator.py
    index 7aa9990..9b78d88 100644
    a b class PariFunctionGenerator(object): 
    233233        D = sorted(D.values(), key=lambda d: d['function'])
    234234        sys.stdout.write("Generating PARI functions:")
    235235
    236         self.gen_file = open(self.gen_filename, 'w')
     236        self.gen_file = open(self.gen_filename + '.tmp', 'w')
    237237        self.gen_file.write(gen_banner)
    238         self.instance_file = open(self.instance_filename, 'w')
     238        self.instance_file = open(self.instance_filename + '.tmp', 'w')
    239239        self.instance_file.write(instance_banner)
    240240
    241241        for v in D:
    class PariFunctionGenerator(object): 
    249249
    250250        self.gen_file.close()
    251251        self.instance_file.close()
     252
     253        # All done? Let's commit.
     254        os.rename(self.gen_filename + '.tmp', self.gen_filename)
     255        os.rename(self.instance_filename + '.tmp', self.instance_filename)

I'm out of my depth at this point. I think I need help from someone that actually know what pari is at all =X

comment:33 Changed 4 years ago by buck

I believe I've fixed the ValueError: unknown prototype character 'b' issue. Pari recently added the concept of "bit precision", which (I believe) is the natural precision format in sage, so needs no conversion:

  • src/sage_setup/autogen/pari/args.py

    diff --git a/src/sage_setup/autogen/pari/args.py b/src/sage_setup/autogen/pari/args.py
    index 57356b4..a2749df 100644
    a b class PariArgumentPrec(PariArgumentClass): 
    254254        s = "        {name} = prec_bits_to_words({name})\n"
    255255        return s.format(name=self.name)
    256256
     257class PariArgumentBitPrec(PariArgumentClass):
     258    def _typerepr(self):
     259        return "bitprec"
     260    def ctype(self):
     261        return "long"
     262    def always_default(self):
     263        return "0"
     264    def get_argument_name(self, namesiter):
     265        return "bitprecision"
     266
    257267class PariArgumentSeriesPrec(PariArgumentClass):
    258268    def _typerepr(self):
    259269        return "serprec"
    pari_arg_types = { 
    277287        'U': PariArgumentULong,
    278288        'n': PariArgumentVariable,
    279289        'p': PariArgumentPrec,
     290        'b': PariArgumentBitPrec,
    280291        'P': PariArgumentSeriesPrec,
    281292
    282293    # Codes which are known but not actually supported for Sage

There's again a further issue that I don't understand yet:

Compiling sage/libs/pari/pari_instance.pyx because it depends on sage/libs/pari/auto_instance.pxi.
[1/2] Cythonizing sage/libs/pari/gen.pyx

Error compiling Cython file:
------------------------------------------------------------
...
            
            [0 0 -2 -1]
        """
        cdef GEN _al = al.g
        pari_catch_sig_on()
        cdef GEN _ret = algmultable(_al)
                                  ^
------------------------------------------------------------

sage/libs/pari/auto_gen.pxi:1558:35: Call with wrong number of arguments (expected 2, got 1)

On brief inspection, I don't see that the number of arguments should have changed.

comment:34 Changed 4 years ago by rws

While this is great work I believe you should open a separate Pari upgrade ticket for this.

comment:35 Changed 4 years ago by buck

@rws Yes I didn't anticipate this would be such work. Will do.

comment:36 Changed 4 years ago by buck

Pari upgrade work migrated to #19905.

comment:37 follow-up: Changed 4 years ago by buck

I was finally able to integrate pari master-branch with sage, and while the problem is much improved, it's not fixed. Here's the graph from the original bug report, using newest pari:



Those discontinuities are false. An example bad value:

sage: gamma(100, 7.01)
-2.71843697211100e151

I will report this upstream, and attach a worksheet showing more detail.

Last edited 4 years ago by buck (previous) (diff)

Changed 4 years ago by buck

Changed 4 years ago by buck

Changed 4 years ago by buck

comment:38 in reply to: ↑ 37 ; follow-up: Changed 3 years ago by slelievre

This ticket seems to be solved now in Sage 7.3.beta9.

Trying the code from the ticket description...

$ ./sage -v
SageMath version 7.3.beta9, Release Date: 2016-07-22
$ ./sage -q
sage: gamma(60, 30).numerical_approx()
1.38682990237888e80
sage: %gp

  --> Switching to PARI/GP interpreter <--

pari: \p 50
realprecision = 57 significant digits (50 digits displayed)
pari: incgam(60,30)
1.3868299023788801161747839921239559320835799009004 E80
pari: \p 30
realprecision = 38 significant digits (30 digits displayed)
pari: incgam(60,30)
1.38682990237888011617478399212 E80
pari: <ctrl-D>

  --> Exiting back to Sage <--

sage: _ = var('R k')
sage: integrated = -(gamma(1/2*k, 1/2*R*k) - gamma(1/2*k))/gamma(1/2*k)
sage: plot(integrated.subs(k=41), (R,0,6))
Launched png viewer for Graphics object consisting of 1 graphics primitive
sage: integrated = -(gamma(1/2*k, 1/2*R*k) - gamma(1/2*k))/gamma(1/2*k)
sage: plot(integrated.subs(k=41.), (R,0,6))
Launched png viewer for Graphics object consisting of 2 graphics primitives
sage: from sage.libs.mpmath import utils as mpmath_utils
sage: import mpmath
sage: mpmath_utils.call(mpmath.gammainc, 60, 0, 30)
1.28307801824120e74

Replying to buck:

I was finally able to integrate pari master-branch with sage, and while the problem is much improved, it's not fixed. Here's the graph from the original bug report, using newest pari:

Trying that now seems to work for me, I get nice graphs with no spikes.

sage: sum(
....:     plot(
....:         integrated(k=k),
....:         (x, 0, 2.5),
....:         ymax=1.6,
....:         color=color,
....:         legend_label='k=%i' % k,
....:         figsize=6,
....:         aspect_ratio=1.0,
....:     )
....:     for k, color in zip(
....:         [21, 22, 23, 26, 31, 41],
....:         ['blue', 'brown', 'red', 'green', 'black', 'yellow', 'orange'],
....:     )
....: )
Launched png viewer for Graphics object consisting of 6 graphics primitives

Those discontinuities are false. An example bad value:

sage: gamma(100, 7.01)
-2.71843697211100e151

I will report this upstream, and attach a worksheet showing more detail.

What is the correct value? I get the following in Sage 7.3.beta9.

sage: gamma(100, 7.01)
9.33262154439442e155

comment:39 in reply to: ↑ 38 Changed 3 years ago by rws

  • Milestone changed from sage-6.7 to sage-duplicate/invalid/wontfix
  • Report Upstream changed from Reported upstream. Developers acknowledge bug. to Fixed upstream, in a later stable release.
  • Status changed from new to needs_review

Replying to slelievre:

This ticket seems to be solved now in Sage 7.3.beta9.

Agree.

What is the correct value? I get the following in Sage 7.3.beta9.

sage: gamma(100, 7.01)
9.33262154439442e155
sage: ComplexBallField(200)(100).gamma(ComplexBallField(200)(701/100))
[9.3326215443944152681699238856266700490715968264381621468593e+155 +/- 4.56e+96]

comment:40 Changed 3 years ago by slelievre

  • Status changed from needs_review to positive_review

comment:41 Changed 3 years ago by vbraun

  • Resolution set to invalid
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.