Description
CSCI 1180 Programming Assignment 1
Problem 1
Euler’s sieve is a souped-up version of the sieve of Eratosthenes, which finds the prime numbers. It works
as follows
L = the list of numbers from 2 to N;
P = 2; /* The first prime */
while (P^{2} < N) {
L1 = the list of all X in L such that P <= X <= N/P.
L2 = P*L1;
delete everything in L2 from L;
P = the next value after P in L;
}
return L;
For instance, for N=27, successive iterations proceed as follows:
Initialization
L = [2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27]
P = 2
First iteration
L1 = [2 3 4 5 6 7 8 9 10 11 12 13]
L2 = [4 6 8 10 12 14 16 18 20 22 24 26]
L = [2 3 5 7 9 11 13 15 17 19 21 23 25 27]
P = 3
Second iteration
L1 = [3 5 7 9]
L2 = [9 15 21 27]
L = [2 3 5 7 11 13 17 19 23 25]
P = 5
Third iteration
L1 = [5]
L2 = [25]
L = [2 3 5 7 11 13 17 19 23]
A. Write a Matlab function EulerSieve1(N) which constructs the Euler sieve, implementing L, L1, L2 as
arrays of integers, as above.
B. Write a Matlab function EulerSieve2(N) which constructs the Euler sieve, implementing L, L1, and
L2 as Boolean arrays, where L[I] = 1 if I is currently in the set L. Thus, the final value returned in the above
example would be the array
[0 1 1 0 1 0 1 0 0 0 1 0 1 0 0 0 1 0 1 0 0 0 1 0 0 0 0]
1
Problem 2:
There is a theorem that states that, if you carry out the following procedure:
P = any polygon (this can be concave or even cross itself).
loop {
compute the midpoint of each side of P
P = the polygon formed by connecting these midpoints in sequence;
}
Then P will converge toward a series of points that lie on an ellipse. Picture on the next page.
A. Assume that the polygon is represented by two vectors: Px, the array of x-coordinates and Py, the array
of y-coordinates. For example, the polygon with vertices at h0, 0i, h2, 8i, h4, 0i, – h−2, 6i, h6, 6i would be
represented as the arrays
Px = [0, 2, 4, -2, 6]
Py = [0, 8, 0, 6, 6]
Write a Matlab function [Qx,Qy]=ConnectMidpoints(P) that, given a polygon Px,Py constructs the coordinates of the midpoints of P in sequence. For instance if P is the matrix above then ConnectMidpoints(P)
would return
Qx = [1, 3, 1, 2, 3]
Qy = [4, 4, 3, 6, 3]
Each element of Qx is the average of two consecutive columns of Py, and likewise with Qy. For instance
Qx(1) = (Px(1)+Px(2))/2.
The last element of Qx is the average of the last and first element of Px e.g. in the above example Qx(5) =
(Px(5)+Px(1))/2.
Your code should of course work for polygons with any number of points, not just polygons with 5 points.
B. Write a Matlab function ConvergingPolygons(P,N) which takes as input a polygon P and a number
N and draws pictures of the first N polygons in this sequence, starting with P. As in the pictures on the
next page, it should alternate pictures that show the polygon with pictures that show the polygon and the
midpoints.
Let Matlab adjust the scale on each successive picture, or the picture will soon become too small to see.
Also, as always with geometric drawings in Matlab, call axis equal to make sure that the x and y axes
have the same scale.
Note: To make multiple plots on a single figure, use hold on and hold off. To make multiple figure, use
figure(). So the code inside the loop that generated the even-numbered pictures on the next page had the
form
In calling the plot function, be sure to repeat the first point at the so that the polygon closes.
2
Output of ConvergingPolygons(P,6) with x- and y- coordinate initialized by the following Matlab calls:
Px = 7+[1,3,5,7,2.5,4,6].*cos(10*pi*(1:7)/7)
Py = 7+[1,3,5,7,2.5,4,6].*sin(10*pi*(1:7)/7)
3
After 9 more iterations, it straightens itself out and then continues to get closer to a very long, thin ellipse.
4
CSCI 1180 Programming Assignment 2
Problem 1 (5 points)
In programming assignment 1, problem 2, you wrote a program that took a polygon P as input
and found a polygon Q that was the midpoints of all the sides of P. In this assignment you will
invert that process; that is, you will write a program that takes the midpoints Q and either returns
a polygon P such that Q is the midpoints of the sides of P, or reports that no such polygon exists.
That is, given a sequence of k points in the plane ~q1 . . . ~qk, find a sequence ~p1 . . . ~pk such that
~q1 is the midpoint of ~p1 and ~p2.
~q2 is the midpoint of ~p2 and ~p3.
. . .
~qk−1 is the midpoint of ~pk−1 and ~vk.
~qk is the midpoint of ~pk and ~p1.
Write a Matlab function InvertMidpoints(Q) which takes as argument a polygon Q, represented
in the form used in programming assignment 1, and returns a pair of values [FOUND, P]. If there is
a solution, then FOUND is true and P is a 2 × k matrix where P[:,i] is the x- and y- coordinates of
the ith vertex of P. If there is no solution then FOUND is false and P is the 2 × k zero matrix.
Thus, for example, if Q is the array
Q =
1 3 1 2 3
4 4 3 6 3
then InvertMidpoint(Q) should return FOUND=true and
P =
0 2 4 −2 6
0 8 0 6 6
If
Q =
1 3 1 −1
4 4 3 3
then InvertMidpoint(Q) should return FOUND=true and
P =
0 2 4 −2
0 8 0 6
1
If
Q =
1 3 1 −1
4 4 3 6
then InvertMidpoint(Q) should return FOUND=false and
P =
0 0 0 0
0 0 0 0
First hint: First, notice that the x and y coordinates in this problem are completely decoupled.
That is, if we write ~qi = hxi
, yii, and ~pi = hai
, bii, then we have
x1 = (a1 + a2)/2.
x2 = (a2 + a3)/2.
. . .
xk−1 = (ak−1 + ak)/2.
xk = (ak + a1)/2.
and the identical equations relates the y’s to the b’s. Thus, for example, for k = 5, the 5-dimensional
vectors ~x, ~y, ~a,~b satisfy the equations ~x = C~a, ~y = C~b where C is the matrix of coefficients,
C =
1/2 1/2 0 0 0
0 1/2 1/2 0 0
0 0 1/2 1/2 0
0 0 0 1/2 1/2
1/2 0 0 0 1/2
So if there exists a solution you can find it by solving this system of equations separately in for the
x- and the y-coordinates and then putting those solutions together.
(Do NOT write code to solve systems of linear equations; use the built in Matlab solver.)
Second hint: As you will see in problem set 3, the matrix C is always of rank k if k is odd and of
rank k − 1 if k is even. So if k is odd then you can just construct the above matrix and solve the
problem. If k is even, then Matlab won’t allow you to do that, even if there is a solution. Rather
you have to (a) throw out one of the equations; (b) solve the remaining set of k − 1 equations in k
unknowns; (c) check to see whether the solution you have found satisfies the last equation (i.e. to
within some error tolerance ǫ which you can take to be 10−8
.)
Problem 2
Note: This problem is long, but the code is actually very short; perhaps 20 lines in total.
Suppose that A and B are electrically charged objects, located at points pA and pB with charges
QA and QB. Then the force F~A(B) that B exerts on A is the vector
F~A(B) = QA · QB
|pA − pB|
2
·
pA − pB
|pA − pB|
2
In the above product, the first factor is the magnitude of the force, which is the product of the
charges divided by the distance squared; the second factor is the direction of the force, which is the
direction from B to A.
If there are several objects B1 . . . Bk exerting a force on A, then the total force on A is the sum of
the forces:
F~A({B1 . . . Bk}) = X
k
i=1
F~A(Bi)
If the charge on A and the position of all the charges is fixed, then the net force is a linear function
of vector of charges hQ~ = Q1 . . . Qki.
For instance, in two dimensions, we could have the following situation, illustrated in the picture.
Object Location Charge |pA − pB| Magnitude of F~A(B) F~A(B)
A h0, 1i 1 — —
B1 h4, 4i 50 5 50/25 = 2 2 · h−4, −3i/5 = h−1.60, −1.20i
B2 h−3, 1i 36 3 36/9 = 4 4 · h3, 0i/3 = h4.00, 0.00i
B3 h1, 0i −6
√
2 −6/2 = −3 −3 · h−1, 1i/
√
2 = h2.12, −2.12i
Total h4.52, −3.32i
A
B1
B2
B3
Total
Force from B2
Force from B3
force from B1
3
Problem 2.A (2.5 points)
Write a function function F = ForceMatrix(PA,PB) where
• PA is a 2-dimensional column vector of the coordinates of object A of charge 1.
• PB is a 2 × k matrix, where the ith column, PB[:,i] is the coordinates of object Bi
.
• F, the value returned, is the 2 × k matrix with the property that for any vector of charges ~Q,
the value F · ~Q is the net force on A.
For instance, in the above example, we could call
> PA = [0;1];
> PB = [4,1,-3; 4,0,1];
> F = ForceMatrix(PA,PB)
F =
-0.0320 -0.3536 0.1111
-0.0240 0.3536 0
> QB = [50; -6; 36];
> F*QB
ans =
4.5213
-3.3213
Problem 2.B: 0.5 points
Write the following two functions: function TF = TotalForce(PA,PB,QB) and C = PossibleCharge(PA,PB,TF).
In both of these PA, PB are the same as in problem 1. In TotalForce, the input QB is a column
vector of the charges on B and the value returned TF is the total force on A, a column vector. In
PossibleCharge, TF is the total force as a column vector and the value returned C is a possible
charge vector that would give rise to that force. If k > 2 then there are multiple possible answers
but your code only has to return one of these. For example, using the same values of PA,PB,QB we
could write,
> TF = TotalForce(PA,PB,QB)
TF =
4.5213
-3.3213
> C = PossibleCharge(PA,PB,TF)
0
-9.3941
10.8000
>> TotalForce(PA,PB,C)
ans =
4.5213
-3.3213
4
Having done problem 2.A, each of these functions should consist of one quite simple line of MATLAB.
The code for TotalForce should always work, unless A is at the same position as one of the Bi
’s.
If all the points, are collinear, then there may not exist any solution for PossibleCharge.
Problem 2.C: 2 points
Suppose as before there are k fixed charges B1 . . . Bk in the plane. You know the locations, but not
the value of the charges, and you want to find out the value of the charges. A way to do this is
as follows: You take an object A with charge 1, you put it at various points in the plane, and you
measure the net force on it.
Write a function function C = FindCharges(PA,PB,TF) where
• PA is a 2 × w matrix, where the ith column, PA[:,i] is the coordinates of the ith placement
of the test charge A. The dimension w is the number of different placements you try.
• PB is the locations of the charges B1 . . . Bk, as above.
• F is a 2 × q matrix, where the ith column F[:,i] is the total force on A in its ith placement.
• The value returned C is the k-dimensional column vector of charges on the Bi
.
Hint: Look up the Matlab reshape function.
For instance, in the above example, we could call
> PA = [0,2;1,0];
> PB = [4,1,-3; 4,0,1];
> TF(:,1) = TotalForce(PA(:,1),PB,QB);
> TF(:,2) = TotalForce(PA(:,2),PB,QB);
> C = FindCharges(PA,PB,TF)
C =
50.0000
-6.0000
36.0000
5
CSCI 1180 Programming Assignment 3
In this assignment, you will use Matlab to show how the appearance of a 3D object changes as a
result of rotation (assuming that the distance to the object is constant, and large as compared to
the size of the object).
We will greatly simplify the problem by limiting it to convex polyhedra. The advantage of using
convex polyhedra is that it greatly simplifies the problem of determining what parts of the surface
are visible; a face of a convex polyhedron is visible if the normal to the face points toward the viewer.
The assigment is to write a function DrawRotatedPolyhedron(M,P). The input parameters are M, a
3×3 rotation (orthogonal) matrix, and P, a data structure representing a convex polyhedron. What
the function does is to draw a two-dimensional picture of P after rotation by M in a form to be
described below.
The input data structure representing of an n-face polyhedron P is a cellular array of size n, where
each cell is a face of the polyhedron. A face of the polyhedron with k vertices is a 3 × k array where
each column of the array contains the coordinates of the vertices. The columns are in counterclockwise order, as viewed by someone outside the solid looking at the face.
For example, the unit cube is represented as a cellular array { X1,X2,Y1,Y2,Z1 Z2 } where X1 is the
face at the low end of the x-axis, X2 is the face at the high end of the x-axis, and so on. Specifically,
X1 =
0 0 0 0
0 0 1 1
0 1 1 0
X2 =
1 1 1 1
0 1 1 0
0 0 1 1
Y 1 =
0 1 1 0
0 0 0 0
0 0 1 1
Y 2 =
0 0 1 1
1 1 1 1
0 1 1 0
Z1 =
0 0 1 1
0 1 1 0
0 0 0 0
Z2 =
0 1 1 0
0 0 1 1
1 1 1 1
This and other shapes are defined as scripts in file polyhedra.m. The function EulerRotation(Psi,Phi,Theta)
generates the 3D rotation with Euler angles Psi, Phi, Theta. These can be found in the course code
library, linked from the course web site.
The picture of the shape will reflect the appearance of the rotated shape as seen from below; thus,
the projection of the faces that are visible below onto the x-y plane (see figure). The picture will
show (1) the visible vertices and edges (2) in the center of each face, the projection of a small circle
around the center (which will project as an ellipse).
Constructing this picture involves the following steps:
• Apply the matrix M to each of the points in P, generating a new polyhedron P
′
. (Since the
absolute position of of the polyhedron does not matter, you can use natural coordinates and
a 3 × 3 rotation matrix rather than homogeneous coordinates and a 4 × 4 matrix.)
• Exclude faces that are not visible from below. To do this, compute the outward normal to
each face, and exclude any face whose normal has a positive z-component. To find the outward
normal to a face, choose any three consecutive vertices of the face a, b, c; the cross-product
(c − b) × (a − b) is the outward normal.
• Project the vertices and edges onto the x−y plane simply by ignoring the z coordinate. Draw
these in the picture.
1
• Compute the central circles on each visible face as follows:
– Compute the center o of the face as the average of the vertices of the face.
– Choose a radius r as half the distance from o to the nearest of the edges (use the formula
for the distance from a point to an edge discussed in class)
– Find two orthogonal unit-length vectors ˆu, vˆ in the plane of the face. (If a, b, c are vertices
in the face, then you can choose dir(b − a) as one and the perpendicular from c to line
ab as the other.)
– Compute points on the circle in 3D as o + r cos(2πt/N)ˆu + r sin(2πt/N)ˆv for t = 0 . . . N.
– Project these points onto the x − y plane by ignoring the z-coordinate, and connect the
points to plot the ellipse.
CSCI 1180 Programming Assignment 4
This assignment is a generalization of the example of the publisher and reviewers, discussed in the
textbook There is a publisher who needs to decide whether to publish a book, and who has the
option of consulting with reviewers. The difference from the example in the book is that, in this
exercise, the publisher may consult with more than one reviewer.
We will assume the same simple model of publishing and of reviewers as in the example in the book.
The outcome of publishing is either Success, with a specified profit, or Fail with a specified loss.
The outcome of not publishing is 0. Consulting with a reviewer costs a specified amount Fee. A
reviewer gives a Boolean answer, Pos or Neg. Reviewers’ opinions are conditionally independent
given the value of Success/Fail. We will assume that all the reviewers are functionally identical;
that is, they all have the same probabilities and they all cost the same. We will also assume that the
publisher initially chooses the number of reviewers they want to consult with (the problem where the
publisher can consult with reviewers sequentially is good deal trickier — if you feel underchallenged
in this course, you might enjoy thinking about it.)
Write the following functions:
Part A
ProbSuccess(ProbSuc,ProbPosSuc,ProbPosFail,NPos,NNeg) where
ProbSuc = P(Success).
ProbPosSuc = P(For | Success).
ProbPosFail = P(For | Fail).
NPos = the number of reviewers who have given a positive review.
NNeg = the number of reviewers who have given a negative review.
This returns the conditional probability of success, given these parameters.
Part B
PubValue1(ValueSuc,ValueFail,Fee,ProbSuc,ProbPosSuc,ProbPosFail,NPos,NNeg) where:
ProbSuc,ProbPosSuc,ProbPosFail,NPos,NNeg are as above.
ValueSuc = Dollar profit if the book succeeds.
CostFail = Cost if the book fails (a positive number).
Fee = Cost of hiring a reviewer (a positive number).
This returns two values: A Boolean indicating whether the publisher should publish and a number
which is the expected return if they follow the advice of ChooseToPub.
Part C
PubValue2(ValueSuc,ValueFail,Fee,ProbSuc,ProbPosSuc,ProbPosFail,N) which returns the expected return if the publisher chooses to consult with N reviewers.
Note: If N reviewers are consulted then P(NPos|Success) and P(NPos|Failure) each follow a
binomial distribution, with p = ProbPosSuccess and p = ProbPosFail respectively. In Matlab, the
function nchoosek(n,k) computes the binomial coefficient C(n, k) = n!/(k!(n − k)!).
1
Part D
OptimalN(ValueSuc,ValueFail,Fee,ProbSuc,ProbPosSuc,ProbPosFail) which returns a pair of
values: The optimal number of reviewers to consult, and the expected profit.
Clearly, the maximum number of reviewers that can possibly be worth consulting is ProbSuc*ValueSuc/Fee;
at that point, you have spent all your expected profit in fees, even if you could completely guarantee a right decision. So just loop through the values of PubValue2 with N going from 0 to
ProbSuc*ValueSuc/Fee, and return the best value and the corresponding N.
2
CSCI 1180 Programming Assignment 5
Consider the following game: You have a collection of coins of weighted coins of m different types.
Specifically, coins of category i come up heads with probability pi
, and you have ki of them. Therefore
the total number of coins C =
Pm
i=1 ki
. You also have an urn, which is initially empty. The coins
are in a pile outside the urn.
You now do the following. If either the urn or the outside is empty, pick a coin at random. If neither
the urn nor the outside is empty, then pick either the urn or the outside with probability 1/2, and
then pick a coin at random from the place that you have chosen. Flip the coin. If it comes up heads,
then put it in the urn; if it comes up tails then put it outside. Keep doing this forever.
This is a Markov process. The state of the system is determined by the number of coins of each
category in the urn. Since there can be between 0 and ki coins of category i in the urn, the total
number of states is q = (1 + k1) · (1 + k2) · . . .(1 + km)
For example, suppose that there is 1 coin of category 1, which is a fair coin, and 2 coins of category
2, which have a 0.8 chance of coming up heads. Then there are 6 states:
State number 1 2 3 4 5 6
# of type 1 in urn 0 0 0 1 1 1
# of type 2 in urn 0 1 2 0 1 2
The probability of a transition from state 2 to state 3 (for example) is calculated as follows: To
go from 2 to 3, you must pick “outside” (probability 1/2); then out of the two coins outside, pick
the one of type 2 (probability 1/2); then flip heads (probability 4/5)); so an overall probability of
(1/2) · (1/2) · (4/5) = 1/5. Computing the other probabilities in a similar way, we end up with the
transition matrix shown below and illustrated in the figure.
The transition matrix is:
T =
0.3 0.1 0 0.25 0 0
0.533 0.575 0.1 0 0.125 0
0 0.2 0.65 0 0 0.167
0.167 0 0 0.35 0.05 0
0. 0.125 0 0.4 0.425 0.133
0 0 0.25 0 0.4 0.7
1
5
0.4
0.42
0.05
0.35
2 3
0.2
0.1
0.65
0.7
6
1
4
0.3 0.57
0.17 0.25 0.13 0.13
0.4
0.13
0.25 0.17
0.53
0.1
Write the following functions:
• States(K) takes a vector K of the number of coins in each category, and returns an m×q array
S where S[I,J] is the number of coins of type I in the urn in state J. In the above example,
the call States([1,2]) should return the array
0 0 0 1 1 1
0 1 2 0 1 2
• Transition(S,H) takes as argument a state array S as output by States and a vector H,
where H[I] is the probability that coins of type I come up heads. It returns the associated
transition matrix between the states. In this example, the call Transition([0,0,0,1,1,1;
0,1,2,0,1,2], [0.5, 0.8]) should return the above transition matrix.
• Stationary(T) computes the stationary distribution for matrix T. This can be computed as
described in section 10.1.1.
• ProbFullUrn(T,N) takes as argument the matrix T and a number N and returns a vector P(I)
of size N of the probabilities that the urn will become full on or before the Ith flip, for I =
1:N. This is easily computed as follows: Modify the transition matrix T so that the final state
(where all the coins are in the urn) can transition only to itself; call this new transition matrix
U. Let S0, the initial probability distribution, be the column vector [1;0; …0], since initially
all the coins are outside. Then P[I] = (U
I
· S0)[Q]). (Q here is q, the number of states).
In the above example ProbFullUrn(T,6) should return [0,0,0.08, 0.184, 0.2888, 0.3859]
• ExpTimeToFill(T,EPS) returns the expected number of flips needed to fill the urn, to accuracy
approximately EPS. For any I, note that the probability that it takes I steps to fill the urn is
equal to (the probability that the urn is full on the Ith step) minus (the probability that it is
full on the I − 1st step). Compute the running sum of the product I times (the probability
that the urn becomes full on the Ith step) until (1−(the probability that the urn is full on the
Ith step) times I is less than EPS; then return the running sum.
For the above transition matrix you get
ExpTimeToFill(T,0.1) = 9.5418
ExpTimeToFill(T,0.01) = 9.6371
ExpTimeToFill(T,0.001) = 9.6468
2
