Re: LogErfc
- To: mathgroup at smc.vnet.net
- Subject: [mg122650] Re: LogErfc
- From: Peter Falloon <pfalloon at gmail.com>
- Date: Fri, 4 Nov 2011 06:02:36 -0500 (EST)
- Delivered-to: l-mathgroup@mail-archive0.wolfram.com
- References: <j8tkfh$eu0$1@smc.vnet.net>
On Nov 3, 7:47 pm, Shawn Garbett <shawn.garb... at gmail.com> wrote: > I need to have a LogErfc function to be able to compute another > function in a numerically stable manner. I found the function in the > GSL, I found it in R, and even in another system, there's a copy of code to do > it floating around. However, I need it in Mathematica and Log[Erfc[x]] > doesn't cut in due to underflow. Cody solved this in 1969, and I could > code his solution, but I was hoping that someone had already solved > this in Mathematica. Anyone? As you probably know, the problem here is that Erfc[x] goes to zero like Exp[-z^2], and so quickly reaches an underflow, whereas Log[Erfc[x]] is still a number of tractable magnitude. One way around this is to find a representation in terms of another special function in which the exponential dependence is explicit: that way, with a bit of simple manipulation you can apply the Log function and extract the troublesome exponential dependence. One such representation is in terms of the hypergeometric U function (http:// functions.wolfram.com/06.27.26.0002.01). For real z>0 the annoying (z/ Sqrt[z^2]) factor in that equation reduces to 1, and you get: Erfc[z] == Exp[-z^2] * HypergeometricU[1/2, 1/2, z^2] / Sqrt[Pi] Unfortunately, in Mathematica HypergeometricU[1/2,1/2,z^2] is converted to Exp[z^2]*Gamma[1/2, z^2] internally before computing, so we're back where we started (since Gamma[1/2,z^2] suffers the same Underflow[] issues as Erfc for large z). However, if you use the trick of converting one of the 1/2 factors to finite precision (i.e. Real instead of Integer) then a different algorithm is used, which doesn't involve the pesky exponential term. So to cut to the chase, here is a simple definition which should solve your problem: Clear[logErfc] logErfc[z_?NumericQ /; Precision[z]<Infinity] := -z^2 - (1/2)Log[Pi] + Log[HypergeometricU[SetPrecision[1/2, Precision[z]], 1/2, z^2]] Cheers, Peter.