Home > matlab, Olympiad > ICPC 2015 World Final Problem D

ICPC 2015 World Final Problem D


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;
   else
      m = mid;
   end
end
fprintf('Cumulative slice %d : %9.6f mm\n',i,mid/1000);
slices(i) = mid;
end
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);
end

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

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

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

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: