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