Sum of the Euler Totient function
Given a positive integer , the Euler totient function is defined as the number of positive integers less than which are coprime with (i.e. they have no common factors with ). There are formulas for computing starting from the factorization of . One such formula is
where the product is made over all primes dividing .
If you have to compute for all numbers less than a threshold then another property could be useful: is arithmetic, that is, whenever . Therefore you could store all values computed until and for computing the value there are two possibilities: is a prime power and then or is composite and with . Then use the stored values to compute .
I now come to the main point of this post: computing the sum of all values of the totient function up to a certain :
One approach is to compute each and sum them. I will call this the bruteforce approach. For all numerical purposes I will use PariGP in this post. On my computer it takes less than a second to compute and about seconds to compute . This is super linear in time, since the algorithm computes the factorization for each and then sums the values. Using the sieve approach could improve the timing a bit, but the algorithm is still super linear.
In some Project Euler problems it is not uncommon to have to compute something like or even larger. Therefore, there must be more efficient ways to compute out there, so let’s study some of the properties of . In another post I dealt with the acceleration of the computation of the sum of the divisor function.
We have which is the number of pairs with such that . It is not difficult to see that the total number of such pairs is . Moreover, the possible values of are . Now, if for we search instead for pairs satisfying then we have with and we get
There fore the number of pairs with gcd equal to is . Now we arrive at an interesting recursive formula:
At a first sight this looks more complicated, but there is a trick to keep in mind whenever you see a summation over of terms of the form : these quantities are constant on large intervals. Indeed,
Therefore we can change the index of summation from to . The range of for which the interval contains more than one integer is of order . Indeed, . Therefore for we should have at least one integer in the interval . The part where is larger than corresponds to smaller than . Therefore, we can split into two sums, each of order . and get that
where in the last sum we must make sure that in order to avoid duplicating terms in the sum.
Therefore we replaced a sum until to two sums with upper bound . The complexity is not , but something like since we have a recursive computation. Nevertheless, with this new formula and using memoization, to keep track of the values of already computed, we can compute very fast:
is computed instantly (vs second with brute force)
takes second (vs seconds with brute force)
takes seconds (vs over minutes with brute force)
takes seconds
takes about minutes
etc. Recall that these computations are done in Pari GP, which is not too fast. If you use C++ you can compute in seconds, in second and in seconds and in under a minute, if you manage to get past overflow errors.
FreeFem++ Tutorial – Part 2
Click here for the first part. Some other posts related to FreeFem: link 1, link 2.
Here are a few tricks, once you know the basics of FreeFem. If your plan is straightforward: define the domain, build the mesh, define the problem, solve the problem, plot the result… then things are rather easy. If you want to do stuff in a loop, like an optimization problem, things may get more complicated, but FreeFem still has lots of tricks up its sleeves. I’ll go through some of them.
 Defining the geometry of the domain using
border
might be tricky at the beginning. Keep in mind that the domain should be on the left side of the curves definining the boundaries. This could be acheived in the parametrization chosen or you could reverse the parametrization of a particular part of the domain in the following way. Suppose you have the pieces of boundaryC1,C2,C3,C4
. You wish to build a mesh with the commandmesh Th = buildmesh(C1(100)+C2(100)+C3(100)+C4(100));
but you get an error concerning a bad sens on one of the boundaries. If you identify, for example, that you need to change the orientation ofC3
, you can acheive this by changing the sign of the integer defining the number of points on that part of the boundary:mesh Th = buildmesh(C1(100)+C2(100)+C3(100)+C4(100));
You could make sure that you have the right orientations and good connectivities for all the boundaries by running aplot
command before trying mesh: something likeplot(C1(100)+C2(100)+C3(100)+C4(100));
should produce a graph of all your boundaries with arrows showing the orientations. This is good as a debug tool when you don’t know where the error comes from when you define the domain.
 Keep in mind that there are also other ways to build a mesh, like
square
, which meshes a rectangular domain with quite a few options andtrunc
which truncates or modifies a mesh following some criteria. Take a look in the documentation to see all the options.  Let’s say that you need to build a complex boundary, but with many components with similar properties. Two examples come to mind: polygons and domains with many circular holes. In order to do this keep in mind that it is possible to define a some kind of “vectorial boundary”. Let’s say that you want to mesh a polygon and you have the coordinates stored in the arrays
xs,ys
. Furthermore, you have another array of integersind
which point to the index of the next vertex. Then the boundary of the polygon could be defined with the following syntax:
Now the mesh could be constructed using a vector of integers
border poly(t=0,1; i){
x=(1t)*xx[i]+t*xx[ind(i)];
y=(1t)*yy[i]+t*yy[ind(i)];
label=i;
}
NC
containing the number of desired points on each of the sides of the polygon:
mesh Th = buildmesh (poly(NC));
 You can change a 2D mesh using the command
adaptmesh
. There are various options and you’ll need to search in the documentation for a complete list. I use it in order to improve a mesh build withbuildmesh
. An example of command which gives good results in some of my codes is:
Th = adaptmesh(Th,0.02,IsMetric=1,nbvx=30000);
The parameters are related to the size of the triangles, the geometric properties of the griangles and the maximal number of vertices you want in your mesh. Experimenting a bit might give you a better idea of how this command works in practice.  There are two ways of defining the problems in FreeFem. One is with
solve
and the other one is withproblem
. If you usesolve
then FreeFem solves the problem where it is defined. If you useproblem
then FreeFem remembers the problem as a variable and will solve it whenever you call this variable. This is useful when solving the same problem multiple times. Note that you can modify the coefficients of the PDE and FreeFem will build the new problem with the updated coefficients.  In a following post I’ll talk about simplifying your code using macros and functions. What is good to keep in mind is that macros are verbatim code replacements, which are quite useful when dealing with complex formulas in your problem definition or elsewhere. Functions allow you to run a part of the code with various parameters (just like functions in other languages like Matlab).
I’ll finish with a code which computes the eigenvalue of the Laplace operator on a square domain with multiple holes. Try and figure out what the commands and parameters do.
int N=5; //number of holes int k=1; //number of eigenvalue int nb = 100; // parameter for mesh size verbosity = 10; // parameter for the infos FreeFem gives back real delta = 0.1; int nbd = floor(1.0*nb/N*delta*2*pi); int bsquare = 0; int bdisk = 1; // vertices of the squares real[int] xs(N^2),ys(N^2); real[int] initx = 0:(N^21); for(int i=0;i<N^2;i++){ xs[i] = floor(i/N)+0.5; ys[i] = (i%N)+0.5; } cout << xs << endl; cout << ys << endl; int[int] Nd(N^2); Nd = 1; Nd = nbd*Nd; cout << Nd << endl; // sides of the square border Cd(t=0,N){x = t;y=0;label=bsquare;} border Cr(t=0,N){x = N;y=t;label=bsquare;} border Cu(t=0,N){x = Nt;y=N;label=bsquare;} border Cl(t=0,N){x = 0;y=Nt;label=bsquare;} border disks(t=0,2*pi; i){ x=xs[i]+delta*cos(t); y=ys[i]+delta*sin(t); label=bdisk; } plot(Cd(nb)+Cr(nb)+Cu(nb)+Cl(nb)+disks(Nd)); mesh Th = buildmesh(Cd(nb)+Cr(nb)+Cu(nb)+Cl(nb)+disks(Nd)); plot(Th); int[int] bc = [0,1]; // Dirichlet boundary conditions load "Element_P3"; fespace Vh(Th,P2); // variables on the mesh Vh u1,u2; // Define the problem in weak form varf a(u1,u2) = int2d(Th) (dx(u1)*dx(u2) + dy(u1)*dy(u2))+on(1,u1=0)//on(C1,C2,C3,C4,u1=1) +on(bc,u1=0); varf b([u1],[u2]) = int2d(Th)( u1*u2 ) ; // define matrices for the eigenvalue problem matrix A= a(Vh,Vh,solver=Crout,factorize=1); matrix B= b(Vh,Vh,solver=CG,eps=1e20); // we are interested only in the first eigenvalue int eigCount = k; real[int] ev(eigCount); // Holds eigenvalues Vh[int] eV(eigCount); // holds eigenfunctions // Solve Ax=lBx int numEigs = EigenValue(A,B,sym=true,sigma=0,value=ev,vector=eV); cout << ev[k1] <<span id="mce_SELREST_start" style="overflow:hidden;lineheight:0;"></span>< endl; plot(eV[k1],fill=1,nbiso=50,value=1);
Here is the mesh and the solution given by the above program:
Using parfor in Matlab
We all know that loops don’t behave well in Matlab. Whenever it is possible to vectorize the code (i.e. use vectors and matrices to do simultaneous operations, instead of one at a time) significant speedup is possible. However, there are complex tasks which cannot be vectorized and loops cannot be avoided. I recently needed to compute eigenvalues for some 10 million domains. Since the computations are independent, they could be run in parallel. Fortunately Matlab offers a simple way to do this, using parfor.
There are some basic rules one need to respect to use parfor efficiently:
 Don’t use parfor if vectorization is possible. If the task is not vectorizable and computations are independent, then parfor is the way to go.
 Variables used in each computation should not overlap between processors. This is for obvious reasons: if two processors try to change the same variable using different values, the computations will be meaningless in the end.
 You can use an array or cell to store the results given by each processor, with the restriction that processors should work on disjoint parts of the array, so there is no overlap.
The most restrictive requirement is the fact that one cannot use the same variables in the computations for different processors. In order to do this, the simplest way I found was to use a function for the body of the loop. When using a matlab function, all variables are local, so when running the same function in parallel, the variables won’t overlap, since they are local to each function.
So instead of doing something like
parfor i = 1:N commands ... array(i) = result end
you can do the following:
parfor i=1:N array(i) = func(i); end function res = func(i) commands...
This should work very well and no conflict between variables will appear. Make sure to initialize the array before running the parfor, a classical Matlab speedup trick: array = zeros(1,N). Of course, you could have multiple outputs and the output array could be a matrix.
There is another trick to remember if the parpool cannot initialize. It seems that the parallel cluster doesn’t like all the things present in the path sometimes. Before running parfor try the commands
c = parcluster('local'); c.parpool
If you recieve an error, then run
restoredefaultpath c = parcluster('local'); c.parpool
and add to path just the right folders for your code to work.
Plotting 3D level sets in Paraview
Surfaces can be represented as certain levelsets for some 3D functions. Given a set of points with values attached to them a level set associated to a certain number will separate the points into two sets: with values higher or lower than the number considered. It is nice to have good looking plots when working with the levelset method in shape optimization. Paraview is a nice, open source framework, which has the right tools in order to produce high quality plots.
I’ll briefly describe below how to use Paraview to make some nice pictures of levelsets. First of all, you’ll need your data in some format that Paraview can understand. I use vtk file format for which there is a nice automated interface in the software I use for the optimization (FreeFem++). In the vtk file you need to have a set of points and a scalar value attached to them.
If you want to create levelsets associated to certain values, follow the steps below:

 Load your file (containing points with at least a scalar value) and click Apply.
 Next, go to Filters/Alphabetical and select “Cell Data to Point Data” (if you forget to do this, you’ll get a rough surface where you see the discretization). Click Apply.
 Then apply the Contour filter (by clicking the button or going in Filters/Alphabetical). You’ll have to select the field for which the contours will be drawn and then put in the values of the levelsets you want to see. Click Apply. An example can be seen below.
If your levelset cuts the boundary of the domain, Paraview will draw a hole there. If you want to have a closed region instead, you need to use the IsoVolume filter instead of the Contour one. The difference is that you need to specify two values and Paraview will draw the surface enclosing the points corresponding to these values. Many other features are directly available: you can color the level set following another scalar value, you can set the lighting, etc. You can also symmetrize your geometry using the Reflect filter. Below you can see a result built from my work.
You can also create animations in a pretty straightforward way. Just go to View and select the Animation box. Then you’ll see the animation options. Add a Camera object with Orbit field selected. You’ll be presented with multiple options, like the center of rotation, direction of the vertical and initial position. Once everything is set, click the play button to see the animation. Then go to File/Save Animation to save it to a file.
I heard that Paraview could to many things when dealing with visualization aspects, but I hesitated to use it until now since the interface is not straightforward. The use of Filters is not clear in the beginning, but after playing with some examples, everything becomes really easy to use. The next step is to automatize all this using scripts.
Happy New Year!
A hint for Project Euler Pb 613
The text for Problem 613 can be found here. The hint is the following picture 🙂
Project Euler – Problem 264
Today I managed to solve problem 264 from Project Euler. This is my highest rating problem until now: 85%. You can click the link for the full text of the problem. The main idea is to find all triangles ABC with vertices having integer coordinates such that
 the circumcenter O of each of the triangles is the origin
 the orthocenter H (the intersection of the heights) is the point of coordinates (0,5)
 the perimeter is lower than a certain bound
I will not give detailed advice or codes. You can already find a program online for this problem (I won’t tell you where) and it can serve to verify the final code, before going for the final result. Anyway, following the hints below may help you get to a solution.
The initial idea has to do with a geometric relation linking the points A, B, C, O and H. Anyone who did some problems with vectors and triangles should have come across the needed relation at some time. If not, just search for important lines in triangles, especially the line passing through O and H (and another important point).
Once you find this vectorial relation, it is possible to translate it in coordinates. The fact that points A, B, C are on a circle centered in O shows that their coordinates satisfy an equation of the form , where is a positive integer, not necessarily a square… It is possible to enumerate all solutions to the following equation for fixed , simply by looping over and . This helps you find all lattice points on the circle of radius .
Once these lattice points are found one needs to check the orthocenter condition. The relations are pretty simple and in the end we have two conditions to check for the sum of the x and y coordinates. The testing procedure is a triple loop. We initially have a list of points on a circle, from the previous step. We loop over them such that we dont count triangles twice: i from 1 to m, j from i+1 to m, k from j+1 to m, etc. Once a suitable solution is found, we compute the perimeter using the classical distance formula between two points given in coordinates. Once the perimeter is computed we add it to the total.
Since the triple loop has cubic complexity, one could turn it in a double loop. Loop over pairs and construct the third point using the orthocenter condition. Then just check if the point is also on the circle. I didn’t manage to make this double loop without overcounting things, so I use it as a test: use double loops to check every family of points on a given circle. If you find something then use a triple loop to count it properly. It turns out that cases where the triple loop is needed are quite rare.
So now you have the ingredients to check if on a circle of given radius there are triangles with the desired properties. Now we just iterate over the square of the radius. The problem is to find the proper upper bound for this radius in order to get all the triangles with perimeter below the bound. It turns out that a simple observation can get you close to a near optimal bound. Since in the end the radii get really large and the size of the triangles gets really large, the segment OH becomes small, being of fixed length 5. When OH is very small, the triangle is almost equilateral. Just use the upper bound for the radius for an equilateral triangle of perimeter equal to the upper bound of 100000 given in the problem.
Using these ideas you can build a bruteforce algorithm. Plotting the values of the radii which give valid triangles will help you find that you only need to loop over a small part of the radii values. Factoring these values will help you reduce even more the search space. I managed to solve the problem in about 5 hours in Pari GP. This means things could be improved. However, having an algorithm which can give the result in “reasonable” time is fine by me.
I hope this will help you get towards the result.