segfault when plotting gamma function

39 views
Skip to first unread message

Buck Evan

unread,
Apr 15, 2015, 12:55:23 PM4/15/15
to sage-s...@googlegroups.com
I'm trying to plot a function involving the gamma function. It's this one:

-(gamma(1/2*k, 1/2*R*k) - gamma(1/2*k))/gamma(1/2*k)


If I try to plot the function with a k value any greater than 40, I get a segfault:

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'],
    )
)
/Users/buck/trees/mine/sage/sage: line 134: 12283 Segmentation fault: 11  "$SAGE_ROOT/src/bin/sage" "$@"
[I 12:15:10.735 NotebookApp] KernelRestarter: restarting kernel (1/5)

This is extremely similar to the behavior I was getting with the pre-built 6.5 sage, so I think it's an old bug rather than a problem with my build.
In 6.5 there was no output (I suspect the Segmentation Fault message was gobbled somehow), so I never tried to enable / debug a core dump in 6.5.

Inspecting the 6.6 core dump, I see that it used 27528, mostly calling back and forth between two functions.
Smells like stack overflow to me.

(lldb) bt                                                                                                                           [27503/47384]
* thread #1: tid = 0x0000, 0x0000000103588007 libpython2.7.dylib`vgetargskeywords + 23, stop reason = signal SIGSTOP
  * frame #0: 0x0000000103588007 libpython2.7.dylib`vgetargskeywords + 23
    frame #1: 0x000000010358890e libpython2.7.dylib`PyArg_ParseTupleAndKeywords + 174
    frame #2: 0x000000010356c40e libpython2.7.dylib`builtin___import__ + 94
    frame #3: 0x00000001034a377f libpython2.7.dylib`PyObject_CallFunctionObjArgs + 367
    frame #4: 0x000000010774c1c4 coerce.so`__Pyx_Import(name=<unavailable>, from_list=0x0000000117da7560, level=<unavailable>) + 180 at coerce.c:
20339
    frame #5: 0x000000010774d08e coerce.so`__pyx_f_4sage_9structure_6coerce_py_scalar_parent(__pyx_v_py_type=<unavailable>, __pyx_skip_dispatch=<
unavailable>) + 510 at coerce.c:2938
    frame #6: 0x0000000107760318 coerce.so`__pyx_f_4sage_9structure_6coerce_24CoercionModel_cache_maps_canonical_coercion(__pyx_v_self=0x00000001
07648d50, __pyx_v_x=0x00007f8fbfc37238, __pyx_v_y=0x0000000118c5f390, __pyx_skip_dispatch=<unavailable>) + 3128 at coerce.c:9585
    frame #7: 0x0000000107763cfe coerce.so`__pyx_pw_4sage_9structure_6coerce_24CoercionModel_cache_maps_25canonical_coercion + 13 at coerce.c:105
64
    frame #8: 0x0000000107763cf1 coerce.so`__pyx_pw_4sage_9structure_6coerce_24CoercionModel_cache_maps_25canonical_coercion(__pyx_v_self=0x00000
00107648d50, __pyx_args=0x0000000117e5f248, __pyx_kwds=0x0000000000000000) + 65
    frame #9: 0x0000000103574f10 libpython2.7.dylib`PyEval_EvalFrameEx + 26208
    frame #10: 0x0000000103575fb0 libpython2.7.dylib`PyEval_EvalCodeEx + 2144
    frame #11: 0x00000001034dae0c libpython2.7.dylib`function_call + 124
    frame #12: 0x000000010d680381 pynac.so`__Pyx_PyObject_Call(func=0x0000000108f0c410, arg=<unavailable>, kw=PyObject * at scalar) + 65 at pynac
.cpp:26785
    frame #13: 0x000000010d6a76e9 pynac.so`py_gcd(__pyx_v_n=0x00007f8fbfc37238, __pyx_v_k=0x0000000118c5f390) + 521 at pynac.cpp:8416
    frame #14: 0x000000010d7d0f28 libpynac.1.dylib`GiNaC::gcd(GiNaC::numeric const&, GiNaC::numeric const&) + 56
    frame #15: 0x000000010d7ad100 libpynac.1.dylib`GiNaC::add::integer_content() const + 496
    frame #16: 0x000000010d7acef0 libpynac.1.dylib`GiNaC::ex::integer_content() const + 16
    frame #17: 0x000000010d79fe01 libpynac.1.dylib`GiNaC::mul::eval(int) const + 1377
    frame #18: 0x000000010d71ba91 libpynac.1.dylib`GiNaC::ex::construct_from_basic(GiNaC::basic const&) + 113
    frame #19: 0x000000010d7a026b libpynac.1.dylib`GiNaC::mul::eval(int) const + 2507
    frame #20: 0x000000010d71ba91 libpynac.1.dylib`GiNaC::ex::construct_from_basic(GiNaC::basic const&) + 113
    frame #21: 0x000000010d7a026b libpynac.1.dylib`GiNaC::mul::eval(int) const + 2507
    frame #22: 0x000000010d71ba91 libpynac.1.dylib`GiNaC::ex::construct_from_basic(GiNaC::basic const&) + 113
    frame #23: 0x000000010d7a026b libpynac.1.dylib`GiNaC::mul::eval(int) const + 2507
    frame #24: 0x000000010d71ba91 libpynac.1.dylib`GiNaC::ex::construct_from_basic(GiNaC::basic const&) + 113
    frame #25: 0x000000010d7a026b libpynac.1.dylib`GiNaC::mul::eval(int) const + 2507
    frame #26: 0x000000010d71ba91 libpynac.1.dylib`GiNaC::ex::construct_from_basic(GiNaC::basic const&) + 113
    frame #27: 0x000000010d7a026b libpynac.1.dylib`GiNaC::mul::eval(int) const + 2507
    frame #28: 0x000000010d71ba91 libpynac.1.dylib`GiNaC::ex::construct_from_basic(GiNaC::basic const&) + 113
    frame #29: 0x000000010d7a026b libpynac.1.dylib`GiNaC::mul::eval(int) const + 2507
...
    frame #27502: 0x000000010d71ba91 libpynac.1.dylib`GiNaC::ex::construct_from_basic(GiNaC::basic const&) + 113
    frame #27503: 0x000000010d7a026b libpynac.1.dylib`GiNaC::mul::eval(int) const + 2507
    frame #27504: 0x000000010d71ba91 libpynac.1.dylib`GiNaC::ex::construct_from_basic(GiNaC::basic const&) + 113
    frame #27505: 0x000000010d7a026b libpynac.1.dylib`GiNaC::mul::eval(int) const + 2507
    frame #27506: 0x000000010d71ba91 libpynac.1.dylib`GiNaC::ex::construct_from_basic(GiNaC::basic const&) + 113
    frame #27507: 0x000000010d7a026b libpynac.1.dylib`GiNaC::mul::eval(int) const + 2507
    frame #27508: 0x000000010d71ba91 libpynac.1.dylib`GiNaC::ex::construct_from_basic(GiNaC::basic const&) + 113
    frame #27509: 0x000000010d7a026b libpynac.1.dylib`GiNaC::mul::eval(int) const + 2507
    frame #27510: 0x000000010d71ba91 libpynac.1.dylib`GiNaC::ex::construct_from_basic(GiNaC::basic const&) + 113
    frame #27511: 0x000000010d7a026b libpynac.1.dylib`GiNaC::mul::eval(int) const + 2507
    frame #27512: 0x000000010d71ba91 libpynac.1.dylib`GiNaC::ex::construct_from_basic(GiNaC::basic const&) + 113
    frame #27513: 0x000000010d7a026b libpynac.1.dylib`GiNaC::mul::eval(int) const + 2507
    frame #27514: 0x000000010d71ba91 libpynac.1.dylib`GiNaC::ex::construct_from_basic(GiNaC::basic const&) + 113
    frame #27515: 0x000000010d7a026b libpynac.1.dylib`GiNaC::mul::eval(int) const + 2507
    frame #27516: 0x000000010d71ba91 libpynac.1.dylib`GiNaC::ex::construct_from_basic(GiNaC::basic const&) + 113
    frame #27517: 0x000000010d7d1ccd libpynac.1.dylib`GiNaC::operator*(GiNaC::ex const&, GiNaC::ex const&) + 77
    frame #27518: 0x000000010da76768 expression.so`__pyx_f_4sage_8symbolic_10expression_10Expression__mul_(__pyx_v_left=0x0000000117db6cb0, __pyx_v_right=0x0000000117db6ea8, __pyx_skip_dispatch=<unavailable>) + 2504 at expression.cpp:18109
    frame #27519: 0x00000001075bed30 element.so`__pyx_pw_4sage_9structure_7element_11RingElement_11__mul__ + 64 at element.c:17301
    frame #27520: 0x00000001075becf0 element.so`__pyx_pw_4sage_9structure_7element_11RingElement_11__mul__(__pyx_v_left=0x0000000117db6cb0, __pyx_v_right=0x0000000117db6ea8) + 14
    frame #27521: 0x000000010349ff8b libpython2.7.dylib`PyNumber_Multiply + 187
    frame #27522: 0x0000000117bceb47 wrapper_py.so`interp_py(args=0x0000000117dbeda8, constants=0x00007f8fbd980f30, stack=<unavailable>, code=0x00007f8fbd9a8eac) + 519 at interp_py.c:110
    frame #27523: 0x0000000117bc1cca wrapper_py.so`__pyx_pw_4sage_3ext_12interpreters_10wrapper_py_10Wrapper_py_5__call__(__pyx_v_self=0x00000001181db380, __pyx_args=0x0000000117dbed90, __pyx_kwds=0x0000000000000000) + 138 at wrapper_py.c:1478
    frame #27524: 0x00000001034a2bd3 libpython2.7.dylib`PyObject_Call + 67
    frame #27525: 0x00000001035727bd libpython2.7.dylib`PyEval_EvalFrameEx + 16141
    frame #27526: 0x0000000103575fb0 libpython2.7.dylib`PyEval_EvalCodeEx + 2144
    frame #27527: 0x000000010357513d libpython2.7.dylib`PyEval_EvalFrameEx + 26765
    frame #27528: 0x000000010e399354

Volker Braun

unread,
Apr 15, 2015, 1:04:17 PM4/15/15
to sage-s...@googlegroups.com
NameError: global name 'integrated' is not defined

Buck Evan

unread,
Apr 15, 2015, 1:56:57 PM4/15/15
to sage-s...@googlegroups.com
It's the function I mentioned:

sage: var('R k')
(R, k)
sage: integrated = -(gamma(1/2*k, 1/2*R*k) - gamma(1/2*k))/gamma(1/2*k)

Vincent Delecroix

unread,
Apr 15, 2015, 2:04:15 PM4/15/15
to sage-s...@googlegroups.com
Already the following fails

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))
/home/vincent/bin/sage: line 134: 13683 Segmentation fault
"$SAGE_ROOT/src/bin/sage" "$@"

Viviane Pons

unread,
Apr 15, 2015, 2:08:28 PM4/15/15
to sage-s...@googlegroups.com
More context about this:

I'm here with Buck at PyCon and he said he actually had this line running a few years ago. I don't know anything about this stuff so it's difficult to track what's wrong now.

But it's really bad, as you can see, it's not just a usual exception, it's actually crashing Sage...

(also, this is on Sage 6.6 and I believe it was there before. Buck had the error on the previous version and updated it to see if was fixed on 6.6)

Best

Viviane

--
You received this message because you are subscribed to the Google Groups "sage-support" group.
To unsubscribe from this group and stop receiving emails from it, send an email to sage-support+unsubscribe@googlegroups.com.
To post to this group, send email to sage-s...@googlegroups.com.
Visit this group at http://groups.google.com/group/sage-support.
For more options, visit https://groups.google.com/d/optout.

Vincent Delecroix

unread,
Apr 15, 2015, 2:09:40 PM4/15/15
to sage-s...@googlegroups.com
Do you know how to open a ticket?
>> email to sage-support...@googlegroups.com.

Viviane Pons

unread,
Apr 15, 2015, 2:12:28 PM4/15/15
to sage-s...@googlegroups.com
I'll do it, but I was hopping to get more context about it, also checking if it wasn't known somewhere.


To post to this group, send email to sage-s...@googlegroups.com.
Visit this group at http://groups.google.com/group/sage-support.
For more options, visit https://groups.google.com/d/optout.

--
You received this message because you are subscribed to the Google Groups "sage-support" group.
To unsubscribe from this group and stop receiving emails from it, send an email to sage-support+unsubscribe@googlegroups.com.

Viviane Pons

unread,
Apr 15, 2015, 2:28:07 PM4/15/15
to sage-s...@googlegroups.com
Ok, just did it. I'm not sure who to CC on this.

Best

Viviane

Vincent Delecroix

unread,
Apr 15, 2015, 2:29:59 PM4/15/15
to sage-s...@googlegroups.com
For the bug report:

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

Vincent
>>>> email to sage-support...@googlegroups.com.
>>>> To post to this group, send email to sage-s...@googlegroups.com.
>>>> Visit this group at http://groups.google.com/group/sage-support.
>>>> For more options, visit https://groups.google.com/d/optout.
>>>>
>>>>
>>>
>> --
>> You received this message because you are subscribed to the Google Groups
>> "sage-support" group.
>> To unsubscribe from this group and stop receiving emails from it, send an
>> email to sage-support...@googlegroups.com.

Volker Braun

unread,
Apr 15, 2015, 2:30:25 PM4/15/15
to sage-s...@googlegroups.com
The segfault is the stack smashed by a ginac infinite recursion, of course.



On Wednesday, April 15, 2015 at 6:55:23 PM UTC+2, Buck Evan wrote:

Vincent Delecroix

unread,
Apr 15, 2015, 2:31:03 PM4/15/15
to sage-s...@googlegroups.com
I would send a message to sage-devel with the ticket number... and you
could also specify it for the readers of this thread.
>>>>> email to sage-support...@googlegroups.com.
>>>>> To post to this group, send email to sage-s...@googlegroups.com.
>>>>> Visit this group at http://groups.google.com/group/sage-support.
>>>>> For more options, visit https://groups.google.com/d/optout.
>>>>>
>>>>>
>>>>
>>> --
>>> You received this message because you are subscribed to the Google Groups
>>> "sage-support" group.
>>> To unsubscribe from this group and stop receiving emails from it, send an
>>> email to sage-support...@googlegroups.com.

Viviane Pons

unread,
Apr 15, 2015, 2:32:11 PM4/15/15
to sage-s...@googlegroups.com
Will do


To post to this group, send email to sage-s...@googlegroups.com.
Visit this group at http://groups.google.com/group/sage-support.
For more options, visit https://groups.google.com/d/optout.


--
You received this message because you are subscribed to the Google Groups
"sage-support" group.
To unsubscribe from this group and stop receiving emails from it, send an

To post to this group, send email to sage-s...@googlegroups.com.
Visit this group at http://groups.google.com/group/sage-support.
For more options, visit https://groups.google.com/d/optout.


--
You received this message because you are subscribed to the Google Groups "sage-support" group.
To unsubscribe from this group and stop receiving emails from it, send an email to sage-support+unsubscribe@googlegroups.com.

Buck Evan

unread,
Apr 15, 2015, 2:37:16 PM4/15/15
to sage-s...@googlegroups.com
It's http://trac.sagemath.org/ticket/18210

We also found that the plot succeeds for 40 and 42 but not 41.
I believe this is because the symbolic form for 40 and 42 are different and simpler than 41:

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(k=41.0), (R, 0, .001))
Launched png viewer for Graphics object consisting of 1 graphics primitive
sage: integrated(k=40)
-1/121645100408832000*gamma(20, 20*R) + 1
sage: integrated(k=41)
1/319830986772877770815625*(319830986772877770815625*sqrt(pi) - 1048576*gamma(41/2, 41/2*R))/sqrt(pi)
sage: integrated(k=42)
-1/2432902008176640000*gamma(21, 21*R) + 1

--
You received this message because you are subscribed to a topic in the Google Groups "sage-support" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/sage-support/7ex6McLE7AA/unsubscribe.
To unsubscribe from this group and all its topics, send an email to sage-support...@googlegroups.com.

Vincent Delecroix

unread,
Apr 15, 2015, 2:44:28 PM4/15/15
to sage-s...@googlegroups.com
41 is odd. What about 43,44,45?
>>>>>>>>> frame #10: 0x0000000103575fb0 libpython2.7.dylib`PyEval_
>>>>>>>>> EvalCodeEx
>>>>>>>>> +
>>>>>>>>> 2144
>>>>>>>>> frame #11: 0x00000001034dae0c libpython2.7.dylib`function_
>>>>>>>>> call
>>>>>>>> email to sage-support...@googlegroups.com.
>>>>>>>> To post to this group, send email to sage-s...@googlegroups.com.
>>>>>>>> Visit this group at http://groups.google.com/group/sage-support.
>>>>>>>> For more options, visit https://groups.google.com/d/optout.
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>> --
>>>>>> You received this message because you are subscribed to the Google
>>>>>> Groups
>>>>>> "sage-support" group.
>>>>>> To unsubscribe from this group and stop receiving emails from it, send
>>>>>> an
>>>>>> email to sage-support...@googlegroups.com.
>>>>>> To post to this group, send email to sage-s...@googlegroups.com.
>>>>>> Visit this group at http://groups.google.com/group/sage-support.
>>>>>> For more options, visit https://groups.google.com/d/optout.
>>>>>>
>>>>>>
>>>>>
>>>>>
>>>>
>>> --
>>> You received this message because you are subscribed to the Google Groups
>>> "sage-support" group.
>>> To unsubscribe from this group and stop receiving emails from it, send an

Buck Evan

unread,
Apr 15, 2015, 2:46:36 PM4/15/15
to sage-s...@googlegroups.com
Similarly, the evens succeed, the odds segfault.
The main difference seems to be the `/sqrt(pi)` as you (Vincent) indicated.


To post to this group, send email to sage-s...@googlegroups.com.
Visit this group at http://groups.google.com/group/sage-support.
For more options, visit https://groups.google.com/d/optout.



  --
You received this message because you are subscribed to the Google
Groups
"sage-support" group.
To unsubscribe from this group and stop receiving emails from it, send
an

To post to this group, send email to sage-s...@googlegroups.com.
Visit this group at http://groups.google.com/group/sage-support.
For more options, visit https://groups.google.com/d/optout.


--
You received this message because you are subscribed to the Google Groups
"sage-support" group.
To unsubscribe from this group and stop receiving emails from it, send an

To post to this group, send email to sage-s...@googlegroups.com.
Visit this group at http://groups.google.com/group/sage-support.
For more options, visit https://groups.google.com/d/optout.


  --
You received this message because you are subscribed to a topic in the
Google Groups "sage-support" group.
To unsubscribe from this topic, visit
https://groups.google.com/d/topic/sage-support/7ex6McLE7AA/unsubscribe.
To unsubscribe from this group and all its topics, send an email to

To post to this group, send email to sage-s...@googlegroups.com.
Visit this group at http://groups.google.com/group/sage-support.
For more options, visit https://groups.google.com/d/optout.


--
You received this message because you are subscribed to a topic in the Google Groups "sage-support" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/sage-support/7ex6McLE7AA/unsubscribe.
To unsubscribe from this group and all its topics, send an email to sage-support+unsubscribe@googlegroups.com.
Reply all
Reply to author
Forward
0 new messages