# Multiprecision numerical evaluation of integral: Sage vs. PARI/GP vs. mpmath

I am trying to compute thousands of integrals of the below type, that comes up in a conformal mapping problem, to as many accurate digits as possible (preferably 50+):

$$\int_{-1}^1\textrm{d}t \frac{\mathcal{Re}\{\log[(\cos{(\pi/130)} - t)]\}}{\sqrt{1 - t^2}}$$

The results from PARI/GP, Sage and Python's mpmath library respectively are:

-2.1770705767584673426214016567105099553,
(-2.1775860588840983, 1.2746272565903925e-05),
-2.1774410877577223893132496923831596284


Clearly $-2.177$ is correct, but what's the best way to find a more accurate answer, accurate to 50+ digits?

I've tried splitting intervals from $[-1, \tau] \cup [\tau, 1]$, where $\tau = \cos{(\pi/130)}$; that doesn't improve things but actually makes it worse. I am working at much higher decimal precision than 3.

UPDATE: Following the suggestion of IgorRivin, the below is the Mathematica attempt:

tau = N[Cos[Pi/130], 50];
epsilon = 2^-10;
limit = N[ArcCos[1 - epsilon], 50];
T1 = NIntegrate[Re[Log[(tau - t)]]/(Sqrt[1 - t^2]), {t, -1, 1 - epsilon},
WorkingPrecision -> 50, AccuracyGoal -> 50]
T2 = NIntegrate[Re[Log[tau - Cos[theta]]], {theta, 0, limit},
WorkingPrecision -> 50, AccuracyGoal -> 50]


gives the results:

Out[1]= -1.7968036050143567231750633164742621583459497767361

During evaluation of In[881]:= NIntegrate::ncvb: NIntegrate failed to converge to prescribed
accuracy after 9 recursive bisections in theta near
{theta} = {0.024170580099064656300274166159163078224443688377727}.
NIntegrate obtained -0.38078136551592029350719560689304843811203900465893
and 6.905015310301168422497720319277708851261355350120750.*^-7 for the
integral and error estimates. >>

Out[2]= -0.38078136551592029350719560689304843811203900465893

Out[3]= -2.1775849705302770166822589233673105964579887813950


As may be seen the error is still in the seventh decimal place; increasing the working precision or accuracy does not really improve things.

UPDATE: The integral is exactly soluble; the result is given in my comment to the accepted answer of Emil Jeřábek, and compares well to the exact value.

• To increase the precision: for mpmath: mpmath.mp.dps=50 for pari \p 50 or default("realprecision",50) – joro Sep 23 '14 at 13:09
• @joro: thanks for the comment, but I am already working at much higher precision than that. e.g. 100 digits for mpmath. Indeed more digits are spewed out for the integral but the error is still in the 4th or 5th decimal; I would like the error to be in 50+th decimal place. – MaviPranav Sep 23 '14 at 13:28
• A, ok. I am not sure any of these can guarantee you bound on the error via numerical methods, though might be wrong. – joro Sep 23 '14 at 13:29
• The integral simplifies to $\int_0^\pi\log\lvert\cos(\pi/130)-\cos\theta\rvert\,d\theta$. Anyway, I'd think your best bet is to first manipulate the integral in some way to get rid of the singularity, which is causing the numerical instability. – Emil Jeřábek Sep 23 '14 at 13:48
• I think there is some serious confusion between WorkingPrecision and AccuracyGoal. If I am not mistaken (I am no Mathematica expert), the former is the machine precision that is used for computations, that is, the local error allowed at each computation. Using WorkingPrecision=50 it is impossible to get 50 significant digits on an ill-conditioned problem. – Federico Poloni Sep 28 '14 at 21:52

Expanding my comment above: putting $\alpha=\pi/130$, the integral equals $$\int_0^\pi\log\left|\cos\alpha-\cos\theta\right|d\theta=\int_0^\pi\log\frac{\cos\alpha-\cos\theta}{(\theta-\alpha)\sin\alpha}\,d\theta+\int_{-\alpha}^{\pi-\alpha}\log\left|\theta\sin\alpha\right|d\theta$$ $${}=\int_0^\pi\log\frac{\cos\alpha-\cos\theta}{(\theta-\alpha)\sin\alpha}\,d\theta+\pi\log\sin\alpha+\alpha(\log\alpha-1)+(\pi-\alpha)(\log(\pi-\alpha)-1).$$ The new integrand is regular on $[0,\pi]$, hence it has a better chance to be accurately approximated by numerical integration. I leave it to someone knowledgeable with such tools to try it.

In case it helps with cancellation errors, one can further write $$\int_0^\pi\log\frac{\cos\alpha-\cos\theta}{(\theta-\alpha)\sin\alpha}\,d\theta=\int_{-\alpha}^{\pi-\alpha}\log\left(\frac{\sin\theta}\theta+\frac{1-\cos\theta}\theta\cot\alpha\right)d\theta.$$

• Great! Using 1000 digit precision mpmath and pari give -2.177586090303602130500688898237613947338583700369286294325795253194308549176741986430328961610663025 ... for the result, with the error from mpmath being 1.0e-717. – MaviPranav Sep 23 '14 at 16:22
• This kind of thing is a great story to tell a student who thinks "the computer can just do it". Having the insight, and the calculus know how, to transform the integrand into something without a singularity, really showcases why we still need "practice integrating". – Steven Gubkin Sep 23 '14 at 16:50
• @StevenGubkin: just to point out that today many of the integration routines are sophisticated enough that very high precision can be obtained for a large class of integrals without further simplification e.g. see pdf of talk by Cohen. Plus I think the need for computing results to 50+ digits might be better appreciated by researchers than by students, especially when most integrals (including this one) can be handled by straightforward software to give 3-4 digits of accuracy. Nevertheless I largely agree with your comment. – MaviPranav Sep 23 '14 at 21:22
• I agree with Steven Gubkin’s comment. Nevertheless, what I’ve done with the integral is not particularly clever, and I can well imagine that a computer algebra system could automate it: in order to compute $\int_a^bf(x)\,dx$, (1) locate the singularities of $f$ in $[a,b]$ and determine their nature, (2) construct a function $g$ with the same singularities whose integral can be easily evaluated (e.g. with explicit antiderivative), (3) compute $\int_a^bf(x)-g(x)\,dx$ numerically. Clearly, current software has room for improvement. – Emil Jeřábek Sep 24 '14 at 11:43
• @EmilJeřábek: The result is actually independent of $\alpha$ and equals $-\pi \log{2}$, as can be obtained from 1.6.10.17, 1.6.10.19 of "Integrals and Series", Vol. 1, by Prudnikov et al.; the numerical result compares very well to the exact result to any desired precision. – MaviPranav Sep 28 '14 at 20:57

Not really appropriate for this site, but Mathematica gives the below (after complaining about convergence problems). On the other hand, integrating your function from $-1$ to $1-\epsilon,$ and using integration by parts (replacing $1/\sqrt{1-x^2}$ by $\arcsin x$) on $(1-\epsilon, 1)$ will work fine.

 NIntegrate[Re[Log[(Cos[Pi/130] - t)]]/(Sqrt[1 - t^2]), {t, -1, 1}, WorkingPrecision -> 1000]

-2.1775780498558904421658718679572230928716272110840621589154768161993\
3732329196521375636873547501076636011363642542650859837723538974595915\
7764788025330214570512607280815043529306101917207642206457155899026608\
2156803691185765082975388097870843163080429578815147518480941010975453\
2564669576068466866559420784664974740875784459505209298049493321604174\
3176415931788039891364286170486010635239758949970468419930675399590464\
0567854619350452364464158872338916765806789170552841942946671912635308\
3221922528590416110327210287192250286801237171954615501203393909462095\
6101020235292043758583706357214383010894691436046713561727357080593836\
8462307036190240989635047153608373293610983565888298857840567497293348\
6324154168093992830510761014856689217520886031941994883181962649076638\
3911888258064741024384159150997182973958068547683198741266508948041775\
2015404055400566473979785323587373843688119801934424131403167705858111\
8040844584952102028151461657122337540572678900659420126432793440968927\
4865463748369420978012

• thanks for the comment and suggestion, but the error is still not in the 1000th digit but at something much less. Moreover I have also tried your substitution suggestion but the error is still in the 7th digit (see update to original question). – MaviPranav Sep 23 '14 at 14:23
• @MaviPranav Well, 7th is better than 5th :) but you can do better yet, I am quite sure... – Igor Rivin Sep 23 '14 at 15:00

To increase the precision: for mpmath: mpmath.mp.dps=50 for pari \p 50 or default("realprecision",50)

Working with precision 1000, mpmath disagrees with Rivin's answer, though the result seems heavily to depend on the precision.

sage: pre=1000;mpmath.mp.dps=pre
sage: time t=mpmath.quad(lambda t:     mpmath.re(mpmath.log(mpmath.cos(mpmath.pi/ mpmath.mpf('130'))-t ))/ mpmath.sqrt(1-t**2),[-1,1])
CPU times: user 34.4 s, sys: 140 ms, total: 34.5 s
Wall time: 34.8 s
sage: t
-2.1775726607115694998616878622345232443775247344010158704666920012976682687035840761365930384169214100410381228168759001755553191218442869332057385292941170091328180080959653384211470618559396443374411645056076930753278504200480523466068297406717463420679346103174787569206292698147279486797487709850966235010193207868896633664816216458559687560722473441911803417313571959410687042439255938037337639274195464203894049531500910228
`