Archive for the ‘Programming’ Category

Sum of the Euler Totient function

May 10, 2018 Leave a comment

Given a positive integer {n}, the Euler totient function {\varphi(n)} is defined as the number of positive integers less than {n} which are co-prime with {n} (i.e. they have no common factors with {n}). There are formulas for computing {\varphi(n)} starting from the factorization of {n}. One such formula is

\displaystyle \varphi(n) = n \prod_{p|n} \left(1-\frac{1}{p}\right),

where the product is made over all primes dividing {n}.

If you have to compute {\varphi(n)} for all numbers less than a threshold then another property could be useful: {\varphi} is arithmetic, that is, {\varphi(mn) = \varphi(m)\varphi(n)} whenever {\gcd(m,n)=1}. Therefore you could store all values computed until {k} and for computing the value {\varphi(k+1)} there are two possibilities: {k+1=p^\alpha} is a prime power and then {\varphi(k+1) = p^{\alpha}-p^{\alpha-1}} or {k+1} is composite and {k+1 = mn} with {m,n\leq k,\ \gcd(m,n)=1}. Then use the stored values to compute {\varphi(k+1)=\varphi(m)\varphi(n)}.

I now come to the main point of this post: computing the sum of all values of the totient function up to a certain {N}:

\displaystyle \text{Compute } S(N) = \sum_{i=1}^N \varphi(i).

One approach is to compute each {\varphi(i)} 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 {S(10^6)} and about {12} seconds to compute {S(10^7)}. This is super linear in time, since the algorithm computes the factorization for each {n} 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 {S(10^{11})} or even larger. Therefore, there must be more efficient ways to compute {S(N)} out there, so let’s study some of the properties of {S(N)}. In another post I dealt with the acceleration of the computation of the sum of the divisor function.

We have {S(N) = \sum_{i=1}^N \varphi(i)} which is the number of pairs {(a,b)} with {1\leq a\leq b \leq N} such that {\gcd(a,b) =1}. It is not difficult to see that the total number of such pairs is {n(n+1)/2}. Moreover, the possible values of {\gcd(a,b)} are {1,2,...,N}. Now, if for {m \leq N} we search instead for pairs satisfying {\gcd(a,b)=m} then we have {a = ma',\ b = mb'} with {\gcd(a',b')=1} and we get

\displaystyle 1 \leq a' \leq b' \leq N/m,\ \gcd(a',b')=1.

There fore the number of pairs with gcd equal to {m} is {S(\lfloor N/m\rfloor )}. Now we arrive at an interesting recursive formula:

\displaystyle S(N) = \frac{n(n+1)}{2} - \sum_{m=2}^N S(\lfloor N/m \rfloor ).

At a first sight this looks more complicated, but there is a trick to keep in mind whenever you see a summation over {m} of terms of the form {\lfloor N/m \rfloor}: these quantities are constant on large intervals. Indeed,

\displaystyle \lfloor N/m \rfloor = d \Leftrightarrow md \leq N < m(d+1)\Leftrightarrow N/(d+1)<m \leq N/d.

Therefore we can change the index of summation from {m} to {d=\lfloor N/m \rfloor}. The range of {d} for which the interval {I_d = [N/(d+1),N/d]} contains more than one integer is of order {\sqrt(N)}. Indeed, {N/d-N/(d+1) = N/(d(d+1))}. Therefore for {d\leq \sqrt(N)} we should have at least one integer in the interval {I_d}. The part where {d} is larger than {\sqrt{N}} corresponds to {m} smaller than {\sqrt{N}}. Therefore, we can split {S(N)} into two sums, each of order {\sqrt{N}}. and get that

\displaystyle S(N) = \frac{n(n+1)}{2}- \sum_{m=2}^{\sqrt{N}} S(\lfloor N/m \rfloor)-\sum_{d=1}^{\sqrt{N}}\left(\lfloor N/d \rfloor-\lfloor N/(d+1)\rfloor\right)S(d),

where in the last sum we must make sure that {d \neq \lfloor N/d \rfloor} in order to avoid duplicating terms in the sum.

Therefore we replaced a sum until {N} to two sums with upper bound {\sqrt{N}}. The complexity is not {\sqrt{N}}, but something like {N^{2/3}} since we have a recursive computation. Nevertheless, with this new formula and using memoization, to keep track of the values of {S} already computed, we can compute {S(N)} very fast:

{S(10^6) = 303963552392} is computed instantly (vs {1} second with brute force)

{S(10^7) = 30396356427242} takes {1} second (vs {12} seconds with brute force)

{S(10^8)} takes {5} seconds (vs over {3} minutes with brute force)

{S(10^9)} takes {30} seconds

{S(10^{11})} takes about {12} minutes

etc. Recall that these computations are done in Pari GP, which is not too fast. If you use C++ you can compute {S(10^8)} in {0.2} seconds, {S(10^9)} in {1} second and {S(10^{10})} in {6} seconds and {S(10^{11})} 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.

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


    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.

  2. 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.
  3. 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:
    border poly(t=0,1; i){
    Now the mesh could be constructed using a vector of integers NC containing the number of desired points on each of the sides of the polygon:

    mesh Th = buildmesh (poly(NC));
  4. 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.
  5. 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.
  6. 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&lt;N^2;i++){
  xs[i] = floor(i/N)+0.5;
  ys[i] = (i%N)+0.5;

cout &lt;&lt; xs &lt;&lt; endl;
cout &lt;&lt; ys &lt;&lt; endl;

int[int] Nd(N^2);

Nd = 1;
Nd = -nbd*Nd;

cout &lt;&lt; Nd &lt;&lt; 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){


mesh Th = buildmesh(Cd(nb)+Cr(nb)+Cu(nb)+Cl(nb)+disks(Nd));


int[int] bc = [0,1]; // Dirichlet boundary conditions

load &quot;Element_P3&quot;;
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)
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 &lt;&lt; ev[k-1] &lt;<span id="mce_SELREST_start" style="overflow:hidden;line-height:0;"></span>&lt; endl;


Here is the mesh and the solution given by the above program:


Using parfor in Matlab

February 27, 2018 Leave a comment

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:

  1. 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.
  2. 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.
  3. 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

you can do the following:

parfor i=1:N
   array(i) = func(i);

function res = func(i)

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');

If you recieve an error, then run

c = parcluster('local');

and add to path just the right folders for your code to work.

Plotting 3D level sets in Paraview

December 30, 2017 1 comment

Surfaces can be represented as certain level-sets 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 level-set 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 level-sets. 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 level-sets associated to certain values, follow the steps below:

    1. Load your file (containing points with at least a scalar value) and click Apply.
    2. 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.
    3. 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 level-sets you want to see. Click Apply. An example can be seen below.lvlview

If your level-set 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

November 18, 2017 Leave a comment

The text for Problem 613 can be found here. The hint is the following picture 🙂


Metapost – or how to code an image

November 9, 2017 Leave a comment

What software do you use when you need to draw a nice image? I often need to draw things related to research and I never managed to efficiently use a software which uses the mouse or touchpad to draw and modify things. Moreover, if you need to add some mathematical text to the figure things get even more complicated. For me it is more natural to use code to generate graphics, if this is possible. Using something like Matlab to draw is possible. The advantage is that once you have a working code which produces what you want, you can easily modify it. If you have an image where you need to repeat things, using loops can facilitate the job (programmers will understand…).

This is were Metapost comes into play. I found this a long time ago while searching for a tool to build nice graphics for my master thesis. Being used to LaTeX for typesetting math, it was a natural way to draw using code. I do not claim it is the best or the simplest way, but for me it works. Most importantly, it allows me to build high quality, vectorized graphics, which are memory efficient. The advantage of vector graphics is that you can zoom as much as you want and you’ll never see the pixels.

There are a bunch of places where you can learn Metapost. I’ll put here some references I use a lot. The easiest way to start drawing in Metapost is to take a look at some existing codes and modify them. You can find lots of examples in this tutorial. If you have a linux system, you can compile metapost .mp files to obtain high quality pdfs using the command mptopdf. If not, then you can use the online Metapost editor found here. In any case, here are some of my codes for some drawings in Metapost.

Here is the code for my website logo: logo If you want to see the lossless vectorized pdf click to see the following file: logo-0    You can see that it has a border which changes from one color to another. This was done using a loop and a parameter to vary the color as we go along the boundary.

% copy from here to use the online previewer
u:=25; % 25 = 25bp = 25 PostScript points = 25/72 in
wi:=10; % width in units u
he:=7; % height in units u
hoehe:=he*u; % height
%for i=0 upto he:
%draw (0, i*u)--(breite, i*u) withcolor .7white;
%for j=0 upto wi:
%draw (j*u, 0)--(j*u, hoehe) withcolor .7white;
path p,q;
pickup pensquare scaled 20;
draw p;
numeric c,d,detail;
color a,b,co;
for i=1 upto (detail/2):
q := subpath (h*(i-1)/detail,h*i/detail) of p;
draw q withcolor (2*i/detail)[a,b];
for i=detail downto (1+detail/2):
q := subpath (h*(i-1)/detail,h*i/detail) of p;
draw q withcolor (2*(detail-i)/detail)[a,b];
label ("B",(1/15)*(2.6u,2.7u)) scaled 15 withcolor co ;
label ("2" infont defaultfont scaled 8,(4.9u,4.1u))
withcolor co;
% end copy here for online previewer

Here’s another example of a figure where I needed a grid of disks aligned on top of another figure. Using loops in Metapost allowed me to get what I wanted: cm

For the high quality pdf click here: cioranescu-murat-0  and see the code below

u:=25; % 25 = 25bp = 25 PostScript points = 25/72 in
wi:=10; % width in units u
he:=7; % height in units u
hoehe:=he*u; % height

path p,pa;
draw p;
fill p withcolor .8white;
pair a,b;
numeric c,d;
a:= point c of p;
b:= point d of p;
pickup pencircle scaled 1.5;
%draw a;
%draw b;
path q,r,s;
q:= subpath (c,d) of p;
%draw q withcolor red;
%r:= b..((a+b)/2 shifted (.5(a-b) rotated -90))..a;
%draw r;
%fill buildcycle(r,q) withcolor .5red;
%fill pa withcolor .9999blue;
%label.rt(btex $\Omega$ etex,.5(u,3u))  scaled 2;

draw (-u,9u)-- (10u,9u)--(10u,-2u)--(-u,-2u)--cycle;
for i=0 upto 9:
for j=-1 upto 8:
    draw fullcircle scaled 0.3u shifted (i*u, j*u);
    fill fullcircle scaled 0.3u shifted (i*u, j*u) withcolor white;

There are many ways today to draw what you want. Metapost is one of the tools available, if you like coding.

Project Euler – Problem 264

July 28, 2017 Leave a comment

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 x^2+y^2=n, where n is a positive integer, not necessarily a square… It is possible to enumerate all solutions to the following equation for fixed n, simply by looping over x and y. This helps you find all lattice points on the circle of radius \sqrt{n}.

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.

%d bloggers like this: