### Archive

Posts Tagged ‘programming’

## Sum of the Euler Totient function

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)

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.

Categories: Number theory, Programming

## FreeFem++ Tutorial – Part 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

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.

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){
x=(1-t)*xx[i]+t*xx[ind(i)];
y=(1-t)*yy[i]+t*yy[ind(i)];
label=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:

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){
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

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

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

## Metapost – or how to code an image

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

prologues:=3;
verbatimtex
%&latex
\documentclass{minimal}
\begin{document}
etex
beginfig(0);
% 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
breite:=wi*u;
%for i=0 upto he:
%draw (0, i*u)--(breite, i*u) withcolor .7white;
%endfor
%for j=0 upto wi:
%draw (j*u, 0)--(j*u, hoehe) withcolor .7white;
%endfor;
path p,q;
p:=(6u,0.5u)--(6u,0)--(0,0)--(0,6u)--(6u,6u)--(6u,5.5u);
pickup pensquare scaled 20;
draw p;
h=length(p);
numeric c,d,detail;
color a,b,co;
a:=(1,0.45,0);
b:=.2white;
co:=.2white;
detail:=500;
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];
endfor;
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];
endfor;
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
endfig;
end;

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:

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

prologues:=3;
verbatimtex
%&latex
\documentclass{minimal}
\begin{document}
etex
beginfig(0);
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
breite:=wi*u;

path p,pa;
p:=(0,4u)..(0,2u)..(2u,2u)..(3u,u)..(3u,0)..(6u,u)..(7u,3u)..(9u,5u)..cycle;
pa:=(3u,5u)..(5u,6u)..(8u,6u)..(8u,3.8u)..(7.5u,3.4u)..(7u,3u)..cycle;
draw p;
fill p withcolor .8white;
pair a,b;
h=length(p);
i=5;
numeric c,d;
c:=0.65*h;
d:=0.7*h;
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;
endfor
endfor;
endfig;
end;

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

Categories: Programming

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

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

October 14, 2016 1 comment

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.