Re: Tables of special functions

*To*: mathgroup <mathgroup at yoda.physics.unc.edu>*Subject*: Re: Tables of special functions*From*: David Jacobson <jacobson at cello.hpl.hp.com>*Date*: Wed, 28 Oct 92 21:30:37 -0800

Simon Arridges writes: > I've got a problem with the dereaded Mathematica precision racheting. Using > version 2.0 on Sparc 2 > > I need to produce large (~2000 terms) tables of special functions, > specifically Bessels and Legendres, at quite high precision. ... > [he says that simple evaluation is too slow]. He goes on to say that, following the advice I gave in a recent article in the Mathematica Journal, he tried using SetPrecision in every iteration. (He could have simply used $MinPrecision, except that he is running version 2.0, and $MinPrecision is only available in 2.1.) > Reading the Mathematica journal I tried this : > Version 3) > SetLegendreTables/: SetLegendreTables[ tab_, M_, x_, P_] := > { > tab = Table[0, {i, M+1}]; > tab[[1]] = 1; > tab[[2]] = N[x, P]; > Do[ > {tab[[k+1]] = ((2 k -1) x tab[[k]] - (k -1)tab[[k-1]])/k ; > tab[[k+1]] = SetPrecision[tab[[k+1]],P] }, > {k,2,M} > ] > } > > I.e. keep x symbolic, reset precision to P on each iteration. This seems to > simply zero pad the entries in tab[[k+1]] to the requested precision, so the > accuracy is still disintegrating. I'm confused by his remark that he is keeping x symbolic, but ignoring that, I tried his code, and got perfect results. Now I must admit that I'm running 2.1, not 2.0, and 2.1 has made some significant improvements, but I don't think that it would matter much here. Here is the output I get, evaluating the LegendreP functions at .2, or more precisely (bad pun) at N[2/10,60]. In[1]:= SetLegendreTables/: SetLegendreTables[ tab_, M_, x_, P_] := { tab = Table[0, {i, M+1}]; tab[[1]] = 1; tab[[2]] = N[x, P]; Do[ {tab[[k+1]] = ((2 k -1) x tab[[k]] - (k -1)tab[[k-1]])/k ; tab[[k+1]] = SetPrecision[tab[[k+1]],P] }, {k,2,M} ] } In[2]:= SetLegendreTables[foo,2000,N[2/10,60],60] Out[2]= {Null} In[3]:= Last[foo] Out[3]= 0.01386867441591621201491576703881968766597169195265886984594469 In[4]:= LegendreP[2000,N[2/10,60]] Out[4]= 0.01386867441591621201491576703881968766597169195265886984594 It looks like it works perfectly to me. I suppose I should have asked for more digits in cell 4 just to see if they were all accurate. (I wouldn't expect them to be.) I thought about the possibility that he really passed a symbol in for x, but then the computation remains symbolic and really gets slow, and he already considered that in a piece of the message I chopped off. Now if he really wants it to go fast, he should avoid storing into array elements, especially of long arrays. The above computation took 23.93 seconds on our HP 9000/877 (75 MIPS snakes processor in a box designed for timesharing use). But if we recode it to not use an array, it really zips. Here is code that is strictly functional, if that matters. In[12]:= NextLegendre[x_,k_,lk_,lkm1_] := SetPrecision[((2k+1)x lk -k lkm1)/(k+1),Precision[x]] In[13]:= NextLegendreStep[{x_,k_,lk_,lkm1_}] := {x,k+1,NextLegendre[x,k,lk,lkm1],lk} In[14]:= makeLegendreTable[x_,n_] := Transpose[ NestList[NextLegendreStep,{x,1,x,1},n]][[4]] In[16]:= Timing[baz=makeLegendreTable[N[2/10,60],2000];] Out[16]= {3.64 Second, Null} Not bad for 2001 nearly 60-digit LegendreP functions! In[26]:= Last[baz] Out[26]= 0.01386867441591621201491576703881968766597169195265886984594469 (Just to show its really working.) -- David Jacobson