# Compiling a recursive formula

My question is related to computing what is called "invariant measure" for a particular well known fractal - the Sierpinski triangle.

We have an array m of four two by two matrices, say

m = {{{1.22951,-0.102459},{0.058548,1.1856}},{{2.34016,2.08943},{-2.31509,1.50439}},
{{1.52091,2.98751},{2.58555,-2.0641}},{{-11857.1,6214.29},{3214.29,-1678.57}}};


and an array a of four vectors, say:

a = {{-0.127561,-0.238217},{-0.322608,0.247723},{-0.601847,-0.165997},{1304.29,-353.571}};


We also have an array c of four coefficients (related to probabilities and determinants of the matrices), for instance:

c = {0.817674, 1.95388, 2.22763, 180.714};


Then we have four functions of a two-dimensional vector x, for i=1,2,3,4 defined as

w[i_,x_] := m[[i]].x + a[[i]]


The recursive formula I am trying to code is

f[ n_, x_ ] := Sum[c[[i]]*f[ n-1, w[i,x] ],{i,1,4}]


with f[ 0 ,x ] = 1 over the unit square 0 < x[[1]] < 1, 0 < x[[2]] < 1, and zero elsewhere:

f[0, x_] := If[x[[1]] > 0 && x[[1]] < 1 && x[[2]] > 0 && x[[2]] < 1, 1, 0]


I would like to plot functions f[ n, x ] for, say, n=1 up to n=10. It takes a lot of recursions, so I would like to find the optimal code that can be compiled for a fast execution. I tried to compile, but all my efforts failed, probably because I do not really understand what can be compiled and how. Without compiling it takes Mathematica 22 seconds to compute just one value::

Timing[f[10, {0.5, 0.5}]]

{22.667, 0.869726}


Can it be done faster with Mathematica? If so, how? Any help will be appreciated

• You might be interested in this question. – J. M. will be back soon Oct 23 '12 at 7:18
• You can already gain some speed with slight modifications to your definitions : w2[x_] := Dot[m, x] + a ; f2[0, x_] := If[x[[1]] > 0 && x[[1]] < 1 && x[[2]] > 0 && x[[2]] < 1, 1, 0] ; f2[n_, x_] := Dot[c, f2[n - 1, #] & /@ w2[x]]. – b.gates.you.know.what Oct 23 '12 at 7:30
• I do not know. Tested again, with a fresh Kernel, got 10 seconds for my original version with w[i,x], 11 seconds with w2[x][[i]]. Thanks for the link, it deals with recursion indeed, but not with compiling. Without compiling, I am afraid, my problem will not fly on Mathematica's wings, I need to cut the time of execution by a factor of at least 100. – arkajad Oct 23 '12 at 9:17

You can very simply compile if you specify explicitly that fc returns a real; this will get rid of the errors. As pointed out by @asim compilation to "C" does not increase speed in this case.

wc = Compile[{{m, _Real, 3}, {a, _Real, 2}, {x, _Real, 1}},
Dot[m, x] + a (*, CompilationTarget -> "C"*)]

fc = Compile[{{m, _Real, 3}, {a, _Real, 2}, {c, _Real, 1}, {n, _Integer}, {x, _Real, 1}},
If[n == 0, If[x[[1]] > 0 && x[[1]] < 1 && x[[2]] > 0 && x[[2]] < 1, 1, 0],
Dot[c, fc[m, a, c, n - 1, #] & /@ wc[m, a, x]]],
{{wc[_, _, _], _Real, 2},  {fc[_, _, _, _, _], _Real}}(*, CompilationTarget -> "C"*)]

fc[m, a, c, 10, {0.5, 0.5}] // AbsoluteTiming

(* {2.567999, 0.869726} *)

• Thanks. But when I am running: wc = Compile[{{m, _Real, 3}, {a, _Real, 2}, {x, _Real, 1}}, Dot[m, x] + a , CompilationTarget -> "C"], I am getting the message: "A library could not be generated from the compiled function". This is what I was receiving in my original attempts. Some functions would compile to "C", but functions like this one would not. What can be the reason? My CCompilers[] show: {{"Name" -> "Visual Studio", "Compiler" -> CCompilerDriverVisualStudioCompilerVisualStudioCompiler, ... "CompilerName" -> Automatic}} – arkajad Oct 23 '12 at 10:06
• @arkajad Please see edit, I got rid of the errors. – b.gates.you.know.what Oct 23 '12 at 19:34
• I find that compilation alone (not to C) is significantly faster. – asim Oct 23 '12 at 21:35
• @asim Good point, thanks. – b.gates.you.know.what Oct 24 '12 at 7:14
• This is confusing: Calling CompilePrint@fc (after Needs@"CompiledFunctionTools") shows that wc and fc still initiate MainEvaluate calls. Is it alright? Is MainEvaluate` something necessary that cannot be avoided here? – István Zachar Feb 21 '13 at 16:22