# 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().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?

edit retag close merge delete

Sort by » oldest newest most voted

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().lfunction() and it only works after changing GP element into RR as shown above.

more