# Issue on numerical_integral with a Pari/GP function

Hello there,

I tried to evaluate the integral of a Pari/GP function but failed so far. Here is an example code on (SageMath version 8.9, Release Date: 2019-09-29 and windows 10):

sage: lchi4 = DirichletGroup(4).list()[1].lfunction(); lchi4
PARI L-function associated to Dirichlet character modulo 4 of conductor 4 mapping 3 |--> -1
sage: numerical_integral(lambda t: lchi4(t).real(), 2,3)
---------------------------------------------------------------------------
SystemError                               Traceback (most recent call last)
/opt/sagemath-8.9/local/lib/python2.7/site-packages/sage/all_cmdline.pyc in <module>()
----> 1 numerical_integral(lambda t: lchi4(t).real(), Integer(2),Integer(3))

/opt/sagemath-8.9/local/lib/python2.7/site-packages/sage/calculus/integration.pyx in
sage.calculus.integration.numerical_integral (build/cythonized/sage/calculus/integration.c:4061)()
353          _b = b
354          W = <gsl_integration_workspace*> gsl_integration_workspace_alloc(n)
--> 355          sig_on()
356          gsl_integration_qag(&F,_a,_b,eps_abs,eps_rel,n,rule,W,&result,&abs_err)
357          sig_off()

SystemError: calling remove_from_pari_stack() inside sig_on()


However, for the Riemann zeta function, numerical_integral works fine as below:

sage: numerical_integral(lambda t: zeta(t), 2,3)
(1.3675256886839795, 1.518258506343328e-14)


It looks like it would work for lchi4 if it can be coerced into a symbolic expression like zeta(x) as shown below:

sage: type(lchi4(3))
<type 'sage.rings.complex_number.ComplexNumber'>
sage: type(zeta(3))
<type 'sage.symbolic.expression.Expression'>


Can you please let me know how to evaluate the integral numerically for a Pari/GP function as lchi4?

Thank you in advance.

edit retag close merge delete

Sort by ยป oldest newest most voted

After searching Google for this issue, I've found myself an answer:

sage: def lchi4(s):
....:     return RR(gp('lfun(lfuncreate(-4), %f)' %s))
....:
sage: lchi4(1)
0.78539816339744830961566084581987572105
sage: numerical_integral(lambda t: lchi4(t), 2,3)
(0.9465779691617932, 1.163268688464501e-09)


I hope this might help for those who have the same problem. Please note that this is not the solution for DirichletGroup(4).list()[1].lfunction() and it only works after changing GP element into RR as shown above.

more