Archive

Posts Tagged ‘algorithm’

Project Euler tips

A few years back I started working on Project Euler problems mainly because it was fun from a mathematical point of view, but also to improve my programming skills. After solving about 120 problems I seem to have hit a wall, because the numbers involved in some of the problems were just too big for my simple brute-force algorithms.

Recently, I decided to try and see if I can do some more of these problems. I cannot say that I’ve acquired some new techniques between 2012-2016 concerning the mathematics involved in these problems. My research topics are usually quite different and my day to day programming routines are more about constructing new stuff which works fast enough than optimizing actual code. Nevertheless, I have some experience coding in Matlab, and I realized that nested loops are to be avoided. Vectorizing the code can speed up things 100 fold.

So the point of Project Euler tasks is making things go well for large numbers. Normally all problems are tested and should run within a minute on a regular machine. This brings us to choosing the right algorithms, the right simplifications and finding the complexity of the algorithms involved.

ICPC 2015 World Final Problem A

This is the solution to Problem A from the International Collegiate Programming Contest. The list of problems can be found here. This first problem consists simply of reading the parameters of the function defined below, and computing its values on the set ${\{1,2,...,n\}}$. Then, you need to find the maximum decrease in the function values.

The inputs are parameters ${p,a,b,c,d,n}$, the function is

$\displaystyle \text{price}(k) =p(\sin(ak+b)+\cos(ck+d)+2).$

and the values to be considered are ${p(1),p(2),...,p(n)}$. Below is a Matlab code which works, at least for the given Sample Inputs. This should be optimized by removing the for loop. It is instant for $n=1000 fg=000000$ but it should be modified to work until $n=10^6 fg=000000$.

function res = Prob1_2015(p,a,b,c,d,n)
dis = 1:n;

vals = p*(sin(a*dis+b)+cos(c*dis+d)+2);

mat = zeros(n,1);
for i = 1:n
mat(i) = vals(i)-min(vals(i:end));
end

res = max(mat)


The sample outputs are:

 >> Prob1_2015(42,1,23,4,8,10)
104.855110477394

 >> Prob1_2015(100,7,615,998,801,3)
0

 >> Prob1_2015(100,432,406,867,60,1000)
399.303812592112


New version which works for large $n$. It suffices to eliminate the computation of the min at each of the phases of the iteration. This solves the problem in 0.2 seconds for $n= 10^6$

function res = Prob1_2015(p,a,b,c,d,n)
dis = 1:n;
vals = p*(sin(a*dis+b)+cos(c*dis+d)+2);
mat       = zeros(n,1);
lastmax   = vals(1);
for i = 1:n
lastmax   = max(lastmax,vals(i));
mat(i)    = lastmax-vals(i);
end
res = max(mat)

Categories: matlab, Olympiad Tags: , ,

Identifying edges and boundary points – 2D Mesh – Matlab

April 21, 2015 1 comment

A triangulation algorithm often gives as output a list of points, and a list of triangle. Each triangle is represented by the indexes of the points which form it. Sometimes we need extra information on the structure of the triangulation, like the list of edges, or the list of boundary points. I will present below two fast algorithms for doing this.

Finding the list of edges is not hard. The idea is to go through each triangle, and extract all edges. The algorithm proposed below creates the adjacency matrix in the following way: for each triangle ${[i,j,k]}$ we set the elements ${a_{ij},a_{jk},a_{ik}}$ (and their symmetric correspondents) to be equal to ${1}$.

In order to find the points on the boundary (in two dimensions), it is enough to look for the edges which are sides to only one triangle. We can do this using the adjacency matrix. Note that if ${A}$ is the adjacency matrix, then ${A^2=(b_{ik})}$ stores the number of paths of length ${2}$ (two sides) between two points of the triangulation. Note that any edge which is not on the boundary will contain the starting and ending point of two paths of length ${2}$. If ${[i,j,k]}$ is a triangle such that points ${i,j}$ are on the boundary, then ${b_{i,j}=1}$ (there is one path of length ${2}$ going through ${i,k,j}$. We also have ${a_{i,j} = 1}$. Conversely, if ${a_{i,j} = 1}$ and ${b_{i,j} = 1}$ then ${i,j}$ are connected, and there is a unique path of length ${2}$ going from ${i}$ to ${j}$. Thus, ${i,j}$ is an edge on the boundary. Therefore, we just need to identify the indexes ${i}$ such that there exists ${j}$ with ${a_{i,j} b_{i,j} = 1}$.

Below are two short Matlab codes doing these two algorithms. I guess they are close to being optimal, since only sparse and vectorized operations are used.


%p is the list of points
%T is the list of triangles, ap is the number of points
A = min(sparse(T(:,1),T(:,2),1,ap,ap)+sparse(T(:,2),T(:,3),1,ap,ap)+sparse(T(:,3),T(:,1),1,ap,ap),1);
A = min(A+A',1);
% this finds the boundary points, whose indexes are stored in Ibord
B = A^2.*A==1;
Ibord = find(sum(B,2)>0);


Distance from a point to an ellipsoid

Here I will discuss the problem of finding the point which realizes the minimal distance from a given ${x_0}$ to an ellipsoid. First I will look at the theoretical properties which will lead us to ${x_0}$ and then I will present a numerical algorithm to find the nearest point.
Let’s begin with the problem, which can be formulated in a general way like this: we are given ${y \in \Bbb{R}^n}$, and a matrix ${A}$, positive semidefinite, which gives an ellipsoid ${x^TAx=1}$. Without loss of generality assume that ${A}$ is diagonal (if not, we diagonalize it using an orthogonal matrix and work in a new coordinate system). The problem becomes
$\displaystyle \min_{x^TAx=1} \|x-y\|^2.$