about Compile and return type

Question

In the mathematica document "tutorial/CompilingWolframLanguageExpressions"

1

It says about how to set return type of a called function you used in Compile.

Clear[com]
com[i_] := Binomial[2 i, i]

test=Compile[{x, {i, _Integer}}, x^com[i], {{com[_], _Integer}}]

From the above example, we can see that because the com function evaluates to integer, so we set {{com[_], _Integer}} to let Compile know the return type of com.

But if you inspect it further

Needs["CompiledFunctionTools`"]
CompilePrint[test]

You can see there is MainEvaluate when calling function com.

So, I don't understand the meaning of this example in the document. If a compiled function has a MainEvaluate process, then I think compiling it is just nonsense, for it won't speed up things, right?

2

Then I came up with another question, it is also mentioned in the same document page. If we compile Sqrt as below

sqrtcom1 = Compile[{{x, _Real}}, Sqrt[x]]

we will run into problems if we evaluate sqrtom1[-1.], it will give errors like

CompiledFunction::cfn: Numerical error encountered at instruction 1; proceeding with uncompiled evaluation. >>

This is because the return type of Sqrt is assumed to be real by default. This can be see from

In[20]:= ToCompiledProcedure[sqrtcom1][[4]]

Out[20]= CompiledResult[Register[Real, 1]]

So, theoretically we could solve this by

sqrtcom2 = Compile[{{x, _Real}}, Sqrt[x], {{Sqrt[_], _Complex}}]

But this is not working!! ToCompiledProcedure[sqrtcom2][[4]] still gives CompiledResult[Register[Real, 1]] and sqrtcom2 still gives errors.

Why is it not working?

Answer 1

Let me counter Daniel Lichtblau's answer

This has zero to do with type inferencing.

by saying, the example in the tutorial you linked is all about type inference. It is not about compiling com to make it faster. It is about helping the compiler to deduce the correct type for the expression.

You have to understand one thing: Highlevel Mathematica language is untyped, which means that it is not known upfront whether Sqrt[x] is an integer, a real, a complex or whether it stays as a general expression. It all depends on the value that x has.

Compiled code is completely different, because all variables will have a type. Either explicitly given by you, or derived/assumed by the compiler.

Therefore, your first question is not about whether or not com inside the compile will be too slow. The question is, can the compiler derive the type of com and therefore assume a correct return type.

Since this example does work correctly even without the explicit type hint, let me give a different example. Here, realf is a function that returns a real number, when the input is an integer. In Mathematica 10.3, this leads to an error message:

realf[i_] := 1.5*i;
f1 = Compile[{{i, _Integer}}, realf[i]]

Calling f1[3] will give you a warning, saying that realf will be the reason that the uncompiled version of f1 is used. If you check the f1 with CompilePrint you will find the line

I1 = MainEvaluate[ Hold[realf][ I0]]

The important part is that I1 means integer register. So because Compile has no information about realf, it assumes this call will be of type integer which is wrong. If we change the definition to

f2 = Compile[{{i, _Integer}}, realf[i], {{realf[_], _Real}}]

and check the compiled code again, we see that now a real register is used for the result of realf.

R0 = MainEvaluate[ Hold[realf][ I0]]

Therfore, f2[3] will run without message since the types are consistent. Nevertheless, realf will be an external call that is evaluated by the kernel.

What Daniel's answer is showing you is, that in the specific example of com being defined as Binomial, you can expand the call and indeed compile all instructions to gain a lot of speed.

As for your second example,

Compile[{{x, _Real}}, Sqrt[x]]

there is one additional thing to note: You call Sqrt[x] where x is the input of type Real. Therefore, the compiler deduces that you want the square-root that works on reals. It seems, not even the type hint at the end of Compile prevents this, but there is a simpler solution:

sqrtHal = Compile[{{x, _Real}}, 
  Module[{xx = 0. I},
   xx = x;
   Sqrt[xx]
  ]
]

Look what we did: we created another variable xx and by giving an initial complex value, we force the type-system to assume xx to be complex. The rest works like in most programming languages. When we assign xx=x a type-conversion takes place and xx is still complex where the imaginary part is zero and the real part is x. Furthermore, the correct complex Sqrt function is selected and therefore

sqrtHal[-1]
(* 0. + 1. I *)

works without problems.

As soon as you have understood the reasons behind this, you easily find another solution:

sqrtHal2 = Compile[{{x, _Real}},
  Sqrt[x + 0. I]
 ]

Answer 2

This has zero to do with type inferencing. There is no Compile support for Binomial. This variant should show the same MainEvaluate invocation:

test2 = Compile[{x, {i, _Integer}}, x^Binomial[2*i, i]]

Let's see if there are evaluator calls:

test2 // InputForm

(* Out[40]//InputForm=
CompiledFunction[{10, 10.3, 5468}, {_Real, _Integer}, 
 {{3, 0, 0}, {2, 0, 0}, {3, 0, 1}}, {{2, {2, 0, 1}}}, {0, 4, 2, 0, 0}, 
 {{15, 1, 0, 2}, {47, Binomial, 2, 0, 2, 2, 0, 0, 2, 0, 3}, 
  {41, 263, 3, 0, 0, 2, 0, 3, 3, 0, 1}, {1}}, 
 Function[{x, i}, x^Binomial[2*i, i]], Evaluate] *)

Same story except slightly more efficient because it is a direct call to Binomial rather than to com. One can evade this, and get a modest speed improvement, by resorting to the FunctionExpand form of Binomial.

test3 = compile[{x, {i, _Integer}}, x^FunctionExpand[com[i]]] /. 
   compile -> Compile;
test3 // InputForm

(* Out[45]//InputForm=
CompiledFunction[{10, 10.3, 5468}, {_Real, _Integer}, 
 {{3, 0, 0}, {2, 0, 0}, {3, 0, 2}}, {{3.141592653589793, {3, 0, 1}}, 
  {2, {2, 0, 1}}, {1, {2, 0, 4}}, {0.5, {3, 0, 3}}}, {0, 6, 7, 0, 0}, 
 {{15, 1, 0, 2}, {41, 263, 2, 0, 1, 2, 0, 2, 2, 0, 3}, 
  {40, 59, 3, 0, 1, 3, 0, 2}, {10, 0, 4}, {13, 3, 4, 5}, 
  {40, 76, 3, 0, 5, 3, 0, 4}, {12, 4, 0, 2}, {40, 76, 2, 0, 2, 2, 0, 5}, 
  {10, 5, 5}, {40, 60, 3, 0, 5, 3, 0, 6}, {10, 3, 5}, {16, 5, 2, 4, 6, 5}, 
  {41, 263, 3, 0, 0, 3, 0, 5, 3, 0, 2}, {1}}, 
 Function[{x, i}, x^((2^(2*i)*Gamma[1/2 + i])/(Sqrt[Pi]*Gamma[1 + i]))], 
 Evaluate] *)

We get modest speed gains as we proceed.

AbsoluteTiming[Do[test[4., 3], {10^6}]]
AbsoluteTiming[Do[test2[4., 3], {10^6}]]
AbsoluteTiming[Do[test3[4., 3], {10^6}]]

(* Out[46]= {1.410492, Null}

Out[47]= {0.904184, Null}

Out[48]= {0.513575, Null} *)

One might get further speedup using Listability and maybe parallelism.

Answer 3

I can't pretend to fully understand this, there are others here like Leonid or Halirutan that talked about this. I notice on this page in the section about "Type Propagation" it says

This shows that for this input, the compiler chooses the result of the square root to also be a real number. This is less general but faster.

So maybe it treats Sqrt specially? This is a workaround,

sqrt[x_] := Sqrt[x];
sqrtcom3 = Compile[{{x, _Real}},
   sqrt[x]
   , {{sqrt[_], _Complex, 0}}];
sqrtcom3[-4.]
(* 0. + 2. I *)