6
$\begingroup$

I recently upgraded to Mathematica 12, and I've been excited to explore new features such as FunctionCompile, which produces llvm code.

I tried compiling the following function:

fc = FunctionCompile[
  Function[Typed[k, "MachineInteger"], Sum[1.0/n, {n, 1, k}]]]

But I got an error message saying "TypeError. Could not find a definition for n." This struck me as odd, because the older Compile function worked fine with the same example, even when compiling to C:

fcc = Compile[{{k, _Integer}}, Sum[1.0/n, {n, 1, k}], 
  CompilationTarget -> "C"]

How do I provide a definition for n to FunctionCompile?

$\endgroup$
  • 5
    $\begingroup$ FunctionCompile is still marked as "experimental" and we have had already other issues with type coercion. Best you tell support about this and wait for one of the upcoming releases. I have to say, I am quite disappointed with FunctionCompile so far; I have not found any useful application in my everyday life that FunctionCompile could hande and Compile could not. And usually, libraries created by FunctionCompile are considerably slower than those created by Compile. But I still have hope that this will improve in upcoming releases. $\endgroup$ – Henrik Schumacher Jul 25 '19 at 17:35
  • 5
    $\begingroup$ Sum is not supported right now. FunctionCompile[Function[Typed[k, "MachineInteger"], Total[Table[1.0/n, {n, 1, k}]]]] should work. $\endgroup$ – ilian Jul 26 '19 at 1:14
5
$\begingroup$

[This is a full rewrite of my original answer]

J. Antonio's question points out how FunctionCompile is indeed still in very early experimental stages of development, that something so simple as computing a finite sum won't compile. But I agree with Henrik that this new compilation engine holds great promise, not just for the greatly expanded number of data types available, but that true functional programming is starting to be supported inside the new CompiledCodeFunction objects that FunctionCompile produces.

That said, I'm going to take a huge leap backwards to suggest that the procedural code shown below will probably be hard to beat in terms of execution time. It works equivalently as a pure function in the Wolfram Engine (albeit slowly), and runs fine in the WVM or C using Compile, and also using the new FunctionCompile mechanism. I'm limited to k=10^7 or so on my little Raspberry Pi 3 Model B owing to its 1GB of RAM and slow ARM processor. But this code should scale quite well with O[n] time and doesn't require building a large table of terms regardless of how large k mioght be.

In[1]:= $Version
Out[1]= "12.0.1 for Linux ARM (32-bit) (June 23, 2019)"

First, the pure function, and its execution time with k=10^7.

In[2]:= f = Function[k,
          Module[{sum = 0.0},
            Do[sum += 1.0/n, {n, k}];
            sum
          ]]

In[3]:= Timing[f[10^7]]
Out[3]= {93.3939, 16.6953}

Not too bad for uncompiled code. But here's where Compile and FunctionCompile will really make it fly.

In[4]:= cf = Compile[{{k, _Integer}},
          Module[{sum = 0.0},
            Do[sum += 1.0/n, {n, k}];
            sum
          ], CompilationTarget -> "C"]

In[5]:= Timing[cf[10^7]]
Out[5]= {0.243868, 16.6953}

In[6]:= ccf = FunctionCompile[
          Function[Typed[k, "MachineInteger"],
            Module[{sum = 0.0},
            Do[sum += 1.0/n, {n, k}];
            sum
          ]]]

In[7]:= Timing[ccf[10^7]]
Out[7]= {1.03342, 16.6953}

In this case, the older Compile engine outperforms FunctionCompile by about 4x. I also experimented with For and While equivalent loops, but they were slower still, and I also looked at Ilian's solution using Total to add up the terms in a Table, but that was slower still and I had to limit myself to k=10^7 as k=10^8 would have required too much memory. The O[n] execution speed of the loop means that I could have pushed the Do version further, though I'd have to wait longer for results with higher values of k, obviously.

You'll note that I was fastidiously declaring the sum value to be local to a Module running within the Function. That's actually quite important for three reasons: 1) global references aren't supported in FunctionCompile, 2) updating the sum value as a global would clobber anything called sum in the current context (usually the global context), and 3) in this simple loop the MathLink overhead of having a CompiledFunction update a global would make the execution time slower, not faster. So it's really important to have sum declared locally to the module for speedy execution.

As an aside, I suspect that Wolfram will eventually make FunctionCompile treat variables and functions operating in the connecting WolframEngine identical to how pure functions and Compile'd functions work. I experimented with using the new KernelFunction as a way to peek and poke variables and functions in the Wolfram Engine (I go into this a bit in the attached notebook), and found it to be a rat's nest and unnecessarily complex when the old way that Compile handled it worked fine (with performance penalties if you abused the MathLink by running looping variables across it, of course).

The reference to n in the iterator, both in my Do-based code, J. Antonio's Sum-based code, and Ilian's Table-based code is worth talking about as well. Strictly speaking from a syntactic perspective, this is a global reference (or a reference to the current context) since there is no declaration of what n is within the Function. However, Sum, Product, ParallelSum, ParallelProduct (none of which are currently supported), Do, Table, and a few other functions make use of iterators which operate by special rules. I won't go into the full details, but if you read Wolfram's documentation on Block you'll see how this has been implemented since the 1990s and before Compile (the WVM-based version) was implemented. It's rather ugly, imho, but suffice to say that you actually are changing a global variable in pure code executing in the Wolfram Engine. When you spot iterators used but not declared in the context of Compile'd or FunctionCompile'd functions, the compilers do indeed properly convert those references into local variables, otherwise there'd be a huge performance penalty in the case of Compile'd functions for having the iterator actually located in the Wolfram Engine and continually accessed over MathLink.

So personally for a long time I've always made a habit of declaring iterators in compiled code as local to a Module, as these two versions below illustrate. They operate at the same speed as the first versions, so declaring the iterator isn't strictly needed, it's just something that gives me a bit of piece of mind, perhaps, knowing that there's a total disconnect from the Wolfram Engine where the iterators are concerned. Plus also, note that the error message J. Antonio received was about the variable n being global because FunctionCompile didn't know Sum was one of the special functions with iterators. I can't help but point out that had n been explicitly declared to be local, then the error message would have been that Sum was not supported.

cf = Compile[{{k, _Integer}},
  Module[{n,sum = 0.0},
    Do[sum += 1.0/n, {n, k}];
      sum
    ], CompilationTarget -> "C"]

ccf = FunctionCompile[
  Function[Typed[k, "MachineInteger"],
    Module[{n,sum = 0.0},
      Do[sum += 1.0/n, {n, k}];
      sum
    ]]]

The final thing I'll point out is that generally speaking I cringe a bit at the overuse of approximate values in Real64 form; I'm an exact math kind of guy. I also noticed that by running the loop in reverse from {n,k,1,-1} instead of {n,1,k} you'll get a bit more accuracy since you won't tend to be adding up numbers with vastly different magnitudes as much. In the k=10^7 case I analyzed, this simple change meant only the two least significant digits were wrong with {n,k,1,-1}, instead of three digits with {n,1,k}. Here is the result of my quick analysis for k=10^7.

In[334]:= ColumnForm[
            {N[maxExactSum, 30],
              NumberForm[maxApproxSum1, Infinity], 
              NumberForm[maxApproxSum2, Infinity]}
          ]
Out[334]= {
            { 16.6953113658598518153991189395},
            {"16.69531136585727"},
            {"16.69531136585996"}
          }

The exact result was computed using this function which cannot be compiled because it relies on arbitrary length integers (another thing I hope Wolfram will eventually support via FunctionCompile). The only difference here from the original is changing the 1.0 to the exact integer 1.

Function[k,Sum[1/n,{n,k}]

That's all for now. I hope this revised answer was a bit more cogent and clear than my first attempt. Happy mathing, all!!

PS -- I just noticed I can't attach files with answers... I'll have to post the notebook some other time somewhere else.

$\endgroup$
  • $\begingroup$ Minor comment. I think you should credit ilian with the Table approach, not Henrik. $\endgroup$ – Carl Woll Nov 1 '19 at 23:59
  • $\begingroup$ Thanks, Carl Woll. I fixed that, and I also incorrectly said For wasn't supported. It turned out my test code was lacking the increment clause of 'n++'. So my use of three arguments as opposed to four didn't match any template FunctionCompile had for For, and I received the same error as if For was not supported. But overall, with new insight into how Compile and FunctionCompile handle global variable references (not just that the latter doesn't), my answer is due for a major overhaul and shortening. I'm likely just to dump this answer and start over again Sunday 11/3/2019. $\endgroup$ – Bert Sierra Nov 2 '19 at 21:00
  • $\begingroup$ Come to think of it, there is also the built-in HarmonicNumber function, which is defined by the identity HarmonicNumber[k] == Sum[1/n,{n,k}]. I don't have my Raspberry Pi with me today or I would run a check to see how fast it evaluates the corresponding rational result. Be sure to use a semicolon so as not to swamp your screen with digits, as in exact=HarmonicNumber[10^7]; See: reference.wolfram.com/language/ref/HarmonicNumber.html $\endgroup$ – Bert Sierra Nov 4 '19 at 16:58
  • 1
    $\begingroup$ Since your comparing machine number results, you can just use HarmonicNumber[1.*^7] which takes less than a millisecond. $\endgroup$ – Carl Woll Nov 4 '19 at 17:16
  • $\begingroup$ Carl -- That's definitely a good way to go for huge values of k; abandon the loops and O[n] execution speed. But I'm still an exact number guy so when I have access to my little RPi tomorrow I'll be checking out HarmonicNumber[10^7]. PS -- Did you mean 10.^7 as your argument? I'm not sure how to read 1.*^7, but maybe that's a way of writing essentially 1.0E7 that I'm not familiar with. $\endgroup$ – Bert Sierra Nov 4 '19 at 17:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.