Today is a LAB day.
We have a LAB #8 to be done during class time, for credit. When you have completed the lab, please let Detelina know so that she can approve and record your work.
LAB#8:
Task 1: to begin our program, we are going to need to pick a
random point (x,y) inside the unit square. Begin by creating a main program.
Now write a function some_point() which calls rand() to get a random
integer, divides it by RAND_MAX using float arithmetic, and returns the
result as a real number between 0 and 1. Remember that using RAND_MAX
means you have to include
Now have your main program call some_point() twice, setting x and y, and print them:
x = some_point ( );
y = some_point ( );
PRINT X AND Y HERE
These should be "random" real numbers between 0 and 1.
Task 2: we are going to need to flip a three-sided coin, that is, to randomly get values 0, 1 or 2. Write a function called flip() which calls rand() to get a random integer, takes its remainder modulo 3, and returns that integer as the "result" of the coin flip. You will have to "declare" int flip(void); at the beginning of your program, also.
Just for a check, call flip() 10 times:
x = some_point ( );
y = some_point ( );
PRINT X AND Y HERE
FOR 10 TIMES
j = flip ( );
PRINT J;
END LOOP
Task 3: Now we imagine that our point (x,y) is inside the triangle whose three corners are (0,0), (1,0) and (0,1). Now we are going to flip the coin and move halfway between (x,y) and corner 0, 1, or 2:
x = some_point ( );
y = some_point ( );
PRINT X AND Y HERE
FOR 10 TIMES
j = flip ( );
PRINT j
if J is 0
x = average of x and 0
y = average of y and 0
but if J is 1
x = average of x and 1
y = average of y and 0
otherwise
x = average of x and 0
y = average of y and 1
PRINT X AND Y HERE
END LOOP
Task 4: Your list of points should stay in the unit square. We'd like to simplify this data by imagining that it occurs on a sort of checkerboard made of "cells". Instead of keeping track of the exact (x,y), we'd just like to keep track of what cell we are in. If we divide the interval [0,1] into n subintervals, than a typical value x corresponds to interval i=(int)(x*n). A point (x,y) will belong in cell (i,j), where i=(int)(x*n) and j=(int)(y*n).
Create a function locate(n,x) which takes as input the number n, and the value x, and returns the interval locate=(int)(x*n). Set the value of n to 10. Use this function to determine the i1 and j1 indices of the cell for each point (x,y).
n = 10; <-- initialize n up here!
x = some_point ( );
y = some_point ( );
PRINT X AND Y HERE
FOR 10 TIMES
j = flip ( );
PRINT j
if J is 0
x = average of x and 0
y = average of y and 0
but if J is 1
x = average of x and 1
y = average of y and 0
otherwise
x = average of x and 0
y = average of y and 1
PRINT X AND Y HERE
i1 = locate(n,x);
j1 = locate(n,y);
PRINT I1 and J1 HERE
END LOOP
Your values I1 and J1 should be between 0 and N-1 = 9.
Task 5: Imagine that instead of looping 10 times, we loop many times. What happens to our point (x,y) as it wanders around? One way to find out is to keep track of how many times it lands in each cell. We can do that by creating a two dimensional array cell[n][n], initializing it to zero, and then setting cell[i1][j1] to 1 if (x,y) lands in the (I1,J1) cell.
n = 10;
int cell[n][n]; <-- remember, we can declare an array this way in C.
SET CELL to 0 <-- Initialize
x = some_point ( );
y = some_point ( );
PRINT X AND Y HERE
FOR 100 TIMES <-- Bump this up to 100
j = flip ( );
PRINT j
if J is 0
x = average of x and 0
y = average of y and 0
but if J is 1
x = average of x and 1
y = average of y and 0
otherwise
x = average of x and 0
y = average of y and 1
PRINT X AND Y HERE
i1 = locate(n,x);
j1 = locate(n,y);
// PRINT I1 and J1 HERE <-- delete or comment this out now.
SET CELL(I1,J1) TO 1 <-- How do you do this?
END LOOP
Task 6: Our array CELL is an N by N table of integers, containing the number of times our "bouncing" value (x,y) ended up in each cell. Rather than printing the set of values, can we try to visualize it, that is, make a picture that will reveal any patterns?
Copy the function bw.c. Its declaration is
void bw ( int m, int n, int a[m][n] );
We will need to use this function in your program, so add this
declaration to the list of function declarations at the beginning of your program.
Increase the value of n to 500. Increase the number of loop iterations to 100,000. Once the loop has been completed, call bw(n,n,cell).
n = 500; <-- Increase
int cell[n][n];
Set CELL to 0
x = some_point ( );
y = some_point ( );
PRINT X AND Y HERE
FOR 100000 TIMES <-- Bump this up to 100,000
j = flip ( );
PRINT j
if J is 0
x = average of x and 0
y = average of y and 0
but if J is 1
x = average of x and 1
y = average of y and 0
otherwise
x = average of x and 0
y = average of y and 1
PRINT X AND Y HERE
i1 = locate(n,x);
j1 = locate(n,y);
SET CELL(I1,J1) to 1
END LOOP
call "bw()" here!
When you call "bw", it creates a file called "picture.pbm" containing graphics information. To see what's in there, type
eog picture.pbm
eog is simply one of many programs that can display graphics files.
If things went correctly, you've drawn a picture of the Sierpinski triangle.
HOMEWORK #8 (must be turned in by next Thursday):
Homework 8.1: the combinatorial coefficient C(n,k) counts the number of ways you can choose k objects, given a choice of n. For instance, if you can order 3 free toppings from a list of 5, the value of C(5,3) is 10, meaning your list could be one of ABC, ABD, ABE, ACD, ACE, ADE, BCD, BCE, BDE, or CDE. There is a formula:
C(n,k) = n! / k! / (n-k)!
(where n! = n * (n-1) * ... * 2 * 1 is the "factorial function")
so that C(5,3) = 5!/3!/2! = 120 / 6 / 2 = 10.
Write a function combin(n,k) which computes the combinatorial coefficient for any n and k. (N and K will always be nonnegative, and K will always be no greater than N.) You should write a function to evaluate the factorial, and combin() should use that function.
Demonstrate that your function is working by computing C(5,3), C(8,4) and C(9,6).
Homework 8.2: the transpose of an M by N matrix A is an N by M matrix B for which B(i,j) is set to A(j,i). Write a function "tranpose()" which takes a 4x5 matrix A and a 5x4 matrix B as input. The function should store into B the transpose of the matrix A. Test your function by using as input the matrix A(i,j) = 10*i+j:
0 1 2 3 4
10 11 12 13 14
20 21 22 23 24
30 31 32 33 34
Your main program should print out the tranpose matrix B after it has
been computed.
Homework 8.3: (Graduate students only!). The program "midpoint_example.c" uses the midpoint method, with n intervals, to estimate the integral of a function from a to b, by calling the function double midpoint(int n, double a, double b). Right now, midpoint() assumes that the function to be integrated is named "poly()". But I want to modify the program so that I can call midpoint once to integrate "poly()" and another time to integrate a function called "wave".
You can make this happen by adding one extra argument to midpoint(), which "declares" the function, (explains how it is used) and gives it a temporary or local name (which might as well be "f").
How can we pass the function f as an argument? If it was a real number, we'd just say "double f". But f is slightly more complicated, it's a function that takes a real input, and returns a real output. In other words, instead of a "double f", it's a "double f ( double x )". That is what your extra argument looks like. You make this change in TWO places: where we declare midpoint() at the very top of the program, and where the midpoint function begins.
Once you've got the extra input argument, modify the line "q = q + dx * ?(x)" by replacing "?" by "f", the local name I assume you chose. Now fix the two calls to midpoint, passing in "poly" the first time and "wave" the second time. Your integral estimates should be about 27 the first time and 3.7 the second.
If you are having trouble with this problem, look at the "euler.c" program from Wednesday.
For each of your homework programs, turn in a copy of the program, and a copy of the output.
Programs we might discuss: