The original idea of this article is from my primitive trial to solve

Problem 143 of Project Euler. I am not presenting a solution to this problem, instead, I want to demonstrate how to easily launch multiple threads from compiled code without using extra licenses.

Briefly speaking, I am looking for the triangle ABC with interger sides AB = c, BC= a and CA = b , provided that the distances between the

first Fermat Point of this triangle and its vertices are integers (distances are denoted as p, q and r) and the sum p+q+r is preset with value s, as well. Finally I will sum up all distinct values of s under 120,000 which yields at least one integer tuple (a, b, c, p, q, r) .

Here is my code:

Because

Compile[] only supports a limited number of Mathematica function, I need to define some internal function by myself. For instance, IntegerQ[Sqrt[<An Integer>]] will not work in the compiled code as the

IntegerQ cannot be compiled. So this perfect square check function is rewrote in simple functions:

perfectSquareQ =

Compile[{{n, _Integer}},

Module[{sq = Floor[Sqrt[Abs[n]]]}, If[sq^2 == n, True, False]]]

I can quickly check this function as Compile returns an functional object which is very similar to

LinearModelFit In[1]:= perfectSquareQ[5]

Out[1]:= False

Then comes the main body of the procedure. I use "cf" to denote my compiled function. This "cf" basically checks if a given integer sum p+q+r can produce integer sides of a trianle. There may be multiple results or no result. "cf" terminates either when it finds the first solution or after it scans every possibilities in vain (a costly situation).

cf = Compile[{{limit, _Integer}},

Module[{limP = Ceiling[limit/3], r = 1, a = 0, b = 0, c = 0,

res = 0},

Do[

If[res > 0, Break[]];

Do[

r = limit - p - q;

a = Abs[(limit - p)^2 - r*q];

b = Abs[(limit - r)^2 - p*q];

c = Abs[(limit - q)^2 - p*r];

(*If[b>510,Print[{a,b}]];*)

If[perfectSquareQ[a] && perfectSquareQ[b] && perfectSquareQ[c],

res = limit; Break[]]

, {q, Ceiling[(limit - p)/2], limit - p - 1}]

, {p, limit - 2, limP, -1}];

res

],

CompilationTarget -> "C",

CompilationOptions -> {"InlineExternalDefinitions" -> True},

RuntimeAttributes -> {Listable},

Parallelization -> True

]

Besides the function code itself, you can see I have several options in the definition of "cf".

- CompilationTarget: This option tells Mathematica to use the compiled C code instead of runing with "WVM" (wolfram virtual machine)
- CompilationOptions: I use this option because the "cf" calls another compile function "perfectSquareQ"
- RuntimeAttributes/Parallelization: These two options are critical to our parallel computation process. The symbol RuntimeAttributes is similar to SetAttributes, here it let the head "cf" take elements of a list automatically, ie. cf[{1,2,3}] -> {cf[1],cf[2],cf[3]}. While Parallelization let cf[1] , cf[2] and cf[3] get executed on different cpu cores automatically. With these two options, cf[Range[1,100]] can be parallelized with no effort.

As you run this code on your machine (with proper C compiler), you can see check the following statistics from the system monitor:

Of course this method still takes a rather long time to find the possibilites of p+q+r < 30,000 if you do not adopt a smart algorithm with number theory thinges. However, you may save a lot of time to do similar procedure programmings with many jumpy cases, which are not suitable for

Map or

Table function in Mathematica.