Re: Re: LegendreP error (bug?) in Mathematica

*To*: mathgroup at smc.vnet.net*Subject*: [mg81333] Re: [mg81290] Re: LegendreP error (bug?) in Mathematica*From*: DrMajorBob <drmajorbob at bigfoot.com>*Date*: Wed, 19 Sep 2007 05:36:15 -0400 (EDT)*References*: <fcb3qh$fk1$1@smc.vnet.net> <fcnle7$sac$1@smc.vnet.net> <21737589.1190175665164.JavaMail.root@m35>*Reply-to*: drmajorbob at bigfoot.com

Just curious, mind you... What are these huge values of LegendreP likely to be used for? Bobby On Tue, 18 Sep 2007 04:54:01 -0500, Roman <rschmied at gmail.com> wrote: > Oleksandr, et al., > > Your compiled code works only for small enough L and M, since the > associated Legendre polynomials quickly exceed the range of machine > numbers. It then defaults to slow evaluation. To speed things up you > can start from a stable recursion relation for the spherical harmonics > (like the GSL people do), something like this for Y[L,M,x]= > SphericalHarmonicY[L,M,ArcCos[x],0]: > > Y[L_, M_, x_] := 0/;M>L > Y[L_, M_, x_] := (-1)^M Y[L,-M,x]/;M<0 > Y[L_, M_, x_] := (-1)^L(2L-1)!!Sqrt[(2L+1)/(4Pi(2L)!)](1-x^2)^(L/ > 2)/;M==L > Y[L_, M_, x_] := x Sqrt[2L+1]Y[L-1,L-1,x]/;M==L-1 > Y[L_, M_, x_] := Nest[{#[[1]]+1,#[[3]], > Sqrt[(2#[[1]]+1)/((#[[1]]+M)(#[[1]]-M))](x > Sqrt[2#[[1]]-1]#[[3]]- > Sqrt[(#[[1]]+M-1)(#[[1]]-M-1)/(2#[[1]]-3)]#[[2]])}&, > {M+2,Y[M,M,x],Y[M+1,M,x]},L-M-1][[3]] > > If you compile this like you did, you'll get a good algorithm for the > spherical harmonics which never overflows the machine numbers. From > those you'd use arbitrary precision math to compute > P[L_, M_, x_] := Sqrt[4Pi(L+M)!/((2L+1)(L-M)!)]Y[L,M,x] > which may easily overflow the machine numbers. > > Roman. > > > On Sep 18, 6:50 am, sashap <pav... at gmail.com> wrote: >> On Sep 17, 2:32 am, Roman <rschm... at gmail.com> wrote: >> >> >> >> > John, Oleksandr, Bob, Andrzej, et al., >> >> > If the problem here is arithmetic cancellations, as Oleksandr claims, >> > then Mathematica is definitely not using the best algorithm out there! >> > It is true that when you write out an associated Legendre polynomial >> > explicitely and then plug in numbers, the evaluation of the polynomial >> > becomes terribly unstable because of extreme cancellations. But that's >> > NOT how you're supposed to compute Legendre polynomials; instead, you >> > can use recursion relations, which are numerically stable. >> >> > Andrzej: Thanks for pointing out the issue with MachinePrecision >> > numbers lacking precision tracking. But again, this just explains >> > Mathematica's behavior but not why such a bad algorithm was chosen in >> > the first place. >> >> > I agree that by upping the WorkingPrecision you can make the >> > Mathematica built-in algorithm work (otherwise something would be >> > seriously amiss), but this does not explain why the built-in algorithm >> > is so bad in the first place that it requires extra precision. Here's >> > a more stable algorithm that I've assembled from equations (66,69,70) >> > of MathWorld's Legendre page (http://mathworld.wolfram.com/ >> > LegendrePolynomial.html): (no guarantees on speed) >> >> > P[L_, M_, x_] := 0 /; M > L >> > P[L_, M_, x_] := (-1)^L (2 L - 1)!! (1 - x^2)^(L/2) /; M == L >> > P[L_, M_, x_] := x (2 L - 1) P[L - 1, L - 1, x] /; M == L - 1 >> > P[L_, M_, x_] := Nest[{#[[1]] + 1, #[[3]], >> > (x (2 #[[1]] + 1) #[[3]] - (#[[1]] + M) #[[2]])/(#[[1]] + 1 - M)} >> > &, >> > {M + 1, P[M, M, x], P[M + 1, M, x]}, L - M - 1][[3]] >> >> > It works only for integer L and M, and does not do any checks on the >> > parameters. You can see that this is more stable by plotting something >> > like >> >> > Plot[P[200, 43, x], {x, -1, 1}] >> > Plot[LegendreP[200, 43, x], {x, -1, 1}] >> >> > >From this you can calculate the spherical harmonics using equation >> (6) >> >> > ofhttp://mathworld.wolfram.com/SphericalHarmonic.html >> >> > Y[L_, M_, th_, ph_] := Sqrt[((2 L + 1)/(4 Pi))*((L - M)!/(L + M)!)] * >> > P[L, M, Cos[th]] E^(I*M*ph) >> >> > Roman. >> >> Roman, >> >> Thank you for a very useful observation. Version 6 is using >> hypergeometric function to compute LegendreP, which is >> suboptimal for machine numbers. >> >> Using Compile makes your evaluation scheme quite fast: >> >> lp = Compile[{{n, _Integer}, {m, _Integer}, {z, _Real}}, >> Block[{P, L, M, x}, >> Which[m > n, 0., m == n, (-1)^m (2 m - 1.)!! (1 - z^2)^(m/2), >> m == n - 1, z (2 n - 1) (-1)^m (2 m - 1.)!! (1 - z^2)^(m/2), >> True, Nest[{#1[[1]] + 1, #1[[3]], ( >> z (2 #1[[1]] + 1) #1[[3]] - (#1[[1]] + m) #1[[2]])/(#1[[ >> 1]] + 1 - m)} &, ({m + 1, #1, z (2 m + 1) #1} &)[(-1)^ >> m (2 m - 1.)!! (1 - z^2)^(m/2)], n - m - 1][[3]]]]]; >> >> Clear[P] >> P[n_Integer?NonNegative, m_Integer?NonNegative, >> z_Real /; MachineNumberQ[z] && -1 <= z <= 1] /; >> Developer`MachineIntegerQ[n] && Developer`MachineIntegerQ[m] := >> lp[n, m, z] >> P[n_, m_, z_] := LegendreP[n, m, z] >> >> z is restricted to -1<=z<=1, because LegendreP grows too fast >> outside of it, so Compile will throw an exception. Try >> >> lp[200, 43, -1.2] >> >> With the definition above Plot[P[200,43,x],{x,-1,1}] takes >> half a second and produces reliable result. >> >> Future version of Mathematica will use this method for >> numerical evaluation of LegendreP polynomials. >> >> Regards, >> Oleksandr Pavlyk >> Special Functions Developer >> Wolfram Research >> >> >> >> > On Sep 15, 10:09 am, sashap <pav... at gmail.com> wrote: >> >> > > On Sep 14, 3:13 am, Roman <rschm... at gmail.com> wrote: >> >> > > > I confirm the problem. Just as an example, >> >> > > > In[1] := LegendreP[200, 43, 4/5] // N >> > > > Out[1] = 2.9256424676613492`*^97 >> >> > > > In[2] := LegendreP[200, 43, 0.8] >> > > > Out[2] = 6.151579920980095`*^118 >> >> > > > give strikingly different results! (The former result is >> accurate.) >> >> > > The cause for such a discrepancy is arithmetic cancellations. >> > > You can see it by subtracting two close and sufficiently large >> > > numbers: >> >> > > In[9]:= N[10^17 + 23] - N[10^17] >> >> > > Out[9]= 16. >> >> > > Mathematica's arbitrary precision arithmetic tells you >> > > that the result is not to be trusted: >> >> > > In[10]:= N[10^17 + 23, 16] - N[10^17, 16] >> >> > > Out[10]= 0.*10^1 >> >> > > Something similar happens inside LegendreP algorithm. >> > > The cure is simple, try using arbitrary precision arithmetic. >> >> > > Compare >> >> > > Plot[LegendreP[200, 43, x], {x, 0, 1}] >> >> > > vs. >> >> > > Plot[LegendreP[200, 43, x], {x, 0, 1}, WorkingPrecision -> 75] >> >> > > You can even do >> >> > > Plot[LegendreP[200, 43, x], {x, 0, 1}, >> > > WorkingPrecision -> Infinity] >> >> > > then all arguments x are rationalized (with SetPrecision[x, >> Infinity]) and >> > > LegendreP is evaluated exactly, and then turned into a number from >> > > plotting. >> >> > > Oleksandr Pavlyk >> > > Special Functions developer >> > > Wolfram Research >> >> > > > It seems that this problem occurs only for the associated Legendre >> > > > polynomials with large m; for m=0 the numerical result is >> accurate. >> > > > MathWorld (http://mathworld.wolfram.com/LegendrePolynomial.html) >> gives >> > > > a recursion relation (Eq. 66) for the associated Legendre >> polynomials, >> > > > and I was under the impression that this gave stable results. >> John, >> > > > maybe you can use this recursion relation to get better results , >> or >> > > > you can call the GNU Scientific Library through MathLink >> (http://www.gnu.org/software/gsl/). Bhuvanesh, I am very curious how you >> > > > explain this behavior. >> >> > > > Roman. > > > > -- DrMajorBob at bigfoot.com

**Re: mutual information of red and green channel**

**Identifying Linearly Dependent Equations**

**Re: LegendreP error (bug?) in Mathematica**

**Re: LegendreP error (bug?) in Mathematica**