MathGroup Archive 2005

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: Re: Why is y a local variable here?

  • To: mathgroup at
  • Subject: [mg63245] Re: [mg63217] Re: Why is y a local variable here?
  • From: Daniel Lichtblau <danl at>
  • Date: Tue, 20 Dec 2005 04:19:54 -0500 (EST)
  • References: <> <> <dnv444$oas$> <do0ll8$43k$> <>
  • Sender: owner-wri-mathgroup at

Steven T. Hatton wrote:
> Peter Pein wrote:
>>>SDP[0, x_] := x
>>>SDP[n_Integer?Positive, x_] := Module[{sd, srp = Sqrt[Prime[n]]},
>>>    sd[y_] = SDP[n - 1, y];
>>>    Expand[sd[x + srp] sd[x - srp]]]
>>>SDP[5, x]
>>>Maeder makes in interesting observation at the bottom of page 221 in
>>>_Programming in Mathematica_.
>>Hi Steven,
>>  if you want to save evaluation time with the sd[]-construct, this will
>>  fail,
>>because sd[] will be redefined with evevery call of SDP[]. In that case
>>you should follow the standard procedure:
>>SDP[0,x_]=x; SDP[n_,x_}:=SDP[n,x]=...
> Here's the original code:
> (* :Copyright: 1989-1996 by Roman E. Maeder *)
> (* :Sources: Roman E. Maeder. Programming in Mathematica, 3rd ed.
> Addison-Wesley, 1996.*)
> (* :Discussion: See Sections 4.7 and 5.5 of "Programming in Mathematica"*)
> BeginPackage["ProgrammingInMathematica`SwinnertonDyer`"]
> SwinnertonDyerP::usage = "SwinnertonDyerP[n, var] gives the minimal
> polynomial of the sum of the square roots of the first n primes."
> Begin["`Private`"]
> SwinnertonDyerP[ 0, x_ ] := x
> SwinnertonDyerP[ n_Integer?Positive, x_ ] :=
>     Module[{sd, srp = Sqrt[Prime[n]]},
>         sd[y_] = SwinnertonDyerP[n-1, y];
>         Expand[ sd[x + srp] sd[x - srp] ] ]
> End[]
> EndPackage[]
> Each evaluation of SwinnertonDyerP[ n_Integer?Positive, x_ ] /; n>1 will
> assign the value of SwinnertonDyerP[n-1, y] to its local sd.  Notice that
> Expand[ sd[x + srp] sd[x - srp] ] will only be evaluated after sd[y_] =
> SwinnertonDyerP[n-1, y] has completed.  That means the local sd[] will
> never hold an unevaluated SwinnertonDyerP[n-1, y].  The
> SwinnertonDyerP[n-1, y] will have been replaced by simple expressions
> involving only times and y$ (or x from the top of the stack).
> This means the stack is only constructed once, and all Expand[ sd[x + srp]
> sd[x - srp] ] ] are simple evaluations as the stack is popped.  The sd[]
> are never, as far as I can see, redefined.  If I use a global sd, the
> algorithm works, and the local variables are not created.  The time of
> evaluation is, however, not affected.
> The majority of the time seems to be spent on the Expand[] evaluations.
> Removing it entirely from the algorithm, and simply returning the product
> sd[x + srp] sd[x - srp] increases the execution time considerably.  But
> evaluating Expand[%] on the result takes far longer than doing it within
> the algorithm.

Presumably you meant "decreases the execution time considerably. But..."

In any case, yes, Expand is the culprit, or at least the main one in 
this tragedy. It was not always so. In particular this took a turn for 
the worse re speed around version 5.1. I tracked it to an unintended 
consequence of some improvements in expanding of things with radicals. 
So I expect the speed will return in a future release. Here is what I 
now obtain in the development Mathematica kernel:

SDP[0, x_] := x
SDP[n_Integer?Positive, x_] := Module[{srp = Sqrt[Prime[n]]},
         sd[y_] = SDP[n - 1, y]; Expand[sd[x + srp]*sd[x - srp]]]

In[3]:= Timing[s8 = SDP[8, x];]
Out[3]= {0.364023 Second, Null}

On the same machine, version 5.2 of Mathematica takes around 83 seconds.

A partial workaround is to preexpand anything that might have radicals. 
The variation below will run about the same as the 0.36 seconds above.

SDP2[0, x_] := x
SDP2[n_Integer?Positive, x_] := Module[{srp = Sqrt[Prime[n]]},
      sd2[y_] = SDP2[n - 1, y];
     Expand[Expand[sd2[x + srp]]*Expand[sd2[x - srp]]]]

I apologize for any inconvenience this might have caused.

Daniel Lichtblau
Wolfram Research

  • Prev by Date: Re: Re: Gray's Differential Geometry error?
  • Next by Date: Re: Replacement equivalence?
  • Previous by thread: Re: Why is y a local variable here?
  • Next by thread: Opener Icons for cell group.