Take the 2-minute tour ×

Consider following simple, illustrating example

cf = Block[{a, x, degree = 3},
  With[{expr = Product[x - a[[i]], {i, degree}]},
   Compile[{{x, _Real, 0}, {a, _Real, 1}}, expr]

This is one of the possible ways to transfer code in the body of a Compile statement. It produces the Part::partd error, since a[[i]] is at the moment of evaluation not a list.

The easy solution is to just ignore this message or turn it off. There are of course other ways around it. For instance one could circumvent the evaluation of a[[i]] by replacing it inside the Compile-body before it is compiled

cf = ReleaseHold[Block[{a, x, degree = 3},
   With[{expr = Product[x - a[i], {i, degree}]},
    Hold[Compile[{{x, _Real, 0}, {a, _Real, 1}}, expr]] /. 
     a[i_] :> a[[i]]]

If the compiled function a large bit of code, the Hold, Release and the replacement at the end goes a bit against my idea of beautiful code. Is there a short and nice solution I have not considered yet?

Answer to the post of Szabolcs

Could you tell me though why you are using With here?

Yes, and it has to do with the reason why I cannot use := here. I use With, to have something like a #define in C, which means a code-replacement at the place I need it. Using := in With delays the evaluation and what the body of Compile sees is not the final piece of code which it is supposed to compile. Therefore,

<< CompiledFunctionTools`
   With[{expr := Product[x - a[[i]], {i, degree}]}, 

shows you, that there is a call to the Mathematica-kernel in the compiled function

I4 = MainEvaluate[ Function[{x, a}, degree][ R0, T(R1)0]]

This is bad because Compile should use only the local variables to calculate the result.


Szabolcs solution works in this case but it leaves the whole expression unevaluated. Let me explain, why it is important that the expression is expanded before it is compiled. I have to admit, my toy-example was not the best. So lets try a better one using With and SetDelayed like in the solution of Szabolcs

Block[{a, x}, With[
  {expr := D[Product[x - a[[i]], {i, 3}], x]}, 

Say I have a polynomial of degree 3 and I need the derivative of it inside the Compile. In the above code I want Mathematica to calculate the derivative for unassigned roots a[[i]] so I can use the formula very often in the compiled code. Looking at the compiled code above one sees, that the D[..] cannot be compiled as nicely as the Product and stays unevaluated

11  R1 = MainEvaluate[ Hold[D][ R5, R0]]

Therefore, my updated question is: Is it possible to evaluate a piece of code without evaluating the Part[]-accesses in it better/nicer than using

Block[{a, x}, With[
  {expr = D[Quiet@Product[x - a[[i]], {i, 3}], x]}, 

Edit: I put the Quiet to the place it belongs. I had it in front of code block to make it visible to everyone that I used Quiet here to suppress the warning. As Ruebenko pointed already out, it should in real code always be as close as possible to where it belongs. With this approach you probably don't miss other important warnings/errors.

Update 2

Since we're moving away from the original question, we should move this discussion maybe to a new thread. I don't know to whom I should give the best answer-award to my question since we discussed Mathematica and Scoping more than how to suppress the a[[i]] issue.

Update 3

To give the final solution: I simply suppress (unfortunately like I did all the time) the a[[i]] warning with Quiet. In a real example below, I have to use Quiet outside the complete Block to suppress the warning.

To inject the required code into the body of Compile I use a pure function and give the code to inline as argument. This is the same approach Michael Trott is using in, e.g. his Numerics book. This is a bit like the where clause in Haskell, where you define stuff you used afterwards.

newtonC = Function[{degree, f, df, colors},
   Compile[{{x0, _Complex, 0}, {a, _Complex, 1}},
    Block[{x = x0, xn = 0.0 + 0.0 I, i = 0, maxiter = 256, 
      eps = 10^(-6.), zeroId = 1, j = 1},
     For[i = 0, i < maxiter, ++i,
      xn = x - f/(df + eps);
      If[Abs[xn - x] < eps,
      x = xn;
     For[j = 1, j <= degree, ++j,
      If[Abs[xn - a[[j]]] < eps*10^2,
        zeroId = j + 1;
     colors[[zeroId]]*(1 - (i/maxiter)^0.3)*1.5
    RuntimeOptions -> "Speed", Parallelization -> True]]@@
    (Quiet@Block[{degree = 3, polynomial, a, x},
     polynomial = HornerForm[Product[x - a[[i]], {i, degree}]];
     {degree, polynomial, HornerForm[D[polynomial, x]], 
      List @@@ (ColorData[52, #] & /@ Range[degree + 1])}])

And this function is now fast enough to calculate the Newton-fractal of a polynomial where the position of the roots is not fixed. Therefore, we can adjust the roots dynamically. Feel free to adjust n. Here it runs up to n=756 fluently

(* ImageSize n*n, Complex plange from -b-I*b to b+I*b *)
With[{n = 256, b = 2.0},
   roots = RandomReal[{-b, b}, {3, 2}],
   raster = Table[x + I y, {y, -b, b, 2 b/n}, {x, -b, b, 2 b/n}]},
       Image[Reverse@newtonC[raster, Complex @@@ roots], "Real"],
       {-b, -b}, {1, 1}, 2 {b, b}]}, PlotRange -> {{-b, b}, {-
         b, b}}, ImageSize -> {n, n}]], {{-b, -b}, {b, b}}, 
   Appearance -> Style["\[Times]", Red, 20]


Dynamic Newton fractal

share|improve this question
Now we have something completely different. Now you are evaluating the Product outside the Compile so you can differentiate it, leaving nothing to do for Compile (in this simple example, of course). Anyway, it looks like to me: 1. first you want to generate code 2. then you want to compile it. In this example you use the result of D to generate code. Let me think for a while. –  Szabolcs Jan 5 '12 at 13:59
But please do check that in your real code there is something left for Compile to do, and it's not merely a plain formula that you pass to it. (In you first example it uses a specialized version of Product so it does do something ...) –  Szabolcs Jan 5 '12 at 14:01
I was wrong here, Compile does speed up simple expressions too, as verified on a long series expansion after turning off auto-compilation by Map. Though vectorization may give better results (if all functions in the expression are Listable, as is the case in a series expansion). –  Szabolcs Jan 5 '12 at 14:09
I'd put the quite as close to the function you want to quiet as you can, Quiet is evil! You will miss real bugs due to it; try to minimize it's effect. –  user1054186 Jan 5 '12 at 14:33
I don't have more time to spend on this today, but I want to point out a caveat with localization! Big warning The fact that With[{y = x}, Compile[{{x, _Real}}, y]] works they way you intend might be an accident. Consider With[{y = x}, Compile[x, y]], outputting CompiledFunction[{x$},x,-CompiledCode-] where x gets localized in Compile. –  Szabolcs Jan 5 '12 at 14:39
show 10 more comments

3 Answers

up vote 2 down vote accepted

One way:

  With[{expr = Quiet[Product[x - a[[i]], {i, degree}]]}, 

be careful though, it you really want this.

share|improve this answer
yes, this is what I meant by ignore this message or turn it off. –  halirutan Jan 5 '12 at 13:07
Here is another way, without Quiet. Undocumented, and actually for a different purpose: cf = Block[{a, x, degree = 3}, With[{expr = D[Product[x - Compile`GetElement[a, i], {i, degree}], x]}, Compile[{{x, _Real, 0}, {a, _Real, 1}}, expr]]]; Compile`GetElement will not check bounds, very careful with this one. –  user1054186 Jan 5 '12 at 14:35
@ruebenko We have a proposal for a mathematica specific site for anything related to mathematica. It would be nice if you could commit (involves creating an Area51 account) to the proposal. We're pretty close to launching (only need about 50 more users). My apologies if you have already done so. –  yoda Jan 5 '12 at 18:56
@yoda, Leonid told me, but I just have not had the time to look at it yet. –  user1054186 Jan 5 '12 at 20:08
@ruebenko, very nice. Who needs boundary checks? They only take valuable cpu-time;-) Thanks. –  halirutan Jan 8 '12 at 2:38
show 2 more comments

EDIT -- Big warning!! Injecting code using With or Function into Compile that uses some of Compile's local variables is not reliable! Consider the following:

In[63]:= With[{y=x},Compile[x,y]]
Out[63]= CompiledFunction[{x$},x,-CompiledCode-]

In[64]:= With[{y=x},Compile[{{x,_Real}},y]]
Out[64]= CompiledFunction[{x},x,-CompiledCode-]

Note the renaming of x to x$ in the first case. I recommend you read about localization here and here. (Yes, this is confusing!) We can guess about why this only happens in the first case and not the second, but my point is that this behaviour might not be intended (call it a bug, dark corner or undefined behaviour), so relying on it is fragile ...

Replace-based solutions, like my withRules function do work though (this was not my intended use for that function, but it fits well here ...)

In[65]:= withRules[{y->x},Compile[x,y]]
Out[65]= CompiledFunction[{x},x,-CompiledCode-]

Original answers

You can use := in With, like so:


It will avoid evaluating expr and the error from Part.

Generally, = and := work as expected in all of With, Module and Block.

Could you tell me though why you are using With here? (I'm sure you have a good reason, I just can't see it from this simplified example.)

Additional answer

Addressing @halirutan's concern about degree not being inlined during compilation

I see this as exactly the same situation as if we had a global variable defined that we use in Compile. Take for example:

In[18]:= global=1
Out[18]= 1

In[19]:= cf2=Compile[{},1+global]
Out[19]= CompiledFunction[{},1+global,-CompiledCode-]

In[20]:= CompilePrint[cf2]
        No argument
        3 Integer registers
        Underflow checking off
        Overflow checking off
        Integer overflow checking on
        RuntimeAttributes -> {}

        I0 = 1
        Result = I2

1   I1 = MainEvaluate[ Function[{}, global][ ]]
2   I2 = I0 + I1
3   Return

This is a common issue. The solution is to tell Compile to inline globals, like so:

     CompilationOptions -> {"InlineExternalDefinitions" -> True}]]];


You can check that now there's no callback to the main evaluator.

Alternatively you can inject the value of degree using an extra layer of With instead of Block. This will make you wish for something like this very much.

Macro expansion in Mathematica

This is somewhat unrelated, but you mention in your post that you use With for macro expansion. Here's my first (possibly buggy) go at implementing macro expansion in Mathematica. This is not well tested, feel free to try to break it and post a comment.

Clear[defineMacro, macros, expandMacros]

macros = Hold[];

SetAttributes[defineMacro, HoldAllComplete]
defineMacro[name_Symbol, value_] := (AppendTo[macros, name]; name := value)

SetAttributes[expandMacros, HoldAllComplete]
expandMacros[expr_] := Unevaluated[expr] //. Join @@ (OwnValues /@ macros)


macros is a (held) list of all symbols to be expanded. defineMacro will make a new macro. expandMacros will expand macro definitions in an expression.

Beware: I didn't implement macro-redefinition, this will not work while expansion is on using $Pre. Also beware of recursive macro definitions and infinite loops.


Do macro expansion on all input by defining $Pre:

$Pre = expandMacros;

Define a to have the value 1:

defineMacro[a, 1]

Set a delayed definition for b:

b := a + 1

Note that the definition of b is not fully evaluated, but a is expanded.




Turn off macro expansion ($Pre can be dangerous if there's a bug in my code):

$Pre =.
share|improve this answer
+1 for := - I had no idea that With can accept the delayed assignment syntax. –  Leonid Shifrin Jan 5 '12 at 12:47
@Leonid That's why I included both Rule and RuleDelayed in withRules :-) Mathematica is just too huge a language, isn't it? No matter how long I've been using it, there's always something seemingly simple that comes up that I didn't know about. Last time it was the TableHeadings option of MatrixForm, which I have painfully reimplemented in the past using Grid (painfully because it drew my mind away from the main task at that time and because not being able to see the vertex ordering in my adjacency matrices was always annoying) –  Szabolcs Jan 5 '12 at 12:54
@halirutan I edited my answer to address this. This is a different issue. The solution is to tell the compiler to inline the definition. –  Szabolcs Jan 5 '12 at 13:14
@Szabolcs, sorry, I gave the wrong reason why := is unfortunately not the solution. I have to examine more carefully. I updated my post. –  halirutan Jan 5 '12 at 13:50
Your macro-expansion translates the problem because I have to use the evaluated form of my expression with defineMacro. –  halirutan Jan 5 '12 at 14:12
show 2 more comments

Ok, here is the very oversimplified version of the code generation framework I am using for various purposes:

SetAttributes[symbolToHideQ, HoldFirst];
symbolToHideQ[s_Symbol, expandedSymbs_] :=! MemberQ[expandedSymbs, Unevaluated[s]];

globalProperties[] := {DownValues, SubValues, UpValues (*,OwnValues*)};

Options[getSymbolsToHide] = {
     Exceptions -> {List, Hold, HoldComplete, 
        HoldForm, HoldPattern, Blank, BlankSequence, BlankNullSequence, 
       Optional, Repeated, Verbatim, Pattern, RuleDelayed, Rule, True, 
       False, Integer, Real, Complex, Alternatives, String, 
       PatternTest,(*Note-  this one is dangerous since it opens a hole 
                    to evaluation leaks. But too good to be ingored *)
       Condition, PatternSequence, Except

getSymbolsToHide[code_Hold, headsToExpand : {___Symbol}, opts : OptionsPattern[]] :=
  Join @@ Complement[
          Flatten[Outer[Compose, globalProperties[], headsToExpand]], code},
            s_Symbol /; symbolToHideQ[s, headsToExpand] :> Hold[s],
            Heads -> True
       Hold /@ OptionValue[Exceptions]];

SetAttributes[makeHidingSymbol, HoldAll];
makeHidingSymbol[s_Symbol] := 
    Unique[hidingSymbol(*Unevaluated[s]*) (*,Attributes[s]*)];

makeHidingRules[symbs : Hold[__Symbol]] :=
     Thread[List @@ Map[HoldPattern, symbs] -> List @@ Map[makeHidingSymbol, symbs]];

reverseHidingRules[rules : {(_Rule | _RuleDelayed) ..}] :=
   rules /. (Rule | RuleDelayed)[Verbatim[HoldPattern][lhs_], rhs_] :> (rhs :> lhs);

FrozenCodeEval[code_Hold, headsToEvaluate : {___Symbol}] :=   
   Module[{symbolsToHide, hidingRules, revHidingRules,  result}, 
      symbolsToHide = getSymbolsToHide[code, headsToEvaluate];
      hidingRules = makeHidingRules[symbolsToHide];
      revHidingRules = reverseHidingRules[hidingRules];
      result = 
         Hold[Evaluate[ReleaseHold[code /. hidingRules]]] /. revHidingRules;
      Apply[Remove, revHidingRules[[All, 1]]];

The code works by temporarily hiding most symbols with some dummy ones, and allow certain symbols evaluate. Here is how this would work here:


  (x-a[[1]]) (x-a[[2]])+(x-a[[1]]) (x-a[[3]])+(x-a[[2]]) (x-a[[3]])]]

So, to use it, you have to wrap your code in Hold and indicate which heads you want to evaluate. What remains here is just to apply ReleseHold to it. Note that the above code just illustrates the ideas, but is still quite limited. The full version of my method involves other steps which make it much more powerful but also more complex.


While the above code is still too limited to accomodate many really interesting cases, here is one additional neat example of what would be rather hard to achieve with the traditional tools of evaluation control:

  Hold[f[x_, y_, z_] := 
    With[Thread[{a, b, c} = Map[Sqrt, {x, y, z}]], 
       a + b + c]], 
  {Thread, Map}]

  f[x_, y_, z_] := 
    With[{a = Sqrt[x], b = Sqrt[y], c = Sqrt[z]}, a + b + c]]
share|improve this answer
I just added globalProperties which I left out when cannibalizing the original code. –  Leonid Shifrin Jan 5 '12 at 17:25
Thanks! Hoping to see a package one day! –  Szabolcs Jan 5 '12 at 17:35
@Szabolcs I hope to see it published one day too! :). Thanks for the upvote. –  Leonid Shifrin Jan 5 '12 at 17:39
This is indeet nice, so +1, but it is of course rather complex. –  halirutan Jan 5 '12 at 20:29
@halirutan Thanks for the upvote. The idea here is, it may be complex, but it is a framework (although tiny), meaning that the tool is rather general. If we were shown the internal implementation of things like With, I am pretty sure that would look complex too. What you requested requires significant deviations from the standard evaluation procedure, and making such things possible in general is a complex task (achieved here only for a special set of use cases) - thus complex code. –  Leonid Shifrin Jan 5 '12 at 20:38
show 1 more comment

Your Answer


