## 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 co-prime 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 *brute-force* approach. For all numerical purposes I will use Pari-GP 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 boundary`C1,C2,C3,C4`

. You wish to build a mesh with the command`mesh 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 of`C3`

, 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 a`plot`

command before trying mesh: something like`plot(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 and`trunc`

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 integers`ind`

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=(1-t)*xx[i]+t*xx[ind(i)];

y=(1-t)*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 with`buildmesh`

. 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 with`problem`

. If you use`solve`

then FreeFem solves the problem where it is defined. If you use`problem`

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^2-1); 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 = N-t;y=N;label=bsquare;} border Cl(t=0,N){x = 0;y=N-t;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=1e-20); // 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[k-1] <<span id="mce_SELREST_start" style="overflow:hidden;line-height:0;"></span>< endl; plot(eV[k-1],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 speed-up 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.

## 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.

## Project Euler 607

If you like solving Project Euler problems you should try Problem number 607. It’s not very hard, as it can be reduced to a small optimization problem. The idea is to find a path which minimizes time, knowing that certain regions correspond to different speeds. A precise statement of the result can be found on the official page. Here’s an image of the path which realizes the shortest time:

## FreeFem to Matlab – fast mesh import

I recently wrote a brief introduction to FreeFem++ in this post. FreeFem is a software designed for the numerical study of partial differential equations. It has the advantage of being able to easily define the geometry of the domain, construct and modify meshes, finite element spaces and solve problems on these meshes.

I use Matlab very often for numerical computations. Most of the numerical stuff I’ve done (take a look here if you want) was based on finite differences methods, fundamental solutions and other classical techniques different from finite elements methods. Once I started using finite elements I quickly realized that Matlab is not that easy to work with if you want some automated quality meshing. PDEtool is good, but defining the geometry is not easy. There is also a simple tool: distmesh which performs a simple mesh construction for simple to state geometries. Nevertheless, once you have to define finite element spaces and solve problems things are not easy…

This brings us to the topic of this post: is it possible to interface Matlab and FreeFem? First, why would someone like to do this? Matlab is easier to code and use than FreeFem (for one who’s not a C expert…), but FreeFem deals better with meshes and solving PDE with finite elements. Since FreeFem can be called using system commands, it is possible to call a static program from Matlab. FreeFem can save meshes and variables to files. Let’s see how can we recover them in Matlab.

There is a tool called “FreeFem to Matlab” developed by Julien Dambrine (link on Mathworks). There’s also a longer explanation in this blog post. I recently tried to use the tool and I quickly found that it is not appropriate for large meshes. It probably scans the mesh file line by line which makes the loading process lengthy for high quality meshes. Fortunately there’s a way to speed up things and I present it below. I will not cover the import of the data (other than meshes) since the function *importdata* from the FreeFem to Matlab tool is fast enough for this.