Re: Help evaluating a logarithmic integral

[Bill Allombert <Bill.Allombert@math.u-bordeaux.fr>]

On Thu, Jun 08, 2023 at 12:42:28PM +0200, kevin lucas wrote:
> Thank you both for responding.
> Even after splitting the domain as suggested in the first reply I still get
> a domain error in log with intnum when f is given a floating-point
> argument.
> As for taking the double of the integral over (0,1/2), the integrand is not
> symmetric about 1/2, there's a division by x. Is there something I am
> missing?

I assume you have got this:

? f(v) = intnum(x=0,1/2, log(abs(x^v - (1 - x)^v))/x)+intnum(x=1/2,1, log(abs(x^v - (1 - x)^v))/x);
? f(2)
%2 = -2.4674011002723396547086227499690377838
? f(2.)
  ***   at top-level: f(2.)
  ***                 ^-----
  ***   in function f: intnum(x=0,1/2,log(abs(x^v-(1-x)^v))/x)+intnum
  ***                                 ^-------------------------------
  *** log: domain error in log: argument = 0

The issue is that the expression x^v - (1 - x)^v induces cancellation when x is close to .5.
(so the result get too close to 0 for the log)

You should replace it by a formula that avoid cancellation.

if we set x=1/2+h then we have
(1/2+h)^v - (1/2 - h)^v

the term 1/2^v cancels out, reducing the precision.

when h is close to zero,
exp(v*log((1/2+h)))-exp(v*log((1/2-h)))
is close to 4*v*h

so you can do

t(x,v)=if(abs(x-1/2)>2^-precision(x),x^v - (1 - x)^v,4*v*(x-1/2))
f(v) = intnum(x=0,1/2, log(abs(t(x,v)))/x)+intnum(x=1/2,1, log(abs(t(x,v)))/x)

which gives something:

? f(2)
%35 = -2.4674011002723395953311558693610942877
? f(2.)
%36 = -2.4674011002723395953311558693610942876

I hope this gives the idea.

Cheers,
Bill.