ICPC 2015 World Final Problem D

May 22, 2015 Leave a comment

This is a solution to problem D from the International Collegiate Programming Contest. The list of the problems can be found here. This is one of the problems which was really appealing for me, from the start. The idea is the following: suppose you have a Swiss cheese cube of {10\times 10\times 10} (the units are centimeters). This cheese cube has multiple holes in it, and for the sake of simplicity, we assume they are spheres. Suppose also that their centers and their radii are known. The problem is how to cut this piece of cheese into horizontal slices having equal weights (that is equal volumes). The holes are supposed to be contained in the cube, and non-intersecting.

The problem is really simple if you don’t have any holes. Just cut at evenly distributed points along the height of the cheese cube. Things change if you add holes, because cutting at evenly distributed points may lead to slices of different weights, since the holes may be not evenly distributed in the cube.

The approach I had in mind was to forget about the slices for a moment. Working with slices means that you may need to cut two caps off a sphere, and that leads to complications. Instead of finding the slices, we could search for the planes dividing the cheese cube into two pieces of given volume ratio. That is, given a proportion of the volume, find the height of the cut, such that that proportion is below the cut. To do this, we only need to be capable to find the volume of the pierced cube which lies under a plane of height {h}.

In order to find this volume, it suffices to subtract from {h\times \text{base area}} the volumes of the portions of the spheres which lie under the plane of height {h}. Thus, we only need the formula for the volume of a spherical cap, which is {V_{cap} = \pi (h^2R - h^3/3 )}, where {R} is the radius of the sphere, and {h} is the height of the cap. Such a function, which computes the volume under a plane of height {h} can be easily vectorized, and implemented in Matlab.

Once we have the function which computes the volume under a certain height, it suffices to do a binary search, in order to find the heights which give volumes {V/s,2V/s,...,(s-1)V/s}. Once we have the heights, we can find the widths of the slices. I implemented the algorithm in Matlab below. I don’t have some complicated test cases, but in the ones I have, the division into {100} slices takes {0.1} seconds.

function Prob3_2015

n = 2; %number of spheres
s = 100; %number of slices

r = [10000 40000]; %radii of the spheres

%coordinates of the centers
%the box is parametrized as 100000^3
%and the real dimensions are 100^3 mm^3

ct = [10000 20000 20000;
      40000 50000 60000];
%total volume of the spheres
totalv = 4*pi/3*sum(r.^3);
%actual volume of the cheese cube
actualv = 100000^3-totalv;
slices = zeros(1,s);

for i = 1:s-1
target  = i*actualv/s;

m = 0;
M = 100000;

while M-m>1e-6
   mid  = (m+M)/2;
   amid = slice_vol(mid,ct,r);
   if amid>target
      M = mid;
      m = mid;
fprintf('Cumulative slice %d : %9.6f mm\n',i,mid/1000);
slices(i) = mid;
slices(end)= 100000;
slices = [0 slices];
disp(' ');
for i = 1:s
  fprintf('Width slice %d : %9.6f mm\n',i,(slices(i+1)-slices(i))/1000);

% function calculating volume under a plane
function res = slice_vol(h,ct,r)

if isempty(ct)
   pars = 0;
   pars = min(max(h-(ct(:,3)-r(:)),0),2*r(:));

res  = h*10^10-sum(pi*pars.^2.*r(:)-pi*pars.^3/3);
Categories: matlab, Olympiad Tags: ,

ICPC 2015 World Final Problem B

May 21, 2015 Leave a comment

This the the solution to Problem B from the International Collegiate Programming Contest. The list of problems can be found here. The idea is to find the maximal area of the intersection of two moving polygons. The inputs give the initial positions of two convex polygons and the vectors giving their speed in the plane.

One idea of a solution is as follows: take a maximal time {T} and a discretization of the interval {[0,T]} by dividing it into {N} parts. Iterate over the translations at times given by this discretization, compute at each step the area of the intersection of the two polygons (if there is any intersection at all), and in the end find the time {t_0} for which this area is maximized. Now, even if the discretization step {T/N} is large (greater than the demanded precision of {10^{-3}}), we can conclude that {t_0} is an approximation of the final time with an error smaller than {T/N}. This is due to the fact that the function representing the area of the intersection has two monotonicity intervals, as a consequence of the fact that the polygons are convex. On the first interval, the intersection area is increasing, on the second one it is decreasing. Thus, once we have a discrete maximum, we are close to the real maximum.

Now, all we are left to do in order to achieve any desired precision is to refine the search near this discrete maximum, find a new, better approximation, and refine again, until we have enough precision. Of course, one first problem with this algorithm is the initial search using a maximal time {T}. It is possible that if {T} or {N} are not large enough, then we do not detect any intersection of the two polygons. Thus, an initial guess, based on a mathematical argument is needed in order to reduce the search of the optimal time to an interval which is small enough to have enough initial precision to detect the discrete maximum.

The algorithm presented below, uses Matlab’s predefined function polybool which can compute the intersection of two polygons. Of course, this is an overkill, since dealing with intersections of convex polygons is not that complicated (but still, I didn’t have enough time to play with the problem, in order to provide a more optimized version). I do not treat the search for the initial time interval. As I think about it, I guess a n argument based on finding some line intersections should give us a narrow enough time interval (with some care for the case when the two speed directions are collinear). The algorithm presented below solves the sample cases, but could fail in more general situations.

function Prob1_2015(p1,p2)

% choose initial polygons; see the text of the problem
p1 = [6 3 2 2 4 3 6 6 6 7 4 6 2 2 2];
p2 = [4 18 5 22 9 26 5 22 1 -2 1];

%p1 = [4 0 0 0 2 2 2 2 0 1 1];
%p2 = [4 10 0 10 2 12 2 12 0 -1 1];

np1 = p1(1);
np2 = p2(1);
cp1 = p1(2:1+2*np1);
cp2 = p2(2:1+2*np2);

vp1 = p1(end-1:end);
vp2 = p2(end-1:end);

cp1 = cp1(:);
cp1 = reshape(cp1,2,np1);

xp1 = cp1(1,:);
yp1 = cp1(2,:);
cp2 = cp2(:);
cp2 = reshape(cp2,2,np2);

xp2 = cp2(1,:);
yp2 = cp2(2,:);

%set precision
prec = 1;

%here there should be a clever initial choice of 
%the starting time and stopping time
%this choice works well in this case
start = 0;
stop  = 5;

while prec>1e-6
n = 100;

timex = linspace(start,stop,n);
areas = zeros(size(timex));

for i = 1:n
   t = timex(i);
   xt1 = xp1+t*vp1(1);
   yt1 = yp1+t*vp1(2);
   xt2 = xp2+t*vp2(1);
   yt2 = yp2+t*vp2(2);
   [x,y] = polybool('intersection',xt1,yt1,xt2,yt2);
   areas(i) = polyarea(x,y);
[m,I] = max(areas);
if m<1e-6
   % if no polygonal intersection is detected
   % there is nothing further to do

tapp = timex(I);
fprintf('Precision %f | Time %f\n',prec,tapp);
start = tapp-prec;
stop  = tapp+prec;
prec  = prec/10;
%this draws the final position of the polygon
%as well as their intersection
axis equal
hold on
axis equal
hold off

Here’s what you get for the first sample test provided in the questions:polygons3

 >> Prob1_2015
    res = 4.193548

ICPC 2015 World Final Problem A

May 21, 2015 Leave a comment

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

res = max(mat)

The sample outputs are:

 >> Prob1_2015(42,1,23,4,8,10)  
 >> Prob1_2015(100,7,615,998,801,3)
 >> Prob1_2015(100,432,406,867,60,1000)

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);
res = max(mat)
Categories: matlab, Olympiad Tags: , ,

Balkan Mathematical Olympiad – 2015 Problems

Problem 1. If {{a, b}} and {c} are positive real numbers, prove that

\displaystyle a ^ 3b ^ 6 + b ^ 3c ^ 6 + c ^ 3a ^ 6 + 3a ^ 3b ^ 3c ^ 3 \ge{ abc \left (a ^ 3b ^ 3 + b ^ 3c ^ 3 + c ^ 3a ^ 3 \right) + a ^ 2b ^ 2c ^ 2 \left (a ^ 3 + b ^ 3 + c ^ 3 \right)}.


Problem 2. Let {\triangle{ABC}} be a scalene triangle with incentre {I} and circumcircle {\omega}. Lines {AI, BI, CI} intersect {\omega} for the second time at points {D, E, F}, respectively. The parallel lines from {I} to the sides {BC, AC, AB} intersect {EF, DF, DE} at points {K, L, M}, respectively. Prove that the points {K, L, M} are collinear. (Cyprus)

Problem 3. A committee of {3366} film critics are voting for the Oscars. Every critic voted just an actor and just one actress. After the voting, it was found that for every positive integer {n \in \left \{1, 2, \ldots, 100 \right \}}, there is some actor or some actress who was voted exactly {n} times. Prove that there are two critics who voted the same actor and the same actress. (Cyprus)

Problem 4. Prove that among {20} consecutive positive integers there is an integer {d} such that for every positive integer {n} the following inequality holds

\displaystyle n \sqrt{d} \left\{n \sqrt {d} \right \} > \dfrac{5}{2}

where by {\left \{x \right \}} denotes the fractional part of the real number {x}. The fractional part of the real number {x} is defined as the difference between the largest integer that is less than or equal to {x} to the actual number {x}. (Serbia)

Source: AoPS

Categories: Olympiad Tags:

Identifying edges and boundary points – 2D Mesh – Matlab

April 21, 2015 Leave a 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
%this computes the adjacency matrix
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);

Simple triangle mesh – Matlab code

April 20, 2015 Leave a comment

There are cases when you need to define a simple, regular mesh on a nice set, like a triangle, or a polygon. There are software which are capable of creating such meshes, but most of them will require some computation time (Delaunay triangulation, plus optimization, edge swapping, etc.). I propose to you a simple algorithm for defining a mesh on a triangle. Of course, this can be extended to polygons or other cases (I used a variant of it to create a mesh of a spherical surface).

The idea is very simple: start with three points {P_1,P_2,P_3}, and define the (only) initial triangle as {T = [1,\ 2,\ 3]}. Then, at each step, take a triangle {t_i} from the list {T}, add the midpoints of {t_i} to the list of points and replace the triangle {t_i} with the four smaller triangles determined by the vertices and midpoints of {t_i}. This is a basic mesh refining strategy. Thus, in short, we start with a triangle and do a number of mesh refinements, until we reach the desired side-length. Of course, at each step we should be careful not to add the midpoint of an edge two times (once for every neighbouring triangle), but in what follows, I will always remove duplicates by bruteforce in the end (not the most economic solution).

Read more…

Categories: Uncategorized Tags: , ,

Vojtech Jarnik Competition 2015 – Problems Category 2

March 27, 2015 Leave a comment

Problem 1. Let {A} and {B} be two {3 \times 3} matrices with real entries. Prove that

\displaystyle A - (A^{-1}+(B^{-1}-A)^{-1})^{-1} = ABA,

provided all the inverses appearing on the left-hand side of the equality exist.

Problem 2. Determine all pairs {(n,m)} of positive integers satisfying the equation

\displaystyle 5^n = 6m^2+1.

Problem 3. Determine the set of real values {x} for which the following series converges, and find its sum:

\displaystyle \sum_{n=1}^\infty \left( \sum_{k_i \geq 0, k_1+2k_2+...+nk_n = n} \frac{(k_1+...+k_n)!}{k_1!...k_n!} x^{k_1+...+k_n}\right).

Problem 4. Find all continuously differentiable functions {f : \Bbb{R} \rightarrow \Bbb{R}}, such that for every {a \geq 0} the following relation holds:

\displaystyle \int_{D(a)} xf\left( \frac{ay}{\sqrt{x^2+y^2}}\right) dxdydz = \frac{\pi a^3}{8}(f(a)+\sin a -1),

where {D(a) = \left\{ (x,y,z) : x^2+y^2+z^2 \leq a^2,\ |y| \leq \frac{x}{\sqrt{3}}\right\}}


Get every new post delivered to your Inbox.

Join 380 other followers

%d bloggers like this: